;+
; NAME:
;     rayleighsommerfeld
;
; PURPOSE:
;     Computes Rayleigh-Sommerfeld back-propagation of
;     a normalized hologram of the type measured by
;     digital video microscopy
;
; CATEGORY:
;     Holographic microscopy
;
; CALLING SEQUENCE:
;     b = rayleighsommerfeld(a,z)
;
; INPUTS:
;     a: hologram recorded as image data normalized
;         by a background image. 
;     z: displacement from the focal plane [pixels]
;         If z is an array of displacements, the field
;         is computed at each plane.
;
; KEYWORDS:
;     lambda: Wavelength of light in medium [micrometers]
;         Default: 0.632 -- HeNe in air
;     mpp: Micrometers per pixel
;         Default: 0.135
;
; KEYWORD_FLAGS:
;     nozphase: By default, a phase factor of kz is removed 
;         from the propagator to eliminate the (distracting)
;         axial phase gradient.  Setting this flag leaves 
;         this factor (and the axial phase ramp) in.
;     hanning: If set, use Hanning window to suppress
;         Gibbs phenomeon.
;
; OUTPUTS:
;     b: complex field at height z above image plane,
;         computed by convolution with
;         the Rayleigh-Sommerfeld propagator.
;
; REFERENCES:
; 1. S. H. Lee and D. G. Grier, 
;   "Holographic microscopy of holographically trapped 
;   three-dimensional structures,"
;   Optics Express, 15, 1505-1512 (2007).
;
; 2. J. W. Goodman, "Introduction to Fourier Optics,"
;    (McGraw-Hill, New York 2005).
;
; 3. G. C. Sherman,
;   "Application of the convolution theory to Rayleigh's integral
;   formulas,"
;   Journal of the Optical Society of America 57, 546-547 (1967).
;
; PROCEDURE:
;     Convolution with Rayleigh-Sommerfeld propagator using
;     Fourier convolution theorem.
;
; EXAMPLES:
; Visualize the intensity as a function of height z above the image
; plane given a normalized holographic image a.

; IDL> lambda = 0.6328 / 1.336 ; HeNe in water
; IDL> for z = 1., 200 do begin $
; IDL>    tvscl, abs(rayleighsommerfeld(a, z, lambda=lambda))^2 & $
; IDL> endfor
;
; Visualize the phase singularity due to a particle located along the
; line y = 230 in the image a, assuming that the particle is no more
; than 200 pixels above the image plane.
;
; IDL> help, a
; A     FLOAT     = Array[640, 480]
; IDL> phi = fltarr(640,200)
; IDL> for z = 1., 200 do begin $
; IDL> phiz = atan(rayleighsommerfeld(a,z,lambda=lambda),/phase) & $
; IDL> phi[*,z-1] = phiz[*,230] & $
; IDL> endfor
; IDL> tvscl, phi
;
; MODIFICATION HISTORY:
; 11/28/2006: Sang-Hyuk Lee, New York University, 
;   Original version (called FRESNELDHM.PRO).
;   Based on fresnel.pro by D. G. Grier
; 10/21/2008: David G. Grier, NYU.  Adapted from FRESNELDHM
;   to provide more uniform interface.
; 11/07/2009: DGG.  Code and documentation clean-up.  Properly
;   handle places where qsq > 1.  Implemented NOZPHASE.
; 11/09/2009: DGG. Initial implementation of simultaneous projection
;   to multiple z planes for more efficient volumetric
;   reconstructions.  Eliminate matrix reorganizations.
; 06/10/2010: DGG. Documentation fixes.  Added COMPILE_OPT.
; 07/27/2010: DGG. Added HANNING.
; 10/20/2010: DGG: Subtract 1 from (normalized) hologram
;   rather than expecting user to do subtraction.
;
; Copyright (c) 2006-2010 Sanghyuk Lee and David G. Grier
;
; UPDATES:
;    The most recent version of this program may be obtained from
;    http://physics.nyu.edu/grierlab/software.html
; 
; LICENSE:
;    This program is free software; you can redistribute it and/or
;    modify it under the terms of the GNU General Public License as
;    published by the Free Software Foundation; either version 2 of the
;    License, or (at your option) any later version.
;
;    This program is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;    General Public License for more details.
;
;    You should have received a copy of the GNU General Public License
;    along with this program; if not, write to the Free Software
;    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;    02111-1307 USA
;
;    If the Internet and WWW are still functional when you are using
;    this, you should be able to access the GPL here: 
;    http://www.gnu.org/copyleft/gpl.html
;-
function rayleighsommerfeld, a, z, $
                             lambda = lambda, $ ; wavelength of light
                             mpp = mpp, $       ; micrometers per pixel
                             nozphase = nozphase, $
                             hanning = hanning 

COMPILE_OPT IDL2

; hologram dimensions
sz = size(a)
ndim = sz[0]
if ndim gt 2 then message, "requires two-dimensional hologram"
nx = float(sz[1])
if ndim eq 1 then $
  ny = 1. $
else $
  ny = float(sz[2])

; volumetric slices
nz = n_elements(z)              ; number of z planes

; parameters
if n_elements(lambda) ne 1 then $
   lambda = 0.632               ; HeNe laser
if n_elements(mpp) ne 1 then $
   mpp = 0.135                  ; Nikon rig

ik = complex(0., 2.*!pi * mpp / lambda)       ; wavenumber in radians/pixel

; phase factor for Rayleigh-Sommerfeld propagator in Fourier space
; Refs. [2] and [3]
qx = findgen(nx) - nx/2.
qx *= lambda / (nx * mpp)
qsq = qx^2

if ndim eq 2 then begin
   qy = findgen(ny) - ny/2.    
   qy *= lambda / (ny * mpp)
   qsq = qsq # replicate(1.,ny) + replicate(1.,nx) # qy^2
endif

qfactor = sqrt(complex(1. - qsq))

if ~keyword_set(nozphase) then $
   qfactor -= 1.

if keyword_set(hanning) then $
   qfactor *= hanning(nx, ny)

E = fft(complex(a - 1), -1, /center) ; Fourier transform of input field
res = complexarr(nx,ny,nz)
for j = 0, nz-1 do begin
   if z[j] gt 0 then $
      Hqz = exp(-ik * float(z[j]) * conj(qfactor)) $
   else $
      Hqz = exp(-ik * float(z[j]) * qfactor)

; Convolve field with Rayleigh-Sommerfeld propagator
   thisE = E * Hqz                            ; Convolve with propagator
   thisE = fft(thisE, 1, /center, /overwrite) ; Transform back to real space

   res[*,*,j] = thisE
endfor

return, res

end