In this post on the principles of scalar diffraction there were quite a number of integral expressions for diffraction of light, starting from the Huygens-Fresnel principle and applying various degrees of approximation. Integral formulations are usually very abstract and not very illustrative. However, with today’s available computer power we can turn these integrals into more descriptive imagery in no time (and then waste a few more hours turning the raw data into nice viewgraphs).

To keep it simple we’ll look at only two different apertures, one circular and one square, illuminated by a plain plane wave. We will compare the exact solution (using Rayleigh-Sommerfeld 1 coefficients), Fresnel’s approximation, and Fraunhofer’s approximation at various distances from the aperture, starting very close to the aperture (at a distance of $10\lambda$ from a circular aperture of radius $10\lambda$). The intensity diffraction patterns obtained with the different methods are these:

As expected, both approximations fail close to the aperture, the Fraunhofer approximation even fails quite epicly. Fresnel fails because $z$ is not much larger than the extent of the aperture and/or screen, and Fraunhofer requires that the aperture be very, very small compared to the distance $z$ to the screen. Let’s see what happens when $z$ is increased tenfold:

Here, Fresnel’s approximation comes very close, but Fraunhofer is still amiss. Moving away a bit further, we have:

When $z$ becomes large enough, the Fraunhofer integral (and thus the Fourier transform) indeed describe the true field quite well.

The results are similar when looking at a square aperture:

Again, Fresnel’s approximation works pretty well at $z=100\lambda$ and Fraunhofer’s does from $1000\lambda$ onwards.

#### Computation

The numerical computations for these images were done with MATLAB. Interestingly, the approximate integrals didn’t consistently compute faster than the exact one – at least not while numerically calculating the integrals directly (instead of using e.g. the Fourier transform for the Fraunhofer formulation). This may be due to my very inelegant and rather brutish-force approach. The code was written within minutes and consists of two files, the integrand function and the integration wrapper:

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 aperture bounds here since we can just as well define the integration boundaries, see below. If the aperture were not as simply shaped we’d define it here in the integrand function. Implementing the Fresnel and/or Fraunhofer integrands is pretty straightforward. The integration wrapper 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 variables `qx`

and `qy`

contain the positions in the plane at `z`

at which we want to determine the field `A`

.

Instead of using the ol’ `dblquad`

(which can be a bit slow) we use Lawrence F. Shampine’s TwoD two-dimensional quadrature code. One advantage of TwoD is its ability to deal with circular integration regions, as shown above. Plus, it’s much faster (in most cases more than an order of magnitude). As with the MATLAB-internal quadrature functions, the integrand is passed via an anonymous function handle that is created on-the-fly – see [1].

For the simple geometries above, it is sufficient to compute only a cross section of the field (leaving out e.g. the interation over `qy`

). Due to the symmetry of the problem we can calculate the complete 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 circular aperture and

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

for the square aperture. The circular aperture version calculates the radius $\rho$ from the origin for each point of the screen and then interpolates the cross section data for each $\rho$. The square version uses the independence of both Cartesian components. This reduces (in the above case) computation time to 1/500.

[1] MATLAB anonymous functions [Nerd Rage]

last posts in diffraction:

I really appreciate the time and energy you put into these posts to make us readers a little brighter each time we read one of those posts!

I don’t know how to get the Irradiance Image. I did:

» imagesc(qx,qy,rho)

However, I got a different result.