FFT problem for Fresnel simulation

[CijiG]

Hi all,

I'm trying to simulate the Fresnel diffraction by using this expersion :
$$ A(x')=\frac{1}{j\lambda z}e^{jk(z+\frac{x'^2}{2z})}F(A^{trans}(x)e^{jk\frac{x^2}{2z}})_{u=\frac{x'}{\lambda z}} $$

So when I use this formula my problem is that I don't know how to take the good frequency, here $\frac{x'}{\lambda z}$

and in my simulation, the result is that the light is not diverging after the rectangular aperture. I don't know how to fix this problem, because I wan't to plot the (x,z) plan so z is not constant. When I plot only one column, which corresponf to a z fix, I've the good profile of intensity but it's steel not diverging when I choose different z value. So here is my Matlab (R2017b) simulation code :

%%%%%%    Simulation of the Fresnel Diffraction, Matlab    %%%%%%%

%% Initialization %%

clear all
close all

%% Grating %%

for i_grating = 1:250            % creation of a rectangular aperture
    A_t(i_grating) = 0;
end
for i_grating = 251:350;
    A_t(i_grating) = 1;
end 
for i_grating = 351:601
    A_t(i_grating) = 0;
end 


%% Parameters %%
     
Lambda = 795e-9;                % wavelength
k = 2*pi/Lambda;                % the wave vector
z = [0.06 :0.01: 6];            % distance from the grating
dx = 20e-6;                     % size of a pixel in micrometer
x =[-6000e-6 :dx: 6000e-6];     % coordinate on the grating plane
x2 =[-6000e-6 :dx: 6000e-6];    % coordinate on the screen plane

plot(A_t)                       % show the amplitude just after the grating(z=0+)

%% Simulation %%

for i_z = 1:length(z)             % first loop on the distance to the grating
    for i_x = 1:length(x)         % loop on all point of the grating
          f(i_x)=A_t(i_x)*exp(1i*k*((x(i_x))^2)/(2*z(i_z)));
    end
    A = fftshift(fft(f));
    for i_x2 = 1:length(x2)         % we do a loop for all point in the observing plane at z = cst
        C =(1/(1i*Lambda*z(i_z)))*exp(1i*k*(z(i_z)+(x2(i_x2))^2)/(2*z(i_z)));
        A_p(i_x2,i_z) = abs(C*A(i_x2));
    end
   
end

%% 2D Map %%

figure;
imagesc(A_p);
colormap jet;         % color representing the light intensity

I used python color to make it easier to read, but the % are not seen as comments sorry. And here you have the results :

On the right, what we see is the light intensity, the x-axis correspond to the vertical axis and the z axis to the horizontal axis.

thanking you ,

CijiG