# Fresnel diffraction propagation code using angular spectrum method

I'm trying to follow along with the calculation of Fresnel diffraction presented in Siegman's Lasers (1986), page 660. It follows the angular propagation of plane waves in the paraxial approximation. I've gathered that calculating the diffraction pattern formed by an illuminated planar aperture on some other parallel plane is a 3-step process:

1. Given a known field over an aperture, perform an inverse Fourier transform to calculate the angular spectrum of the field at that plane.

(Eq. 89) $$U_{pw}(s_x,s_y,z_0) = \int\int{u(x,y,z_0) \times e^{+j2\pi(s_xx+s_yy)}}dxdy$$

1. Propagate the angular spectrum to the other plane using the paraxial approximation.

(Eq. 90) $$U_{pw}(s_x,s_y,z) = U_{pw}(s_x,s_y,z_0) \times e^{-jk(z-z_0)+j\pi\lambda(s_x^2+s_y^2)(z-z_0)}$$

1. Retrieve the field at the new plane by performing a Fourier transform.

(Eq. 91) $$u(x,y,z) = \int\int{U_{pw}(x,y,z) \times e^{-j2\pi(s_xx+s_yy)}}ds_xds_y$$

I've written some MATLAB code to try and calculate the diffraction pattern given an illuminated circular aperture. I can predict the Airy disc diameter, but I can't get a matching solution by using transforms. I think my issue may be related to how I am taking samples for the transforms?

%% Example: Circular aperture.

lambda    = 632.8e-9; % 632.8nm wavelength light.
diameter  = 100e-6;   % 100um circular aperture diameter.
aperturez = 0;        % Aperture is located at z=0.
imagez    = 100e-2;   % 100cm from the aperture to the image plane.

aperture_sample_width = 5e-4; % Sample the aperture for a width of 0.5mm.
% image_sample_width    = 2e-2; % Sample the image plane for a width of 2cm.

nsamples = 1024; % Sample count (prefer power of 2).

predict_airy_disc_radius = 1.22 * lambda * (imagez - aperturez) / diameter;
fprintf('Predicted Airy disc radius: %0.2f mm\n', ...

figure

% Plot the magnitude of the field at the (circular) aperture plane.
aperture_field = @(u,v) (sqrt((u.^2 + v.^2)) < (diameter/2));
aperture_grid  = linspace(-aperture_sample_width/2, ...
aperture_sample_width/2, nsamples);
[ugrid, vgrid] = meshgrid(aperture_grid);
subplot(2,2,1)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
aperture_field(ugrid,vgrid))
axis square, xlabel 'mm', ylabel 'mm'

% Plot the magnitude of the angular spectrum at the aperture.
aperture_spectrum = ifft2(aperture_field(ugrid,vgrid));
aperture_spectrum_spacing = 1/aperture_sample_width;
aperture_spectrum_grid = ((1:nsamples) - (nsamples/2)) * ...
aperture_spectrum_spacing;
[sx,sy] = meshgrid(aperture_spectrum_grid);
subplot(2,2,2)
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
ifftshift(aperture_spectrum .* conj(aperture_spectrum)))
axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'

% Plot the magnitude of the angular spectrum at the image plane.
image_spectrum = aperture_spectrum .* ...
exp((-1j * (2*pi/lambda) * (imagez - aperturez)) + ...
( 1j * pi*lambda * (sx.^2 + sy.^2) * (imagez - aperturez)));
subplot(2,2,4)
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
ifftshift(image_spectrum .* conj(image_spectrum)))
axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'

% Plot the magnitude of the field at the image plane.
image_field = fft2(image_spectrum);
image_field_mag = image_field .* conj(image_field);
subplot(2,2,3)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
image_field_mag)
axis square, xlabel 'mm', ylabel 'mm'

% Plot the point spread function (PSF).
figure
plot(aperture_grid, ...
(image_field_mag(nsamples/2,:) / max(image_field_mag(nsamples/2,:))))
axis square, xlabel 'mm', ylabel 'normalized intensity'


All help is appreciated!

I made a couple of minor tweaks to your code and it seems to have "cleaned things up" a bit. The main thing I did was to shift the (sx,sy) by 1 sample (note the added -1 in the definition of aperture_spectrum_grid below) which seemed to have the most impact. The other change was to shift the aperture grid so that zero is actually in the sample.

When I did some diagnostics, I noticed that (using say 32 samples instead of 1024) the center of the spectrum was significantly offset from (0,0) so I suspected that perhaps getting the spectrum to be centered would at least get things more symmetric.

NEW: Added a second version attempting to align the aperture spectrum with the sx, sy grid. Please see below original code for 2nd code listing.

In any case I hope this helps.

clear all;
close all;
%% Example: Circular aperture.

lambda    = 632.8e-9; % 632.8nm wavelength light.
diameter  = 100e-6;   % 100um circular aperture diameter.
aperturez = 0;        % Aperture is located at z=0.
imagez    = 100e-2;   % 100cm from the aperture to the image plane.
imagez    = 0.001;  % try no-shift for sanity check

aperture_sample_width = 5e-4; % Sample the aperture for a width of 0.5mm.
% image_sample_width    = 2e-2; % Sample the image plane for a width of 2cm.

%nsamples = 1024; % Sample count (prefer power of 2).
nsamples = 128;

predict_airy_disc_radius = 1.22 * lambda * (imagez - aperturez) / diameter;
fprintf('Predicted Airy disc radius: %0.2f mm\n', ...

figure

% Plot the magnitude of the field at the (circular) aperture plane.
aperture_field = @(u,v) (sqrt((u.^2 + v.^2)) < (diameter/2));
delta_spacing = aperture_sample_width/nsamples
aperture_grid  = linspace(-aperture_sample_width/2+delta_spacing, ...
aperture_sample_width/2, nsamples);
[ugrid, vgrid] = meshgrid(aperture_grid);
subplot(2,2,1)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
aperture_field(ugrid,vgrid))
axis square, xlabel 'mm', ylabel 'mm'

% Plot the magnitude of the angular spectrum at the aperture.
aperture_spectrum = ifft2(aperture_field(ugrid,vgrid));
aperture_spectrum_spacing = 1/aperture_sample_width;
aperture_spectrum_grid = ((1:nsamples) - (nsamples/2)-1) * ...
aperture_spectrum_spacing;
[sx,sy] = meshgrid(aperture_spectrum_grid);
subplot(2,2,2)
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
ifftshift(aperture_spectrum .* conj(aperture_spectrum)))
axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'
%imagesc(ifftshift(aperture_spectrum .* conj(aperture_spectrum)))

% Plot the magnitude of the angular spectrum at the image plane.
image_spectrum = aperture_spectrum .* ...
exp((-1j * (2*pi/lambda) * (imagez - aperturez)) + ...
( 1j * pi*lambda * (sx.^2 + sy.^2) * (imagez - aperturez)));
subplot(2,2,4)
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
ifftshift(image_spectrum .* conj(image_spectrum)))
axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'

% Plot the magnitude of the field at the image plane.
image_field = fft2(image_spectrum);
image_field_mag = image_field .* conj(image_field);
subplot(2,2,3)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
image_field_mag)
axis square, xlabel 'mm', ylabel 'mm'

% Plot the point spread function (PSF).
figure
plot(aperture_grid, ...
(image_field_mag(nsamples/2,:) / max(image_field_mag(nsamples/2,:))))
axis square, xlabel 'mm', ylabel 'normalized intensity'


Created a new version which attempts to shift the aperture spectrum to align with the sx,sy grid, which I believe is centered in the middle of the matrix. This one seems to correctly align with the predicted Airy disk.

I hope this helps.

clear all;
close all;
%% Example: Circular aperture.

lambda    = 632.8e-9; % 632.8nm wavelength light.
diameter  = 100.0e-6;   % 100um circular aperture diameter.
aperturez = 0;        % Aperture is located at z=0.
imagez    = 10.0e-2;   % 100cm from the aperture to the image plane.

aperture_sample_width = 5000.0e-6; % Sample the aperture for a width of 0.5mm.
% image_sample_width  = 2e-2; % Sample the image plane for a width of 2cm.

%nsamples = 1024; % Sample count (prefer power of 2).
nsamples =1024;

predict_airy_disc_radius = 1.22 * lambda * (imagez - aperturez) / diameter;
fprintf('Predicted Airy disc radius: %0.2f mm\n', ...

figure
colormap(gray)

% Plot the magnitude of the field at the (circular) aperture plane.
aperture_field = @(u,v) (sqrt((u.^2 + v.^2)) < (diameter/2));
delta_spacing = aperture_sample_width/nsamples
aperture_grid  = linspace(-aperture_sample_width/2-delta_spacing, ...
aperture_sample_width/2, nsamples);
[ugrid, vgrid] = meshgrid(aperture_grid);
subplot(2,2,1)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
aperture_field(ugrid,vgrid))
axis square, xlabel 'mm', ylabel 'mm'

% Plot the magnitude of the angular spectrum at the aperture.
aperture_spectrum = ifft2(aperture_field(ugrid,vgrid));
aperture_spectrum_spacing = 1/aperture_sample_width;
aperture_spectrum_grid = ((1:nsamples) - (nsamples/2)-1) * ...
aperture_spectrum_spacing;
[sx,sy] = meshgrid(aperture_spectrum_grid);
subplot(2,2,2)
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
ifftshift(aperture_spectrum .* conj(aperture_spectrum)))
axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'
%imagesc(ifftshift(aperture_spectrum .* conj(aperture_spectrum)))

% Plot the magnitude of the angular spectrum at the image plane.
% Shift this aperture spectrum to align with sx, sy grid.
image_spectrum_2 = fftshift(aperture_spectrum).* ...
exp((-1j * (2*pi/lambda) * (imagez - aperturez)) + ...
( 1j * pi*lambda * (sx.^2 + sy.^2) * (imagez - aperturez)));

subplot(2,2,4)
%imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
%        aperture_spectrum_grid([1,end]) / 1e3, ...
%        ifftshift(image_spectrum .* conj(image_spectrum)))
imagesc(aperture_spectrum_grid([1,end]) / 1e3, ...
aperture_spectrum_grid([1,end]) / 1e3, ...
image_spectrum_2 .* conj(image_spectrum_2))

axis square, xlabel 'cycles / mm', ylabel 'cycles / mm'

% Plot the magnitude of the field at the image plane.
image_spectrum = fftshift(image_spectrum_2); % shift back to default
image_field = fft2(image_spectrum);
image_field_mag = image_field .* conj(image_field);
subplot(2,2,3)
imagesc(aperture_grid([1,end]) * 1e3, ...
aperture_grid([1,end]) * 1e3, ...
image_field_mag)
axis square, xlabel 'mm', ylabel 'mm'

% Plot the point spread function (PSF).
figure
plot(aperture_grid* 1e3, ...
(image_field_mag(nsamples/2,:) / max(image_field_mag(nsamples/2,:))))
axis square, xlabel 'mm', ylabel 'normalized intensity'

• Great, this works. I did later realize on my own that 1. the aperture grid needs to include (0,0) and 2. the aperture_spectrum needs to be shifted for the phase propagation multiplication. The book "Numerical Simulation of Optical Wave Propagation with Examples in MATLAB" by Jason Schmidt was also instrumental in helping me understand the constraints on the FFT sampling to avoid spatial frequency aliasing. – aosborne Jan 25 '20 at 16:09