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.

Answer 1

[Sandeep Mishra]

Hi Sohail,

Upon executing the code snippet in MATLAB R2024a, I encountered the same outcome. During the debugging process, I identified that the 'rectangular_slit' function within the 'Sheet' class is not updating the 'E' variable as expected.

To resolve this issue, you should modify the 'rectangular_slit' function to return the modified object. Additionally, ensure that the 'Sheet' object is updated accordingly after invoking the 'rectangular_slit' function in your main script.

Please update the "rectangular_slit" method of "Sheet" class as follows:

% Updated function for sheet.m file 
function obj = 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 

Please update the call to "rectangular_slit" in the main code as follows:

% Updation function call for main file 
% Add two slits to the sheet 
sheet = sheet.rectangular_slit(-D/2, 0, 22 * mm, 88 * mm); 
sheet = sheet.rectangular_slit(D/2, 0, 22 * mm, 88 * mm); 

For more information, refer to the following MathWorks documentation:

1. ‘Method Syntax’: https://www.mathworks.com/help/releases/R2024a/matlab/matlab_oop/specifying-methods-and-functions.html
2. ‘Using a Class to Display Graphics’: https://www.mathworks.com/help/releases/R2024a/matlab/matlab_oop/a-class-code-listing.html#buq2wcy

I hope this helps.

[Sohail Ahmed]

Thank you so much. I got it ...!!