Nerd Rage

In this post on the prin­ci­ples of scalar dif­frac­tion there were quite a num­ber of inte­gral expres­sions for dif­frac­tion of light, start­ing from the Huygens-Fresnel prin­ci­ple and apply­ing var­i­ous degrees of approx­i­ma­tion. Inte­gral for­mu­la­tions are usu­ally very abstract and not very illus­tra­tive. How­ever, with today’s avail­able com­puter power we can turn these inte­grals into more descrip­tive imagery in no time (and then waste a few more hours turn­ing the raw data into nice view­graphs).

To keep it sim­ple we’ll look at only two dif­fer­ent aper­tures, one cir­cu­lar and one square, illu­mi­nated by a plain plane wave. We will com­pare the exact solu­tion (using Rayleigh-Sommerfeld 1 coef­fi­cients), Fresnel’s approx­i­ma­tion, and Fraunhofer’s approx­i­ma­tion at var­i­ous dis­tances from the aper­ture, start­ing very close to the aper­ture (at a dis­tance of $10\lambda$ from a cir­cu­lar aper­ture of radius $10\lambda$). The inten­sity dif­frac­tion pat­terns obtained with the dif­fer­ent meth­ods are these:

As expected, both approx­i­ma­tions fail close to the aper­ture, the Fraun­hofer approx­i­ma­tion even fails quite epicly. Fres­nel fails because $z$ is not much larger than the extent of the aper­ture and/or screen, and Fraun­hofer requires that the aper­ture be very, very small com­pared to the dis­tance $z$ to the screen. Let’s see what hap­pens when $z$ is increased ten­fold:

Here, Fresnel’s approx­i­ma­tion comes very close, but Fraun­hofer is still amiss. Mov­ing away a bit fur­ther, we have:

When $z$ becomes large enough, the Fraun­hofer inte­gral (and thus the Fourier trans­form) indeed describe the true field quite well.

The results are sim­i­lar when look­ing at a square aper­ture:

Again, Fresnel’s approx­i­ma­tion works pretty well at $z=100\lambda$ and Fraunhofer’s does from $1000\lambda$ onwards.

Com­pu­ta­tion

The numer­i­cal com­pu­ta­tions for these images were done with MATLAB. Inter­est­ingly, the approx­i­mate inte­grals didn’t con­sis­tently com­pute faster than the exact one – at least not while numer­i­cally cal­cu­lat­ing the inte­grals directly (instead of using e.g. the Fourier trans­form for the Fraun­hofer for­mu­la­tion). This may be due to my very inel­e­gant and rather brutish-force approach. The code was writ­ten within min­utes and con­sists of two files, the inte­grand func­tion and the inte­gra­tion wrap­per:

function out = integrand_exact(px, py, qx, qy, z, k)

qminusp = sqrt((qx-px).^2 + (qy-py).^2 + z.^2);
kappa = z ./ qminusp / 1i * k / (2*pi); % cos alpha / i / lambda

out = 1 .* kappa ./ qminusp .* exp(-1i*k*qminusp);

end

We don’t worry about the aper­ture bounds here since we can just as well define the inte­gra­tion bound­aries, see below. If the aper­ture were not as sim­ply shaped we’d define it here in the inte­grand func­tion. Imple­ment­ing the Fres­nel and/or Fraun­hofer inte­grands is pretty straight­for­ward. The inte­gra­tion wrap­per looks like this:

k = 2*pi;
z = 1000;

qx = linspace(-200, 200, 501);
qy = linspace(-200, 200, 501);

for ii=1:length(qx)
    for jj=1:length(qy)

        A(ii,jj) = TwoD(@(x, y) (integrand_exact(x, y, qx(ii), qy(jj), z, k)),...
            0, 2*pi, 0, 10, 'Sector', true, 'AbsTol', 1e-5);

    end
end

The vari­ables qx and qy con­tain the posi­tions in the plane at z at which we want to deter­mine the field A.

Instead of using the ol’ dblquad (which can be a bit slow) we use Lawrence F. Shampine’s TwoD two-dimensional quad­ra­ture code. One advan­tage of TwoD is its abil­ity to deal with cir­cu­lar inte­gra­tion regions, as shown above. Plus, it’s much faster (in most cases more than an order of mag­ni­tude). As with the MATLAB-internal quad­ra­ture func­tions, the inte­grand is passed via an anony­mous func­tion han­dle that is cre­ated on-the-fly – see [1].

For the sim­ple geome­tries above, it is suf­fi­cient to com­pute only a cross sec­tion of the field (leav­ing out e.g. the inter­a­tion over qy). Due to the sym­me­try of the prob­lem we can cal­cu­late the com­plete field on the screen by

rho = repmat(qx, [length(qy), 1]).^2 + repmat(qy.', [1, length(qx)]).^2;
rho = sqrt(rho);
B = interp1(qx, A, rho, 'pchip');

for the cir­cu­lar aper­ture and

B = repmat(A.', [length(A), 1]) .* repmat(A, [1, length(A)]);

for the square aper­ture. The cir­cu­lar aper­ture ver­sion cal­cu­lates the radius $\rho$ from the ori­gin for each point of the screen and then inter­po­lates the cross sec­tion data for each $\rho$. The square ver­sion uses the inde­pen­dence of both Carte­sian com­po­nents. This reduces (in the above case) com­pu­ta­tion time to 1/500.

---

[1] MATLAB anony­mous func­tions [Nerd Rage]

---

last posts in diffraction:

• Paraxiality
• Scalar Diffraction - Fourier
• Scalar Diffraction - Huygens, Fresnel, Fraunhofer