Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB FFT problem for Fresnel simulation

  1. May 29, 2018 #1
    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 :
    Code (Python):
    %%%%%%    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;
    for i_grating = 251:350;
        A_t(i_grating) = 1;
    for i_grating = 351:601
        A_t(i_grating) = 0;

    %% 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
        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));

    %% 2D Map %%

    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 ,


    Attached Files:

  2. jcsd
  3. Jun 3, 2018 #2
    Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?