clear all
close all
clc

E0=double(imread('rectangle.bmp','bmp'));
E0=E0/max(max(E0));             % normalisation
Wx=1e-3;                        % définition des dimensions latérales
Wy=1e-3;                        % 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
zmin=1e-4;                      % distance minimale de propagation
zmax=0.18;                      % distance maximale de propagation
Nframe=100;                     % nombre d'intervalles

figure('Renderer','zbuffer');
set(gca,'NextPlot','replaceChildren');

% affichange de la distribution à z=0
imagesc(x*1e3,y*1e3,E0)
colormap bone
colorbar
axis square
axis([-Wx/2 Wx/2 -Wy/2 Wy/2]*1e3)
xlabel('x (mm)')
ylabel('y (mm)')
title('distribution à z=0')

%pause

F(1)=getframe;


for l=1:Nframe
    z=zmin+(zmax-zmin)*l/Nframe;
    % 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
    imagesc(xi*1e3,eta*1e3,abs(E1).^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= ' num2str(z*1e3) ' mm'])
    F(l)=getframe;
end

movie(F,20)