function [Eout, deltaxout, deltayout] = Fresnelprog2D_TFD_unique (Ein,deltax,deltay,lambda,z)

% Calcul de la propagation de Fresnel selon l'algorithme Huygens-Fresnel /
% TFD unique

[N,M]=size(Ein);
p=[(0:(M/2)-1) (-M/2:-1)];
q=[(0:(N/2)-1) (-N/2:-1)]';
[pp,qq]=meshgrid(p,q);
Eout=deltax*...
    exp(1i*2*pi/lambda*z)/sqrt(1i*lambda*z)*...
    exp(1i*pi*lambda*z*(pp.^2/(M^2*deltax^2)+qq.^2/(N^2*deltax^2))).*...
    fft2(Ein.*exp(1i*pi*((pp*deltax).^2+(qq*deltax).^2)/(lambda*z)));

deltaxout=lambda*z/((M-1)*deltax);
deltayout=lambda*z/((N-1)*deltay);