## Scalar Diffraction - Imagery

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 Huy­gens-Fres­nel prin­ci­ple and apply­ing var­i­ous degrees of approx­i­ma­tion. Inte­gral for­mu­la­tions are usu­al­ly very abstract and not very illus­tra­tive. How­ev­er, with today’s avail­able com­put­er pow­er 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­nat­ed by a plain plane wave. We will com­pare the exact solu­tion (using Rayleigh-Som­mer­feld 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­si­ty dif­frac­tion pat­terns obtained with the dif­fer­ent meth­ods are these:

Fig­ure 1: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 10λ behind a cir­cu­lar aper­ture of 10λ radius which is illu­mi­nat­ed by a plane wave, com­put­ed using the exact solu­tion, Fres­nel and Fraun­hofer approx­i­ma­tions. Red cir­cle shows size of aper­ture for com­par­i­son, bot­tom graphs shows radi­al cross sec­tions.

As expect­ed, 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 larg­er 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:

Fig­ure 2: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 100λ behind the cir­cu­lar aper­ture of Fig. 1.

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:

Fig­ure 3: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 1000λ behind the cir­cu­lar aper­ture of Fig. 1.

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

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

Fig­ure 4: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 10λ behind a square aper­ture of 20λ edge length which is illu­mi­nat­ed by a plane wave, com­put­ed using the exact solu­tion, Fres­nel and Fraun­hofer approx­i­ma­tions. Red square shows size of aper­ture for com­par­i­son, bot­tom graphs shows cross sec­tions.

Fig­ure 5: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 100λ behind the square aper­ture of Fig. 4.

Fig­ure 5: light inten­si­ty dis­tri­b­u­tion in a dis­tance of 1000λ behind the square aper­ture of Fig. 4.

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

#### Computation

The numer­i­cal com­pu­ta­tions for these images were done with MATLAB. Inter­est­ing­ly, the approx­i­mate inte­grals didn’t con­sis­tent­ly com­pute faster than the exact one – at least not while numer­i­cal­ly cal­cu­lat­ing the inte­grals direct­ly (instead of using e.g. the Fouri­er 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 with­in 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 wor­ry 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 pret­ty 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-dimen­sion­al quad­ra­ture code. One advan­tage of TwoD is its abil­i­ty to deal with cir­cu­lar inte­gra­tion regions, as shown above. Plus, it’s much faster (in most cas­es more than an order of mag­ni­tude). As with the MAT­LAB-inter­nal quad­ra­ture func­tions, the inte­grand is passed via an anony­mous func­tion han­dle that is cre­at­ed 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: