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=100e-3;                 % distance

subplot(1,3,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);
%*********************************************************************
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(1,3,2),pcolor(xi*1e3,eta*1e3, abs(E1).^2), shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['z=' num2str(z*1e3) ' mm'])
subplot(1,3,3),pcolor(xi*1e3,eta*1e3,angle(E1).^2),shading interp
axis square
xlabel('xi  (mm)')
ylabel('eta (mm)')
title(['z=' num2str(z*1e3) ' mm'])
