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:
- 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$
- 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)}$
- 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', ...
1e3 * predict_airy_disc_radius);
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!
