% Fresnel convolution
% 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);
    % Fourier transform, and shift, wave field
    U = fftshift(fft2(U));
    % Compute the kernel
    T = exp(-pi*i*l*d* ((X/dpx/Nx).^2 + (Y/dpy/Ny).^2));
    % Multiply by kernel and constants, invers FFT
    W = ifft2(U.*T.*exp(2*pi*i*d/l)) / (i*l*d);
return
