;+
; NAME:
;	LG
;
; PURPOSE:
;	For creating an approximated, generalized (p not necessarily zero!) 
;	Laguerre-Gaussian mode out of a fundamental Gaussian input beam --
;	and, if desired, to shift it away from the undiffracted central
;	spot of the diffractive optic
;
; CATEGORY:
;	Hologram/DOE Generator
;
; CALLING SEQUENCE:
;	LGmode = lg(l, p)
;
; INPUTS:
;	l: the multiple of 2pi associated with traversal about the center 
;	 ("azimuthal number")
;
;	p: the radial index for the LG mode ("radial number")
;
; OPTIONAL INPUTS (KEYWORD PARAMETERS):
;	xpix: number of pixels on a horizontal side, 
;		default setting is 1920 (for JDC EDK SLM)
;		(Set to 1280 for our Cambridge Correlators SLM)
;		(Set to 1024 or 768 for our Hamamatsu SLM)
;
;	ypix: number of pixels on a horizontal side, 
;		default setting is 1080 (for JDC EDK SLM)
;		(Set to 1024 for our Cambridge Correlators SLM)
;		(Set to 768 for our Hamamatsu SLM)
;
;	noscl: unless this is set to 1, the output is subjected to
;	     the modulo procedure and SCALED from 2pi to a grayscale
; 
;	maxGSL: the grayscale level corresponding to 2pi of phase, 
;		default setting is 2^nbit-1
;
;	nbit: provides 2^nbit levels of phase, 
;		default setting is 8 (good for most SLMs)
;
;	rad: the radius of a circular discontinuity,
;		whenever the radial mode number is non-zero
;		(must fit on the SLM or DOE, for the beam to sample it)
;
;	blaze: the number of 2pi blazings across the field 
;	     (we've been using 120 as a default)
;	  -- this is effectively adding a PRISM to refract the beam 
;	     away from the central spot
;	  -- increasing 'blaze' is like increasing the angle of the exit
;	     face of the prism
;
;	noblaze: if no blazing is desired, this keyword must be set 
;		(see RESTRICTIONS below)
;
; OUTPUTS:
;	LGmode: a xpix by ypix real array giving a spiral phase plate, 
;		(with, for non-zero p, radial phase shifts added)
;
; SIDE EFFECTS:
;	For l = 0, p = 0, this creates a blazed diffraction grating 
;		(the DOE equivalent of a prism)
;
; REQUIRES:
;	Calibration of the phase levels of the SLM used, for the desired 
;	operating conditions, in terms of the 'greyscale' levels used here.
;	An appropriate look-up table (LUT) would allow the images generated 
;	by this program to be used directly by the SLM.
;
; RESTRICTIONS:
;	The program does not directly recognize BLAZE = 0, as a
;	consequence of how keyword_set works. My 'workaround' was to 
;	add the argument /NOBLAZE whenever no blazing is desired.
;	(Note that '/nob' will work as an abbreviation for '/NOBLAZE', since
;	in general a slash preceding a keyword sets its value to one, and
;	keywords can be abbreviated to their shortest UNIQUE length. 
;	For example, a keyword called XSTYLE could be abbreviated to XST)
;
;	Extremely large phase gradients will be "aliased" due to the
;	discretization of the SLM (the number of phase levels and the size
;	of the pixels).
;
; MODIFICATION HISTORY:
;	Created by Gabe Spalding, October, 2002
;	Modified June, 2010 to ensure that the modulo function deals well
;	with negative l-values.
;	Modified April 2016 to allow for an unscaled output that is
;	not (yet) subjected to the modulo procedure
;	Modified Sept. 2021 for use with SLMs with rectangular active areas
;-

FUNCTION lg,l, p, blaze=blaze, noblaze=noblaze, rad=rad, nbit=nbit, $
	xpix=xpix, ypix=ypix, maxGSL=maxGSL, circ=circ, noscl=noscl

 IF NOT keyword_set(blaze) THEN blaze = 120
 IF keyword_set(noblaze) THEN blaze = 0
 IF NOT keyword_set(rad) THEN rad = 100
 IF NOT keyword_set(nbit) THEN nbit = 8
 IF NOT keyword_set(maxGSL) THEN maxGSL = 2^nbit-1
 IF NOT keyword_set(xpix) THEN xpix = 1920
 IF NOT keyword_set(ypix) THEN ypix = 1080
 IF NOT keyword_set(noscl) THEN noscl = 0

 LGmode=fltarr(xpix,ypix)
 xoffset = long(0.5 * xpix)
 yoffset = long(0.5 * ypix)
 FOR i=0,xpix-1 DO BEGIN
   FOR j=0,ypix-1 DO BEGIN

     ; First, create a phase blazing to shift everything away
     ; from the undiffracted spot:
     ; (Whenever no blazing is desired, add the argument /noblaze.)
     LGmode(i,j) =  i * (2*!pi) * blaze/(xpix-1)

     ; Next, create a p=0 LG mode 
     ; via addition of a sprial phase profile:
     ; The azimuthal angle is given by
     phi = ATAN(j-yoffset, i-xoffset)+!pi
       ; A term containing the azimuthal angle multiplied by
       ; the azimuthal number creates the spiral phase profile:
    LGmode(i,j) =  LGmode(i,j) + l * phi

       ; Now, find how far the present pixel is from the center
       ; of the diffractive optical element (DOE)
    rho = SQRT((i-xoffset)^2.0 + (j-yoffset)^2.0)

       ; Use this radius as a criteria for adding further phase shifts that 
       ; depend upon the radial number:
    IF p GT 0 THEN BEGIN
	    FOR pindex = 1, p DO BEGIN
  	      ; For non-zero radial number, the DOE will contain p circles,
  	      ; with sizes given by:
          circsize = pindex*rad
          ; Beyond each of those circles, we add a pi phase shift,
          IF rho GE circsize THEN LGmode(i,j) = (LGmode(i,j) + !pi)
      ENDFOR
    ENDIF

  ENDFOR
 ENDFOR

modulo = 2*!pi
;print, 'Pre-Modulo phase range = [', min(LGmode), ',', max(LGmode),']'
WHILE ( min(LGmode) LT 0 ) DO LGmode = LGmode + !pi

IF keyword_set(noscl) THEN BEGIN
  print, 'Output not scaled'
ENDIF ELSE BEGIN  
  LGmode = LGmode - (modulo * long(LGMode/modulo))
; At this point, the range of "LGmode" extends over [ 0, +2pi ] ??
;print, 'Post-modulo phase range = [', min(LGmode), ',', max(LGmode),']'

   ; Finally, scale everything from units of phase to units of 
   ; 'image level'
LGmode = LGmode*((maxGSL)/(2*!pi))
   ; The extra parentheses above improve efficiency: 
   ; otherwise the precedence of multiplication over division would mean 
   ; that the array would be operated on twice.
print, 'Greyscale range:', min(LGmode), max(LGmode)
ENDELSE

IF keyword_set(circ) THEN BEGIN
  rho_max = xoffset < yoffset
  print, 'rho_max:', rho_max
  FOR i=0,xpix-1 DO BEGIN
   FOR j=0,ypix-1 DO BEGIN
   ; The radial coordinate, from the center of the DOE, is given by:
     rho = SQRT((i-xoffset)^2 + (j-yoffset)^2)
     IF rho GT rho_max THEN LGmode(i,j) = 0
   ENDFOR
  ENDFOR
ENDIF

return,LGmode
END