clear all
close all
clc

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;                  % longueur d'onde
z=1e-1;                         % distance de la propagation
colormap bone

% affichange de la distribution à z=0
subplot(2,2,1), imagesc(x*1e3,y*1e3,E0)
axis square
xlabel('x (mm)')
ylabel('y (mm)')
title('distribution à z=0')

% calcul de la propagation
xi =x;
eta=y;
%[E1,deltaxi,deltaeta]=Fresnelprog2D_TFD_unique (fftshift(E0),deltax,deltay,lambda,z);
%xi =((0:M-1)-(M-1)/2)*deltaxi;
%eta=((0:N-1)-(N-1)/2)*deltaeta;

%[E1]=Fresnelprog2D_3TFD(fftshift(E0),deltax,deltay,lambda,z);

[E1]=Fresnelprog2D_spectre_angulaire(fftshift(E0),deltax,deltay,lambda,z);

E1=fftshift(E1);

% affichage du résultat
subplot(2,2,2), imagesc(xi*1e3,eta*1e3,abs(E1).^2)
axis square
axis([-Wx/2 Wx/2 -Wy/2 Wy/2]*1e3)
xlabel('x (mm)')
ylabel('y (mm)')
title(['intensité à z= ' num2str(z*1e3) ' mm'])

subplot(2,2,4), imagesc(xi*1e3,eta*1e3,angle(E1))
axis square
axis([-Wx/2 Wx/2 -Wy/2 Wy/2]*1e3)
xlabel('x (mm)')
ylabel('y (mm)')
title(['phase à z= ' num2str(z*1e3) ' mm'])
