clc
clear all
close all

% Propagation de Fresnel 1D
% Selon Verrier et Atlan, Applied Optics 2011
% Calcule et affiche le front d'onde derrière une fente pour différentes
% distances
        

lambda=633e-9;
N=1024;                 % nombre d'éléments dans le vecteur
W=1000e-6;              % largeur du champ
x=(0:N-1)/(N-1)*W-W/2;  % vecteur x
deltax=W/(N-1);         %
w=50e-6;                % mi-largeur de la fente

minz=100e-6;
maxz=100e-6+minz;         % distance maximale
maxk=N;               % nombre de positions en z

E0=fftshift((abs(x)<w)*1.0);      % distribution du champ à z=0 (rectangle)
%E0=(abs(x)<w).*(w-abs(x))/w; % distribution du champ à z=0 (triangle)
%E0=fftshift(exp(-(x/w).^2));  % distribution du champ à z=0 (gaussienne)
matE=zeros(N,maxk);

for k=1:maxk
    z=minz+(maxz-minz)*k;
    %*********************************************************************
    % Huygens-Fresnel (TFD unique)
    %[E1,deltaxi]=Fresnelprog1D_TFD_unique(E0,deltax,lambda,z);
    %E1=fftshift(E1);
    % Convolution (trois TFD)
    E2=Fresnelprog1D_3TFD(E0,deltax,lambda,z);
    E2=fftshift(E2);
    % Spectre angulaire
    %E3=Fresnelprog1D_spectre_angulaire (E0,deltax,lambda,z);
    %E3=fftshift(E3);
    % Transformée de Fourier
    %EF=deltax*...
    %    exp(1i*2*pi/lambda*z)/sqrt(1i*lambda*z)*...
    %    fft(E0);
    %EF=fftshift(EF);
    %*********************************************************************
    %xi=((0:N-1)-N/2)/N*W/deltax*deltaxi;
    matE(:,k)=E2;
end
zax=(1:maxk)/maxk*(maxz-minz)+minz;
pcolor(zax*1e3,x*1e3,abs(matE.^2)), shading interp
%imagesc(zax*1e3,x*1e3,abs(matE.^2))
colorbar
ylabel('position latérale x (mm)')
xlabel('distance z (mm)')