This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

The diffraction integral for a circular aperture is [✶]

ℰ({x,y,z};a)=(1/(2π))(∫_0^(2π))(∫_0^a)ρ exp(i((ρ^3)cos(3α)-z (ρ^2)-x ρ cos(α)-y ρ sin(α)))dρ dφ

Make an animation of the curves Im(ℰ({x,y,z};a))=0 as a function of a for a small z for (x^2)+(y^2)≲10.

ℐN[{x_, y_, z_}, a_] := 
NIntegrate[ρ Exp[I (ρ^3 Cos[3 α] - z ρ^2 - x ρ Cos[α] - y ρ Sin[α])],
           {α, 0, Pi, 2Pi}, {ρ, 0, a}]

ℐN[{1, 2, 1}, 2] // Timing

Because we are interested in the values of ℰ({x,y,z};a) as a function of a, we use NDSolve with respect to a instead. This gives as values of ℰ({x,y,z};a) for all values of a at once and we only have to carry out a 1D numerical integration. We calculate the values of ℰ({x,y,z};a) on a x,y-grid in a polar coordinate system and store interpolating functions ipo. We make use of the rotational symmetry of ℰ({x,y,z};a) to calculate the values of ℰ({x,y,z};a) only in a sector.

Off[NIntegrate::"ploss"];
(* 1D integrand *)
𝒾𝓃𝓉[{x_, y_, z_}, ρ_?NumberQ] := 
NIntegrate[ρ Exp[I (ρ^3 Cos[3 α] - z ρ^2 - x ρ Cos[α] - y ρ Sin[α])],
           {α, 0, 2Pi}, PrecisionGoal -> 7,
           (* periodic integrand *) Method -> Trapezoidal]

(* parameters *)
R = 12; z = 1/10; ppφ = 16; ppr = 48; 

(* use polar coordinates to easily make use of rotational symmetry;
   x == r Cos[φ], y == r Sin[φ] *)
Do[Do[ipo[{r, φ, z}] = (* solve ODE with respect to a == ρ *)
      NDSolve[{ℰ'[ρ] == 𝒾𝓃𝓉[{r Cos[φ], r Sin[φ], z}, ρ], ℰ[0] == 0},
              ℰ, {ρ, 0, 2}, PrecisionGoal -> 4][[1, 1, 2]],
      {φ, 0, Pi/3, Pi/3/ppφ}], 
   {r, 0, R, R/ppr}]

For a given value of a, the function imZeroLinesPlot generates a plot of the lines Im(ℰ({x,y,z};a))=0.

imZeroLinesPlot[a_] :=
Module[{data, lcp, zeroImLines, (* rotation matrix *)
        ℛ = {{-1, Sqrt[3]}, {-Sqrt[3], -1}}/2.}, 
 (* data in polar coordinates *)
 data = Table[ipo[{r, φ, z}][a],
              {φ, 0, Pi/3, 2Pi/3/ppφ}, {r, 0, R, R/ppr}];
 (* contour plot of vanishing imaginary part *)
 lcp = ListContourPlot[Im[data], PlotRange -> All,
           Axes -> True, ContourShading -> False,
           DisplayFunction -> Identity,  Contours -> {0},
           MeshRange -> {{0, R}, {0, Pi/3}}];
 (* the actual contour lines *)
 zeroImLines = Cases[Graphics[lcp], _Line, Infinity];
 If[zeroImLines =!= {}, 
    (* final graphic *)
    Graphics[{{Thickness[0.002], Hue[0.78 a/2],
     (* generate the five equivalent parts *)
     {Map[ℛ.ℛ.#&, #, {-2}], Map[ℛ.#&, #, {-2}], #}&[
     {Map[{1, -1} #&, #, {-2}], #}&[
      Apply[Function[{r, φ}, {r Cos[φ], r Sin[φ]}], 
            zeroImLines, {-2}]]]},
    (* outer circle *)
    {Thickness[0.01], GrayLevel[0.5], Circle[{0, 0}, R]}},
    Frame -> True, PlotRange -> All, AspectRatio -> Automatic], 
     {}]]

Here are the curves Im(ℰ({x,y,z};a))=0 for various a shown in different colors.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
{(* 2D graphic *)
 Show[Table[imZeroLinesPlot[a], {a, 2, 0, -2/12}]],
 (* 3D graphic *)
 Show[Graphics3D[
 Table[If[# === {}, {}, #[[1, 1]] /. {t_, Hue[x_], lines__} :>
           {t, Hue[x], Map[Append[#, 2x/0.78]&, {lines}, {-2}]}]&[
      imZeroLinesPlot[a]], 
      {a, 2, 0, -2/18}]], BoxRatios -> {1, 1, 3/2}, Axes -> True,
     AxesLabel -> {"x", "y", "a"}, ViewPoint -> {-1, -2, -3/4}]}]]]

And here is the corresponding animation. For an interpretation of the changing curve shape, see [✶].

Do[If[# =!= {}, Show[#, Frame -> False]]&[imZeroLinesPlot[a]],
   {a, 0, 2, 2/48}]