Circular Aperture Diffraction Pattern

Question

I am unsure if this is even possible in Mathematica, but my optics professor assigned a project to me on circular apertures and Fourier transforms. I have found plenty on the Airy function, but I am wondering if there is any way to produce an image like the one below in Mathematica:

Answer 1

Physicist chiming in - Hi!.

I believe there has been some confusion here. It seems to me that OP is meaning to plot an Airy disk which was studied by G.B. Airy but is not given by the Airy function. It is given by the Fourier transform of the indicator function of the unit circle, which actually happens to be a Bessel function (see e.g. wikipedia).

If I understood this right, then the correct solution is as follows:

DensityPlot[BesselJ[1, Sqrt[x^2 + y^2]]/Sqrt[
            x^2 + y^2], {x, -60, 60}, {y, -60, 60}, PlotPoints -> 100, 
            PlotRange -> All]

You can play around with the options of DensityPlot to increase the contrast, add a legend, or change the colour scheme into something more similar to your intended image. I leave this to you.

--

For the mathematically inclined: we are dealing with what we physicists call Fraunhofer diffraction. Given a profile $f(x_1,x_2)$, the corresponding diffraction pattern is proportional to $\tilde f(\xi_1,\xi_2)$. In our case, the profile is a solid disk, so $f(x_1,x_2)=1_{D(1)}=\theta(1-\sqrt{x_1^2+x_2^2})$, with $\theta$ the step-function.

Unfortunately,

FourierTransform[HeavisideTheta[1 - Sqrt[x^2 + y^2]], {x, y}, {a, b}]

is returned unevaluated. Fortunately, this Fourier integral is a well-known result, the quoted Bessel function. It would be nice to know why MMA is unable to evaluate this classic integral. Oh well.

Answer 2

I feel a bit bad about sliding in so lately. @Nasser's answer was already very good. However, the function is very oscillatory, so I thought it might be a better idea to use a Texture to render it. I use ParametricPlot for that along with a suitable TextureCoordinateFunction. Since the function is radially symmetric, the texture needs only be produced for a very narrow strip. This is done by the function tex.

R = 20;(*half box size*)
n = 1000;(*resolution of texture*)
x0 = N[Range[0, R, R/n]];
tex = cut \[Function] Module[{c, cols}, 
    c = Map[ 
      r \[Function] If[Abs[r] <= 10^-12, 0., Max[cut, Log10[Abs[2.0 BesselJ[1, Pi r]/(Pi r)]]]], 
      Sqrt[2.] x0
      ];
    cols = ColorData["SunsetColors"] /@ Rescale[c];
    Image[ConstantArray[Developer`ToPackedArray[List @@@ (cols)], 10]]];
texfun = {x, y, u, v} \[Function] Evaluate[{Sqrt[(u^2 + v^2)]/R/Sqrt[2], 0.5}];
plot = cut \[Function] ParametricPlot[{x, y}, {x, -R, R}, {y, -R, R},
    TextureCoordinateFunction -> texfun, 
    TextureCoordinateScaling -> False,
    PlotStyle -> {Opacity[1], Texture[tex[cut]]},
    BoundaryStyle -> None,
    PlotPoints -> 100,
    PlotRangePadding -> 0, Axes -> False];

Moreover, tex produces a logarithmic texture where all values below the truncation argument cut are replaced by cut. The effect of different truncatations can be seen below:

g = GraphicsGrid[Partition[
   Table[
    plot[cut],
    {cut, {-1., -1.5, -2., -2.5, -3., -3.5}}
    ],
   3], ImageSize -> Full]

Answer 3

The Fourier transforms of Airy function is given by what is called Sombrero function

So here is something to get you started. This finds the radius to each pixel over some square region and applies the definition of Sombrero and plots the result using ListDensityPlot

meshgrid[x_List,y_List]:={ConstantArray[x,Length[x]],Transpose@ConstantArray[y,Length[y]]};

n       = 50; (* grid size *)
z       = 1/2; (*adjust as needed*)
x0      = Range[1,n,z]-n/2; (*shift center*)
{xx,yy} = meshgrid[x0,x0]; 
r       = Sqrt[xx^2+yy^2];  (*radius*)
sombro  = Map[ If[#==0,0,2.0 BesselJ[1,Pi #]/(Pi #)]&,r,{2}];

ListDensityPlot[sombro,ColorFunction->"SunsetColors"]

You can play with different plots or option to make it closer to what you showed. For example, using z=1/4 or different n or other scaling one can do.

Answer 4

I can't improve on the great plots already presented, but MMA can indeed solve the problem. Start with an energy wave from a circular aperture

Eeq=dE==A dx dy Sin[2Pi f t-(2\[Pi] r)/lambda]

Point P on the observation screen has coordinates X, Y, Z, with the origin at the center of the aperture. Z is the distance from the center of the aperture to the center of the screen. Distance from center of aperture to point P on screen r0 is

r0=Sqrt[X^2+Y^2+Z^2]

Dist r, from a point x,y,0 on the aperture to Point P on the screen is

r=Sqrt[(X-x)^2+(Y-y)^2+Z^2]

Use the small aperture approximation to make simplifications

rhs=(r^2-r0^2//ExpandAll)/.{x^2->0,y^2->0}

Clear[r,r0]

r^2-r0^2//Factor
(*(r-r0) (r+r0)*)

lhs=(%/.r+r0->2r0)

Solve[lhs==rhs,r]//Flatten

Rrule=r->r0-(x X+y Y)/r0

Eeq=Eeq/.Rrule 
 (* dE==A dx dy Sin[2 Pi f t-(2 Pi (r0-(x X+y Y)/r0))/lambda] *)

Convert from x, y on the aperture to rho, phi and X, Y on the screen to R Phi.

x=rho Cos[phi];
y=rho Sin[phi];
X=R Cos[Phi];
Y=R Sin[Phi];

Eeq
(* dE==A dx dy Sin[2 \[Pi] f t-(2 \[Pi] (r0-(R rho Sin[phi] Sin[Phi]+R rho Cos[phi] Cos[Phi])/r0))/lambda] *)

Since the diffraction pattern is symmetric about Phi, simplify by setting Phi to 0.

Eeq=Eeq/.Phi->0//ExpandAll   

(*-A dx dy Sin[-(2 Pi f t)-(2 Pi R rho Cos[phi])/(lambda r0)+(2 Pi r0)/lambda]*) 

Eeq=A dx dy Sin[2 Pi f t+(2 Pi R rho Cos[phi])/(lambda r0)-(2 Pi r0)/lambda]//.{-2 f Pi t+(2 Pi r0)/lambda->alpha,-((2 Pi R rho  Cos[phi])/(lambda r0))->beta}//TrigExpand;

Eeq=%/.{alpha->-2 f Pi t+(2 Pi r0)/lambda,beta->-((2 Pi R rho  Cos[phi])/(lambda r0))}

Integrate over the aperture of radius a in cylindrical coordinates, dx dy -> rho drho dphi.

Eeq=%/.dx dy->rho drho dphi

(*A dphi drho rho Cos[(2 Pi r0)/lambda-2 \[Pi] f t] Sin[(2 Pi R rho Cos[phi])/(lambda r0)]-A dphi drho rho Sin[(2 Pi r0)/lambda-2 Pi f t] Cos[(2 Pi R rho Cos[phi])/(lambda r0)]*)

$Assumptions=a>0&&R>0&&lambda>0&&r0>0

Works in the following form:

En=A  Cos[(2 Pi r0)/lambda-2 Pi f t]Integrate[rho Sin[(2 Pi R rho Cos[phi])/(lambda r0)],{phi,0,2 Pi},{rho,0,a}]-A Sin[(2 Pi r0)/lambda-2 Pi f t] Integrate[rho Cos[(2 Pi R rho Cos[phi])/(lambda r0)],{phi,0,2 Pi},{rho,0,a}]

(*-((a A lambda r0 BesselJ[1,(2 a Pi R)/(lambda r0)] Sin[(2 Pi r0)/lambda-2 Pi f t])/R)*)

Remember that R is the radial distance from the center of the screen to the observation pt. Get to the familiar form by

En/.R->(lambda r0 xx)/(2 a Pi)

(*-((2 Pi a^2 A BesselJ[1,xx] Sin[(2 Pi r0)/lambda-2 Pi f t])/xx)*)