close all; 
clear; 
clc ;

%% Synthesized modulus data: m=abs(Fresnel Propagator*Psi)
% Define circular aperture function
close all; clear; clc

% Aperture function definition
circ=@(r) (abs(r)<=1.) ;

%% Global parameters
N=512;
L=2048*1.e-6 ;        % 2048-um simulation size at z=0 cm
dx=L/N;                    % dx (sampling interval, m) 
x=-L/2:dx:L/2-dx; 
y=x; 

ndx=1.52;              % index of refraction
R=100*1.e-3 ;       % curvature radius of a plano-convex lens
fl=R/(ndx-1);       % Effective focal length of the lens
a = L/4;                   % aperture radius of the lens

%% Prepare field amplitude profile
lambda=600*1e-9 ;    % wavelength (m)
k=2.0*pi/lambda;      % wavenumber (rad/m)

%% Prepare phase profile of a thin lens
rn=zeros(N,N);
aperture=zeros(N,N);
for i =1:N
    rn(:,i)=sqrt(x(i).^2+y(:).^2)/a;
end
aperture = double(circ(rn));
Phi = zeros(N,N);
Phi =  k*(a*rn).^2/2/fl ;
Phi = aperture.*Phi;
field0 = aperture.*exp(-1j*Phi);

propz=1.5*fl;
dz=propz/40;
field1 = zeros(size(field0)); field2 = zeros(size(field0));
frame=1;
field1=propagate(field0,lambda,propz,dx);

for ii=1:40
    z=propz-ii*dz;
    field2=propagate(field1,lambda,-dz,dx);
    dp=(abs(field2)).^2;    %get obj hologram
    imagesc(dp);
    colormap gray; 
    axis off; 
    colorbar
    title(sprintf('z=%1.1f cm',100*z))
    drawnow; pause(0.1); 

    F(frame)=getframe;
    frame=frame+1;
    field1=field2;
end
% VideoWriter('holoZs.avi');