function whol = asm(img,out=0,z0=0,z1=1,N=201)
% Put default values
z0
z1
N
out

if ischar(out)
  name=out;
elseif out == 0
  continue;
else
  error('asm: the last term shoud be string');
end

% Calculte domins
lambda = 630E-9
[y,x] = size(img);
p_size = 10E-6;  % pixel size
x_axis = x * p_size;  % ccd size (N* delta lambda)
y_axis = y * p_size; 
x_a = linspace(-round(x/2),round(x/2),x);
y_a = linspace(-round(y/2),round(y/2),y);
fx = linspace(-0.5/p_size,0.5/p_size,x);
fy = linspace(-0.5/p_size,0.5/p_size,y);
[Fx, Fy] = meshgrid(fx,fy);
p_x2= round(x/2);
p_y2= round(y/2);
p  = min(x,y)-2;

%% set phase retardation
d = linspace(z0,z1,N);      % d : reconstruction distance
dd= length(d);
% FFT of impulse response function   G
G=zeros(x,y);
fft_image=fftshift(fft2(img));

for ix = 1:length(d)
    G = (fft_image).*exp(1i*2*pi*d(ix)*sqrt(1/lambda^2-Fx.^2-Fy.^2)); % asm method
    recon = ifft2(ifftshift(G));
    shot= abs(recon); 
    if out == 0;
      continue;
    else
      figure(1);
      imagesc(x_a,y_a,shot);
      title(['z = ' num2str(d(ix))],"fontsize", 15);
      axis tight;
      fram = gcf;
      fname = [name "/out" num2str(ix,"%03d") ".png"];
      print (fram, fname);
    end
end
whol = recon;
