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:

enter image description here

  • 2
    $\begingroup$ 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. $\endgroup$
    – Dunlop
    Commented Nov 29, 2017 at 7:59

4 Answers 4


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]

enter image description here

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.


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.

  • 3
    $\begingroup$ @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... $\endgroup$ Commented Nov 29, 2017 at 10:19
  • $\begingroup$ @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 $\endgroup$ Commented Nov 29, 2017 at 10:27
  • 4
    $\begingroup$ Appropriate username... $\endgroup$
    – ciao
    Commented Nov 29, 2017 at 12:23
  • 1
    $\begingroup$ 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]. $\endgroup$
    – march
    Commented Nov 29, 2017 at 18:21
  • 1
    $\begingroup$ Please, revise your notes in Optics, and for the God sake! As a physicist, do not reference Wikipedia :))... $\endgroup$ Commented Apr 30, 2018 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[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[
    {cut, {-1., -1.5, -2., -2.5, -3., -3.5}}
   3], ImageSize -> Full]

enter image description here

  • 3
    $\begingroup$ 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]} ] $\endgroup$
    – Kuba
    Commented Nov 29, 2017 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


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


Mathematica graphics

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.

  • $\begingroup$ 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. :)) $\endgroup$ Commented Nov 29, 2017 at 13:40
  • 1
    $\begingroup$ @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) $\endgroup$ Commented Nov 29, 2017 at 18:31
  • 1
    $\begingroup$ hats off to potential. $\endgroup$
    – Sumit
    Commented Nov 30, 2017 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


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


Use the small aperture approximation to make simplifications



(*(r-r0) (r+r0)*)



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

 (* 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];

(* 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.


(*-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)]*)


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.