Diffraction Limit.sagews

x,y = var('x','y')
#bessel_J(1,x)
#plot(bessel_J(0,x), (x,0,6*pi))
@interact
def _( undouble_offset=(0,2*pi) ):
#    offset = 1.22*pi / 2
    offset = undouble_offset / 2
    total_plot = plot((bessel_J(0,x-offset)+bessel_J(0,x+offset))^2, (x,-4*pi,4*pi), figsize=[10,10], color='red')
    total_plot += density_plot(log(1 + (bessel_J(0,sqrt((x-offset)^2+y^2)) + bessel_J(0,sqrt((x+offset)^2+y^2)))^2), (x,-4*pi, 4*pi), (y,-4*pi, 4*pi), figsize=[10,10], dpi=200 )
    show( total_plot )

answer = (derivative(bessel_J(0,x),x)==0).find_root(0.1,2*pi)
show( N(answer/pi,256) ) # Solve up to 256 bits

Interact: please open in CoCalc

$\displaystyle 1.219669891266504503529861345818056010952193017369386984539883844012879946397$