# Circular Aperture Diffraction Pattern

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:

• Also check out the Wolfram Demonstration here. This gives a variety of examples of how you can use 2D Fourier Transforms to calculate diffraction patterns with different shapes, and sizes of apertures. You can also download the notebook to see how it is coded. – Dunlop Nov 29 '17 at 7:59

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.

• @Nasser I just would like to point out that the scale of the OP's plot is logarithmic. I am not sure, but maybe, we need the plot of Log[Abs[BesselJ[1, Sqrt[x^2 + y^2]]/Sqrt[x^2 + y^2]]]? This might also increase the contrast if neagtive values get truncated somewhere... – Henrik Schumacher Nov 29 '17 at 10:19
• @HenrikSchumacher I think it's better to plot Log[Abs[...]+1] so that $0$ is mapped to $0$. This produces a nicer plot than that without the +1 – AccidentalFourierTransform Nov 29 '17 at 10:27
• Appropriate username... – ciao Nov 29 '17 at 12:23
• That FourierTransform returns unevaluated, but you can always do it in spherical coordinates: Integrate[r Exp[I k r Cos[f]], {r, 0, 1}, {f, 0, 2 \[Pi]}, Assumptions -> k > 0]. – march Nov 29 '17 at 18:21
• Please, revise your notes in Optics, and for the God sake! As a physicist, do not reference Wikipedia :))... – José Antonio Díaz Navas Apr 30 '18 at 18:25

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[DeveloperToPackedArray[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]


• Here's another one abusing symmetry:ParametricPlot[r {Cos[f], Sin[f]}, {r, .01, 5}, {f, 0, 2 Pi}, ColorFunction -> Function[{x, y, r, f}, ColorData["SunsetColors"][2 Abs[BesselJ[1, Pi r]/(Pi r)]]] , ColorFunctionScaling -> False , Axes -> False , BoundaryStyle -> None , Epilog -> {Thick, Circle[{0, 0}, 5]} ] – Kuba Nov 29 '17 at 13:37

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];
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.

• The only sombrero (hat) function I know is the Mexican hat, i.e., the Laplacian of a 2D Gaussian function. I never heard about this. :)) – José Antonio Díaz Navas Nov 29 '17 at 13:40
• @JoséAntonioDíazNavas Nah, the one and only Mexican hat is $f(x,y,z) = a r^4 - b r^2$. Also known as Landau potential... ;o) – Henrik Schumacher Nov 29 '17 at 18:31
• hats off to potential. – Sumit Nov 30 '17 at 8:10

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)*)
`