clear all
close all
clc

E0=(double(imread('Fourier.bmp','bmp')));
E0=E0/max(max(E0));             % normalisation
Wx=2e-3;                        % définition des dimensions latérales
Wy=Wx;                          % 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

% affichage de la distribution initiale
imagesc(x*1e3,y*1e3,abs(E0).^2)
colormap bone
colorbar
%caxis([0 4])
axis square
axis([-Wx/2 Wx/2 -Wy/2 Wy/2]*1e3)
xlabel('x (mm)')
ylabel('y (mm)')
title(['intensité à z= 0 mm'])
pause
close


%*****************************************
%***** étape 1 - propagation de 3 cm *****
%*****************************************

% calcul de la propagation

z=3e-2;     % distance 3 cm

xi =x;              % Pour le cas où les coordonnées sont conservées
eta=y;              %
deltaxi=deltax;     %
deltaeta=deltay;    %

%[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,1,1), 
imagesc(xi*1e3,eta*1e3,abs(E1).^2)
colormap bone
colorbar
axis square
axis([min(xi) max(xi) min(eta) max(eta)]*1e3)
xlabel('xi (mm)')
ylabel('eta (mm)')
title(['intensité à z= ' num2str(z*1e3) ' mm'])

subplot(2,1,2), imagesc(xi*1e3,eta*1e3,angle(E1))
colormap bone
colorbar
axis square
axis([min(xi) max(xi) min(eta) max(eta)]*1e3)
xlabel('xi (mm)')
ylabel('eta (mm)')
title(['phase à z= ' num2str(z*1e3) ' mm'])

pause
close


%**********************************************************
%***** étape 2 - superposition de l'onde de référence *****
%**********************************************************

dec=1/4;                % décalage en fréquence, exprimé en
                        % tant que fraction de la largeur du spectre
dfx=dec/deltaxi;        % décalage en fréquence
dfy=dec/deltaeta;
thetax=asin(lambda*dfx);% calcul des angles d'inclinaison
thetay=asin(lambda*dfx);%     du faisceau de référence
xx=meshgrid(xi);
yy=meshgrid(eta)';
ER=exp(1i*2*pi/lambda*(sin(thetax)*xx+sin(thetay)*yy));

imagesc(xi*1e3,eta*1e3,angle(ER))
colormap bone
colorbar
axis square
axis([min(xi) max(xi) min(eta) max(eta)]*1e3)
xlabel('xi (mm)')
ylabel('eta (mm)')
title('phase de l''onde de référence')

pause
close

E1=E1/(max(max(abs(E1))));

E2=E1+ER;

% affichage du résultat
imagesc(xi*1e3,eta*1e3,abs(E2).^2)
colormap bone
colorbar
axis square
axis([min(xi) max(xi) min(eta) max(eta)]*1e3)
xlabel('xi (mm)')
ylabel('eta (mm)')
title(['hologramme enregistré à z= ' num2str(z*1e3) ' mm (log)'])

pause
close


%******************************************************
%***** étape 3 - calcul de la FFT de l'hologramme *****
% et sélection de la zone d'intérêt               *****
%******************************************************

HLG = fftshift(fft2(fftshift(abs(E2.^2))));   % FFT de l'intensité
fx=((0:M-1)/(M-1)-0.5)/deltaxi;
fy=((0:N-1)/(N-1)-0.5)/deltaeta;
imagesc(fx,fy,log10(abs(HLG).^2))
colormap bone
colorbar
axis square
axis([min(fx) max(fx) min(fy) max(fy)])
xlabel('fx (1/m)')
ylabel('fy (1/m)')
title('FFT de l''hologramme (log)')

pause
close

%**************************************
%***** choix de la zone d'intérêt *****
%**************************************

ROI=HLG( (M*0.5*dec-1:M*1.5*dec) , (N*0.5*dec-1:N*1.5*dec));
imagesc(log10(abs(ROI)))
colormap bone
colorbar
title('zone d''intérêt (log)')

pause
close

%*****************************************************
%****** FFT inverse et reconstruction de l'objet *****
%*****************************************************

IMR=fftshift(ifft2(fftshift(ROI))); % calcul de l'image réelle

%Oreconstr=fftshift(Fresnelprog2D_3TFD(fftshift(IMR),deltaxi*g,deltaeta*g,lambda,-z));

% rétropropagation vers le plan objet
Oreconstr=fftshift(Fresnelprog2D_spectre_angulaire ...
    (fftshift(IMR),deltaxi/dec,deltaeta/dec,lambda,-z));
subplot(2,1,1),imagesc(abs(IMR).^2);
axis square
colormap bone
title('front d''onde reconstruit - intensité')

subplot(2,1,2),imagesc(angle(IMR));
axis square
colormap bone
title('front d''onde reconstruit - phase')

pause
close all

subplot(2,1,1),imagesc(abs(Oreconstr).^2);
axis square
colormap bone
title('front d''onde repropagé vers le plan objet - intensité')

subplot(2,1,2),imagesc(angle(Oreconstr));
axis square
colormap bone
title('front d''onde repropagé vers le plan objet - phase')