I am writing a code for generating Fraunhofer diffraction pattern from a circular slit by making use of 2D FFT. I am not getting Airy disk pattern as the output. I could not find what is wrong with my code?

The code is given below.

 clc
close all
clear all
format long
lambda = 500e-3;
w = 0.1;% radius of aperture
z = 4 * w^2 / lambda; % distance of the image plane from the aperture
start = -10.0;
stop = 10.0;
midPoint = (start + stop )/2;
stepSize = lambda * z * w;
x = start:stepSize:stop;
y = x;
[X, Y] = meshgrid(x,y);
E = zeros(length(x), length(y));
E(sqrt((X - midPoint).^2 + (Y - midPoint).^2) <= w) = 1;
G = fft2(E);
G = fftshift(G);
G = abs(G).^2;
G = G ./ max(max(G));
imagesc(G);
title('Fraunhofer Diffraction of circular aperture','fontsize',14)
colormap gray; colorbar; axis equal;
clear all

Answer 1

[David Goodmanson]

Hello Athul, It appears that the code is working more or less as it should. I don't understand where the step size expression comes from but it does all right in this case. If you cheat and plot abs(G) rather than the intensity abs(G).^2, you will see the pattern. Using the default colors rather than going to greyscale makes it easier to see, and looking at log10(abs(G).^2) gives you more diffraction rings than you will know what to do with.

You are taking abs(G) here, but it's better to shift the origin of coordinates down to where fft2 thinks it is with

 G = fft2(ifftshift(E))

so as to not introduce extraneous phase shifts into G. Incidentally there are 5001 = 3 x 1667 points in this fft. A number that is more composite would speed things up.

[Athul Shaji]

Hello David. Thanks for answering my question. There is an issue with the contrast. The low contrast is reason for the non visibility of the rings. See the 23rd line of the code ( imadjust) for the adjustment that I have made to make the rings visible.

 clc
close all
clear all
format long

 lambda = 500e-6;
w = 0.01;% radius of aperture
% z = 4 * w^2 / lambda;
z = 10;
start = -0.25;
stop = 0.25;
midPoint = (start + stop )/2;
% stepSize = lambda * z * w;
stepSize = (1/(2^12 - 1));
x = start:stepSize:stop;
y = x;
[X, Y] = meshgrid(x,y);
E = zeros(length(x), length(y));
E(sqrt((X - midPoint).^2 + (Y - midPoint).^2) <= w) = 1;
G = abs(fftshift(fft2(E))).^2;
G = G/max(max(G));
figure(2)
imagesc(imadjust(G, [0;.085],[0;1]));
title('Fraunhofer Diffraction of circular aperture','fontsize',14)
colormap gray; colorbar; axis equal;

[David Goodmanson]

Hello Athul, that's a good solution, although I don't have the right toolbox to see it myself.