% Fresel direct method
% W = Fresnel_convolution(U,l,d,dpx,dpy)
% U - wave field
% l - wave length
% d - distance
% dpx - sample size horizontal
% dpy - sample size vertical
function W = Fresnel_convolution(U,l,d,dpx,dpy)
    % Number of columns (horizontal)
    Nx = size(U,2);
    % Number of rows (vertical)
    Ny = size(U,1);
    % Create coordinate maps
    [X,Y] = meshgrid(-Nx/2:Nx/2-1, -Ny/2:Ny/2-1);
    % xy dependent part
    A = exp(pi*i/d/l * ( (X.*dpx).^2 + (Y.*dpy).^2));
    % st dependent part
    B = exp(2*pi*i*d/l + pi*i*l*d*( (X/Nx/dpx).^2 + (Y/Ny/dpy).^2));
    % Perform calculation, plus an fftshift
    W = B .*fftshift(fft2( U .* A)) / (l*i*d);
    
return
