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érents
% nombres de Fresnel (différentes distances) en utilisant trois algorithmes
        

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

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)

% calcul de la propagation pour différentes valeurs du nombre de Fresnel

Nf=[5 2 1 0.5 0.2 0.1 0.05 0.02 0.01 0.005 0.002 0.001 ...
    0.0005 0.0002 0.0001];

for k=1:length(Nf)
    z=w^2/(lambda*Nf(k));          % calcul de la distance z
                                   % à partir du nombre de Fresnel
    %*********************************************************************
    % Huygens-Fresnel (TFD unique)
    [E1,deltaxi]=Fresnelprog1D_TFD_unique(E0,deltax,lambda,z);
    % Convolution (trois TFD)
    E2=Fresnelprog1D_3TFD(E0,deltax,lambda,z);
    % Spectre angulaire
    E3=Fresnelprog1D_spectre_angulaire (E0,deltax,lambda,z);
    % Transformée de Fourier
    EF=deltax*...
        exp(1i*2*pi/lambda*z)/sqrt(1i*lambda*z)*...
        fft(E0);
    %*********************************************************************
    E1=fftshift(E1);
    E2=fftshift(E2);
    E3=fftshift(E3);
    EF=fftshift(EF);
    xi=((0:N-1)-N/2)/N*W/deltax*deltaxi;
    % tracer l'amplitude
    figure('units','normalized','outerposition',[0 0 1 1])
    subplot(2,1,1),plot(x*1e3,abs(fftshift(E0)).^2,'k-')
    hold on
    subplot(2,1,1),plot(xi*1e3,abs(E1).^2,'r-')
    subplot(2,1,1),plot(x*1e3,abs(E2).^2,'b-')
    subplot(2,1,1),plot(x*1e3,abs(E3).^2,'g-')
    subplot(2,1,1),plot(xi*1e3,abs(EF).^2,'k--')
    axis([-0.5 0.5 0 2])
    grid on
    xlabel('x, xi (mm)')
    ylabel('E')
    title(['z=' num2str(z*1e3) ' mm, Nf=' num2str(Nf(k))])
    legend('distribution en z=0', ...
        'Huygens-Fresnel / 1TFD', 'Convolution / 3TFD', ...
        'Spectre angulaire', 'Fourier')
    grid on
    % tracer la phase
    subplot(2,1,2),plot(x*1e3,angle(fftshift(E0))/pi,'k-')
    hold on
    subplot(2,1,2),plot(xi*1e3,angle(E1)/pi,'r-')
    subplot(2,1,2),plot(x*1e3,angle(E2)/pi,'b-')
    subplot(2,1,2),plot(x*1e3,angle(E3)/pi,'g-')
    %subplot(2,1,2),plot(xi*1e3,angle(EF)/pi,'k--')
    axis([-0.5 0.5 -1 1])
    grid on
    xlabel('x, xi (mm)')
    ylabel('phase / pi')
    title(['z=' num2str(z*1e3) ' mm'])
    grid on
    %norm=sum(abs(Exi).^2)*deltaxi/(sum(abs(Ex).^2)*deltax)
    pause
    close
end