How to compute diffraction integral using fast Fourier transform (fft2) ?

Hi there !

The following code computes Fresnel diffraction integral using fft2. It works for (obj.E = ones(Ny, Nx);) and performs the fft2 analysis, but when the input is (obj.E = zeroes(Ny, Nx);), it plots nothing. I checked the same code on python with zeros and it worked but in MATLAB it is not performing fft2 for zeros (gives only zero values). Kindly help me resolve this issue.

The Sheet.m (code):

classdef Sheet

properties

dx

dy

x

y

xx

yy

Nx

Ny

E

end

methods

function obj = Sheet(extent_x, extent_y, Nx, Ny)

obj.dx = extent_x / Nx;

obj.dy = extent_y / Ny;

obj.x = obj.dx * ((0:Nx-1) - floor(Nx/2));

obj.y = obj.dy * ((0:Ny-1) - floor(Ny/2));

[obj.xx, obj.yy] = meshgrid(obj.x, obj.y);

obj.Nx = Nx;

obj.Ny = Ny;

obj.E = zeros(Ny, Nx);

end

function rectangular_slit(obj, x0, y0, lx, ly)

% Creates a slit centered at the point (x0, y0) with width lx and height ly

mask = ((obj.xx > (x0 - lx/2)) & (obj.xx < (x0 + lx/2))) & ((obj.yy > (y0 - ly/2)) & (obj.yy < (y0 + ly/2)));

obj.E = obj.E + mask;

end

end

end

The main code :

% Simulation input

Lx = 1.4;

Ly = 0.4;

Nx = 2500;

Ny = 1500;

sheet = Sheet(2*Lx, 2*Ly, Nx, Ny);

% Slit separation

mm = 1e-3;

D = 128 * mm;

sheet.rectangular_slit(-D/2, 0, 22 * mm, 88 * mm);

sheet.rectangular_slit(D/2, 0, 22 * mm, 88 * mm);

% Distance from slit to the screen (mm)

z = 5000;

% Wavelength (mm)

lambda = 18.5*1e-7;

k = 2*pi / lambda;

% Initialize complex array for phase information

fft_c = fft2(sheet.E .* exp(1i * k/(2*z) * (sheet.xx.^2 + sheet.yy.^2)));

c = fftshift(fft_c);

% Plot with MATLAB

abs_c = abs(c);

% Screen size (mm)

dx_screen = z*lambda/(2*Lx);

dy_screen = z*lambda/(2*Ly);

x_screen = dx_screen * ((1:Nx) - Nx/2);

y_screen = dy_screen * ((1:Ny) - Ny/2);

figure;

subplot(2, 1, 1);

imagesc(x_screen([1 end]), y_screen([1 end]), abs_c);

colormap(gray);

axis equal;

xlabel('x (mm)');

ylabel('y (mm)');

title('Probability Density |\psi|^2');

xlim([-2, 2]);

ylim([-1, 1]);

subplot(2, 1, 2);

plot(x_screen, abs(c(Ny/2, :)).^2, 'LineWidth', 1);

xlabel('x (mm)');

ylabel('Probability Density |\psi|^2');

title('Probability Density |\psi|^2 along y=0');

xlim([-2, 2]);

[Sohail Ahmed]

Anyone here to answer ??

[David Goodmanson]

Hi SK. if the input is all zeros, is there a reason that you expect that the output would not be all zeros? If anyting I might be questioning the python code that gives nonzero output for zero input.

[Sohail Ahmed]

Yeah, there is a reason. The issue arises with the zero padding of the Fresnel inetgral while applying fft2.