;+
; NAME:
;	DOE_ADD
;
; PURPOSE:
;	For combining phase profiles, modulo 2pi.
;
; CATEGORY:
;	Hologram/DOE Generator
;
; CALLING SEQUENCE:
;	result = doe_add(phase1, phase2) -- *_NOT_YET_TESTED_*
;
; INPUTS:
;	phase1: a phase profile created specifically for the particular SLM at hand
;
;	phase2: a second phase profile, to be added to the first
;
; OPTIONAL INPUTS (KEYWORD PARAMETERS):
;	circ: Turning on this keyword means that the "corners" of
;		the SLM will NOT be used (adding that light to the zeroth
;		order), to avoid breaking the circular symmetry of the beam
;
;	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)
;
; OUTPUTS:
;	result: a real array giving a "combined" phase plate
;
; SIDE EFFECTS:
;
; REQUIRES:
;	Calibration of the phase levels of the device used, for the desired 
;	operating conditions. (Is what we call 2pi calibrated to BE 2pi ?)
;	An appropriate look-up table (LUT) would allow the images generated 
;	by this program to be used directly by the device.
;
; RESTRICTIONS:
;	Extremely large phase gradients will be "aliased" due to the
;	discretization of the device (the number of phase levels and the size
;	of the pixels).
;
; MODIFICATION HISTORY:
;	Created by Gabe Spalding, June, 2010
;	Modified Oct. 2021 for use with SLMs with rectangular active areas
;-

FUNCTION DOE_add, phase1, phase2, circ=circ, nbit=nbit,$
	xpix=xpix, ypix=ypix, maxGSL=maxGSL

 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

result = phase1 + phase2

; I want to restrict my phase modulations
; to the range [ 0, maxGSL ]
; Below is my own modulo function, written to avoid using mod function
; on REAL arguments (See Gumley, Practical IDL Programming, p. 42)
modulo= maxGSL
print, 'Pre-Modulo phase range = [', min(result), ',', max(result),']'
result = result - (modulo * long(result/modulo))

IF keyword_set(circ) THEN BEGIN
 sz = size(result, /dimensions)
 xpix = sz[0]
 xoffset = long(0.5 * xpix)
 yoffset = long(0.5 * ypix)
 rho_max = SQRT( (xoffset)^2 + (yoffset)^2 )
 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 result(i,j) = 0
  ENDFOR
 ENDFOR
ENDIF

; At this point, the phase values of "result" extends 
; over the following range
print, 'Final phase range = [', min(result), ',', max(result),']'

return, result
END