#!/usr/bin/env Rscript
# Diffraction pattern for two airy disks, illustrating the Rayleigh criterion
# Released to Public Domain by Spencer Bliven
x = seq(-5,5,.01)
y = seq(-5,5,.01)

airy <- function(xx,yy,x1,y1) {
    r<-sqrt((xx-x1)**2+(yy-y1)**2)*1.22*pi
    z = (2*besselJ(r,1)/r)**2
    z
}
diffractionPlot = function(spacing=1,useRaster=F) {
    z1 = outer(x, y, FUN=function(xx,yy){airy(xx,yy,-spacing/2,0)})
    z2 = outer(x, y, FUN=function(xx,yy){airy(xx,yy,spacing/2,0)})
    z1[is.na(z1)]=1.0
    z2[is.na(z2)]=1.0
    z = z1+z2
    grad = grey((0:(2**16)/2**16))
    par(pty="s")
    # Non-linear z scale to increase contrast
    image(x,y,z**(1/4),zlim=c(0,1),col=grad,axes=F,useRaster=useRaster)
}

png("airy_0.5.png",bg="#FFFFFF00")
diffractionPlot(.5)
dev.off()
png("airy_1.0.png",bg="#FFFFFF00")
diffractionPlot(1)
dev.off()
png("airy_2.0.png",bg="#FFFFFF00")
diffractionPlot(2)
dev.off()