function [u2] = propASM(u1,pitch,lambda,z)
%========================================================================
%Copyright(c) 2017 EmergIMG project
%Universidade da Beira Interior - Instituto de Telecomunicações, Portugal
%All Rights Reserved.
%
%This is an implementation of the Angular Spectrum Reconstruction Method for holographic data.
%Transfer function approach based on Rayleigh-Sommerfeld formula. 
%Assumes x and y side length and uniform sampling.
%Described in the following paper:
%
%[1] Marco V. Bernardo, Pedro Fernandes, Angelo Arrifano, Marc Antonini, Elsa Fonseca, Paulo T. Fiadeiro, António M. G. Pinheiro and Manuela Pereira,
% "Holographic Representation: Hologram Plane vs. Object Plane", submitted, 2018
%
%Input : u1     - Hologram Plane Field
%        pitch  - Pixel Pitch
%        lambda - Wavelength
%        z      - Propagation Distance
%Output: u2 - Object Plane Field
%
%Example for 2D Dice Synthetic Hologram:
% load('Hol_2D_dice.mat')
%
%
% [u2] = propASM(u1,pitch,lambda,z);
% imshow(abs(u2),[])
%
%========================================================================


[M,N] = size(u1);
dx = pitch;
dy = pitch;
Lx = N*pitch;
Ly = M*pitch;

fx = -1/(2*dx):1/Lx:1/(2*dx)-1/Lx;          % freq coordinates (x)
fy = -1/(2*dy):1/Ly:1/(2*dy)-1/Ly;          % freq coordinates (y)

[FX,FY] = meshgrid(fx,fy);
H = exp( 2i * pi * z*((1/lambda)^2-FX.^2 - FY.^2).^0.5);	% trans function

U1 = fftshift(ifft2(fftshift(u1)));       % shift src, fft src field
U2 = U1.*conj(H);                         % multiply
u2 = fftshift(fft2(fftshift(U2)));        % inv fft, center obs field

end