clc
clear all
close all

% Propagation de Fresnel 1D
% Selon Verrier et Atlan, Applied Optics 2011
        
lambda=633e-9;
N=1024*8;               % 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 ou de la gaussienne
Nf=1;                 % nombre de Fresnel

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

% calcul de la propagation d'un front d'onde
% en fonction du nombre de Fresnel


%*********************************************************************
z=w^2/(lambda*Nf);             % calcul de la distance z
%*********************************************************************
% Huygens-Fresnel (TFD unique)
[E1,deltaxi]=Fresnelprog1D_TFD_unique(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);
EF=fftshift(EF);
xi=((0:N-1)-N/2)/N*W/deltax*deltaxi;
% tracer l'amplitude
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(xi*1e3,abs(EF).^2,'k--')
axis([-1 1 0 2])
grid on
xlabel('x, xi (mm)')
ylabel('E')
title(['z=' num2str(z*1e3) ' mm, Nf=' num2str(Nf)])
legend('distribution en z=0', 'Huygens-Fresnel / 1TFD','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(xi*1e3,angle(EF)/pi,'k--')
axis([-1 1 0 2])
grid on
xlabel('x, xi (mm)')
ylabel('phase / pi')
title(['z=' num2str(z*1e3) ' mm'])
grid on