function fastphase, p, subsample = subsample, undersample = undersample
;+
; NAME:
;      FASTPHASE
;
; PURPOSE:
;      Calculates the phase hologram encoding a desired trapping
;      pattern by superposing fields as fast as possible.
;
; CATEGORY:
;      Computed holography
;
; CALLING SEQUENCE:
;      phi = fastphase(p)
;
; INPUTS:
;      p: [2,npts] or [3,npts] array of (x,y) coordinates of points in the plane,
;         relative to the center of the field of view.
;         Coordinates are measured in arbitrary pixel units, unless
;         the CAL keyword is set.
;
; KEYWORDS:
;      subsample: If set, calculate half-sized hologram for half requested 
;          trap spacing, then tile to create full-sized hologram.
;          Much faster.  Lower spatial resolution.
;
;      undersample: If set, calculate half-sized hologram and then
;          rebin to full size.  Fastest of all.  Crummy holograms.
;           
; OUTPUTS:
;      phi: phase pattern encoding pattern of traps described by
;           P with values ranging from [0,2 pi].
;
; RESTRICTIONS:
;      Can be very memory intensive for large numbers of points.
;
; PROCEDURE:
;      Initial estimate is obtained by superposing the fields of the
;      specified beams.
;
; NOTES:
;      If we calibrate the SLM's phase transfer function, then we
;      should be able to pass the lookup table to this function, and
;      return a hologram of indices into the look-up table.
;
; MODIFICATION HISTORY:
; Created by David G. Grier, New York University, 12/12/2004.
; 1/19/2005: DGG.  Implemented 3D.
; 1/29/2005: DGG.  Major code clean-up leading to improved speed.
;                  Implemented SUBSAMPLE for major speedup.
; 9/29/2006: DGG.  Modified to use calibration constants in
;     HOLO_COMMON.  Removed DIM and CAL keywords.
;
; Copyright (c) 2006 David G. Grier
; 
; 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 shold be able to access the GPL here: 
;    http://www.gnu.org/copyleft/gpl.html
;-

common holo_common, cal

w = cal.doe_w
h = cal.doe_h

; calibration constants
xc = cal.xc                     ; center of phase mask on SLM
yc = cal.yc
xfac = cal.xscale               ; scale factor for square pixels
rfac = cal.scale                ; projection scale factor
thetac = cal.theta              ; orientation


sp = size(p)
ndim = sp[1] ; number of dimensions
npts = sp[2] ; number of points

 subsample = keyword_set(subsample) and ndim eq 2
 if subsample then begin
     w /= 2
     h /= 2
     rfac /= 2
 endif

 undersample = keyword_set(undersample)
 if undersample then begin
     w /= 2
     h /= 2
 endif

twopi = 2.D * !dpi

; trap locations in trapping plane
; rotate points to account for relative SLM-CCD orientation
qx = p[0, *] * cos(thetac) + p[1, *] * sin(thetac)
qy = p[1, *] * cos(thetac) - p[0, *] * sin(thetac)
; wavevectors associated with trap position (times i)
ikx = dcomplex(0., (twopi/w) * reform(qx))
iky = dcomplex(0., (twopi/h) * reform(qy))

if ndim gt 2 then $
   ikz = dcomplex(0., twopi * cal.zfactor * reform(p[2,*]))

; coordinates in SLM plane (row vectors)
x = xfac * rfac * dindgen(w)
y = rfac * dindgen(h)
if ndim gt 2 then begin
    xsq = (x - w/2.D - xc)^2
    ysq = (y - h/2.D - yc)^2
endif

iphase = dcomplex(0.D, twopi * randomu(seed, npts)) ; relative phases

psi = complexarr(w, h)
for n = 0, npts-1 do begin
    ; i \vec{k}_n \cdot \vec{r}
    ikxx = ikx[n] * x + iphase[n]
    ikyy = iky[n] * y
    if ndim gt 2 then begin
        ikxx += ikz[n] * xsq
        ikyy += ikz[n] * ysq
    endif
    ex = exp(ikxx)
    ey = exp(ikyy)
    ; \psi += exp( i \phi_n(\vec{r}) )
    psi += ex # ey
endfor

phi = atan(imaginary(psi), double(psi)) + !dpi

 if subsample then $
   phi = [[phi, phi], [phi, phi]]

 if undersample then $
   phi = rebin(phi,2.*w,2.*h)

return, phi

end