function output = Propagate(input,propagator,dx,lamda,z)
% 波前传播
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input: the wavefront to propagate
% propagator: one of 'fourier','fresnel' or 'angular spectrum'
% dx: the pixel spacing of input wavefront
% wavelength: the wavelength of the illumination
% z: the distance to propagate
% output: the distance to propagate
% 频域坐标
k = 2*pi/lamda;
[ysize, xsize] = size(input);
x = -xsize/2:xsize/2-1;
y = -ysize/2:ysize/2-1;
fx = x./(dx*xsize);
fy = y./(dx*ysize);
[fx, fy] = meshgrid(fx,fy);
switch propagator
case 'fourier'
if z>0
output = fftshift(fft2(fftshift(input)));
% output = fft2(input);
else
output = ifftshift(ifft2(ifftshift(input)));
% output = ifft2(input);
end
case 'angular spectrum'
% Calculate phase distribution for each plane wave component
w = sqrt(1/lamda^2 - fx.^2 - fy.^2);
% exclude evanescent waves
notEvanescent = imag(w) == 0;
% Compute FFT of input
F = fftshift(fft2(fftshift(input)));
% multiply FFT by phase-shift and inverse transform
output = ifftshift(ifft2(ifftshift(F.*exp(2i*pi*z*w).*notEvanescent)));
case 'fresnel'
% Calculate approx phase distribution for each plane wave component
w = fx.^2 + fy.^2;
% Comput FFT
F = fftshift(fft2(fftshift(input)));
% multiply by phase-shift and inverse transform
output = ifftshift(ifft2(ifftshift(F.*exp(-1i*pi*z*lamda*w))));
end
end
