function [out, TF] = propfield(in,dx,lambda,z,type)
%PROPFIELD Propagate complex light field
%
% out = PROPFIELD(in, dx, lambda, z, type);
% [out, TF] = PROPFIELD(in, dx, lambda, z, type);
%
% Explanation of variables (all SI units):
% in - source plane field (2D or 1D)
% dx - pixel size (assuming square pixels for 2D)
% lambda - wavelength
% z - propagation distance
% type - propagation type:
%   'FTF' : Fresnel transfer function
%   'RSTF': Rayleigh-Sommerfeld transfer function
%   'RSIR': Rayleigh-Sommerfeld impulse response
%   'KIR' : Kirchhoff transfer function
%   'AS'  : Angular spectrum approach (TF)
% out - observation plane field
% TF - transfer function for reuse
%
% Example:
%  x = (-499:500)*100e-9; [nd1, nd2] = ndgrid(x);     %grid: 100 nm    
%  in = double(sqrt(nd1.^2+nd2.^2)<=5e-6);     %input: circular aperture, 5?m dia
%  out = propfield(in, 100e-9, 500e-9, 10e-6, 'AS');     %propagating 10?m distance (for 500 nm light)
%  figure, subplot(2,1,1), imagesc(x, x, abs(in).^2, [0 2]), axis equal tight, colormap gray, title('input')
%  subplot(2,1,2), imagesc(x, x, abs(out).^2, [0 2]), axis equal tight, colormap gray, title('output')
%
%
% out = PROPFIELD(in, TF);
%   Reuse a previously output TF (for faster calculation)

% See http://dx.doi.org/10.1016/j.optcom.2007.08.053

% 130329 - created, BJ
% 131205 - corrected major bug for TF approaches using nonsquare fields
% 150505 - corrected bug for odd sized TF (now using 'ceil')


if nargin > 2
    
    if z == 0 
        out = in;
        return
    end
    
    [M,N] = size(in); 
    if isa(in, 'gpuArray'), dx = complex(gpuArray(dx)); end
    %[fM,fN] = ndgrid((-M/2:M/2-1)./(M*dx), (-N/2:N/2-1)./(N*dx)); %bug fixed: M to N in ".../(M*dx));"
    [fM,fN] = ndgrid((ceil(-M/2):ceil(M/2)-1)./(M*dx), (ceil(-N/2):ceil(N/2)-1)./(N*dx)); % bug fixed: ceil as opposed to (-M/2:M/2-1)
    [ndM, ndN] = ndgrid((ceil(-M/2):ceil(M/2)-1)*dx, (ceil(-N/2):ceil(N/2)-1)*dx);
    k = 2*pi/lambda;
    signZ = sign(z); z = abs(z);

    switch(type)

        case 'FTF' % Fresnel transfer function
            TF = exp(1i*k*z*( 1-((fM.^2+fN.^2)*lambda^2/2) )); TF = ifftshift(TF);

        case 'RSTF' % Rayleigh-Sommerfeld transfer function
            TF = exp(1i*k*z).*exp(-1i*pi*lambda*z.*(fN.^2+fM.^2)); TF = ifftshift(TF);

        case 'RSIR' % Rayleigh-Sommerfeld impulse response
            r = sqrt(ndM.^2+ndN.^2+z.^2);
            IR = z.*exp(1i*k*r)./(r.^2*1i*lambda); 
            TF = fft2(ifftshift(IR))*dx^2;
            
        case 'RSIR1D'
            r = sqrt(ndM.^2+ndN.^2+z.^2);
            IR = exp(1i*k*r)./(sqrt(r)*1i*lambda).*(z./r);  
            TF = fft2(ifftshift(IR))*dx^2;
            
        case 'HUYGENS'
            r = sqrt(ndM.^2+ndN.^2+z.^2);
            IR = exp(1i*k*r)./r;
            TF = fft2(ifftshift(IR))*dx^2;

        case 'KIR' %Kirchhoff transfer function
            r = sqrt(ndM.^2+ndN.^2+z.^2);
            IR = exp(1i*k*r).*(z+r)./(1i*2*lambda*r.^2); 
            if sign(z)==-1, IR = conj(IR); end
            TF = fft2(ifftshift(IR))*dx^2;

        case 'AS' % Angular spectrum approach
            %TF = exp(1i*k*z*sqrt( 1-((fM.^2+fN.^2)*lambda^2) )); TF = ifftshift(TF);
            TF = exp(1i*k*z*real(sqrt( 1-((fM.^2+fN.^2)*lambda^2) ))); TF = ifftshift(TF);
    end

    if signZ == -1
        TF = conj(TF);
    end
    
elseif nargin == 2 % if reusing TF
    TF = dx; %copy second argument (see function syntax above)
end %if

ftIn = fft2(ifftshift(in));
out = fftshift(ifft2(TF .* ftIn)); 


%[fM,fN] = ndgrid( ((1:M)-(M+1)/2)/(M*dx) , ((1:N)-(N+1)/2)/(N*dx) );
%fM=-1/(2*dx):1/(M*dx):1/(2*dx)-1/(M*dx); %freq coords
%fN=-1/(2*dx):1/(N*dx):1/(2*dx)-1/(N*dx); %freq coords
%TF = exp(-1i*pi*lambda*z*(fM.^2+fN.^2)); TF = fftshift(TF);
%TF = exp(1i*k*z.*sqrt(1-(lambda*fM).^2-(lambda*fN).^2)); TF = fftshift(TF);
