clc
clear all
close all

% Propagation de Fresnel 1D
% Selon Verrier et Atlan, Applied Optics 2011

E0=rot90(rot90(double(imread('Fourier.bmp','bmp'))));
%E0=(double(imread('trou.bmp','bmp')));
E0=E0/max(max(E0));             % normalisation
Wx=1e-2;                        % définition des dimensions latérales
Wy=1e-2;                        % physiques du champs
[N,M]=size(E0);                 % dimensions de la matrice
deltax=Wx/M;                    % échantillonnage
deltay=Wy/N;                    % échantillonnage
x=((0:M-1)-(M-1)/2)*deltax;     % vecteur x
y=((0:N-1)-(N-1)/2)*deltay;     % vecteur y

lambda=633e-9;

%z=50e-3;                 % distance
Nf=100;
z=Wx^2/(Nf*lambda);

subplot(2,4,1),pcolor(x*1e3,y*1e3, E0), shading interp
axis square
xlabel('x (mm)')
ylabel('y (mm)')
title('z=0')

% calcul de la propagation d'un front d'onde
% en fonction de la distance

%*********************************************************************
% Huygens-Fresnel (TFD unique)
[E1,deltaxi,deltaeta]=Fresnelprog2D_TFD_unique(fftshift(E0),deltax,deltay,lambda,z);
% trois TFD
[E2]=Fresnelprog2D_3TFD(fftshift(E0),deltax,deltay,lambda,z);
% spectre angulaire
[E3]=Fresnelprog2D_spectre_angulaire(fftshift(E0),deltax,deltay,lambda,z);
%*********************************************************************
E1=fftshift(E1);
xi =((0:N-1)-N/2)/N*Wx/deltax*deltaxi;
eta=((0:N-1)-N/2)/N*Wy/deltay*deltaeta;
% afficher l'amplitude et la phase
subplot(2,4,2),pcolor(xi*1e3,eta*1e3, abs(E1).^2), shading interp
axis square
axis([min(x) max(x) min(y) max(y)]*1e3)
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['Huygens-Fresnel, z=' num2str(round(z*1e3,1)) ' mm'])
subplot(2,4,6),pcolor(xi*1e3,eta*1e3,angle(E1)   ),shading interp
axis square
axis([min(x) max(x) min(y) max(y)]*1e3)
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['z=' num2str(round(z*1e3,1)) ' mm'])
%*********************************************************************
E2=fftshift(E2);
% afficher l'amplitude et la phase
subplot(2,4,3),pcolor(x*1e3,y*1e3, abs(E2).^2), shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['3 TFD, z=' num2str(round(z*1e3,1)) ' mm'])
subplot(2,4,7),pcolor(x*1e3,y*1e3,angle(E2)   ),shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['z=' num2str(round(z*1e3,1)) ' mm'])
%*********************************************************************
E3=fftshift(E3);
% afficher l'amplitude et la phase
subplot(2,4,4),pcolor(x*1e3,y*1e3, abs(E2).^2), shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['spectre angulaire, z=' num2str(round(z*1e3,1)) ' mm'])
subplot(2,4,8),pcolor(x*1e3,y*1e3,angle(E3)   ),shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['z=' num2str(round(z*1e3,1)) ' mm'])