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;