;+
; NAME:
;	BLAZE
;
; PURPOSE:
;	To create the DOE equivalent of a prism (namely, a blazed 
; 	diffraction grating), for shifting (refracting) a beam
;	away from the undiffracted central spot. This program should
;	provide easy adjustment of both the degree of beam deflection
;	and the direction of deflection.
;
; For purposes of maintaining one up-to-date version, other programs that 
; will use BLAZE.pro (such as BESSEL.pro, etc.) should merely call 
; this program (and should list it in their headers under 'REQUIRES:')
;
; CATEGORY:
;	Hologram/DOE Generator
;
; CALLING SEQUENCE:
;	prism = blaze()
;
; INPUTS:
;
; OPTIONAL INPUTS (KEYWORD PARAMETERS, which can be abbreviated
;		 to the shortest unique form):
;
;	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)
;
;	maxGSL: the grayscale level corresponding to 2pi of phase, 
;		default setting is 2^nbit-1
;
;	fullblaze: this is the number of 2pi blazings 
;		across a HORIZONTAL field (we've been using 120)
;		-- this is effectively adding a prism to refract the beam
;		 away from the central spot
;		-- increasing 'fullblaze' is like increasing the angle of
;		the exit face of the prism
;
;	dir: the direction of the blazing, an angle (in degrees)
;		(this might be needed when used in conjunction with
;		some spatial filters that block out light along the
;		directions associated with the pixelation noise of the SLM)
;
;	nbit: provides 2^nbit levels of phase, 
;		default setting is 8 (appropriate for newer SLMs)
;
; OUTPUTS:
;	PRISM: a xpix by ypix real array giving a wedge phase plate
;
; SIDE EFFECTS:
;
; 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:
;	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, July, 2010, to bypass scaling to 2pi in favor of
;	directly scaling to bytscl (old form left, but commented out)
;	Modified Sept. 2021 for use with SLMs with rectangular active areas
;-

FUNCTION blaze,fullblaze=fullblaze, dir=dir, nbit=nbit, $
	xpix=xpix, ypix=ypix, maxGSL=maxGSL

 IF NOT keyword_set(fullblaze) THEN fullblaze = 120
 IF NOT keyword_set(dir) THEN dir = 0
 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

     ; First, convert the direction specified by dir to radians:
  dir_rad = dir*(2*!pi)/360.
  theta_corner = ATAN(ypix/xpix)

     ; The maximum distance along dir is:
  IF dir LE theta_corner THEN $
  	max_dist = (xpix - 1.)/COS(dir_rad) $
  ELSE $
  	max_dist = (ypix - 1.)/SIN(dir_rad)

 prism=fltarr(xpix, ypix)

 FOR i=0,xpix-1 DO BEGIN
  FOR j=0,ypix-1 DO BEGIN

       ; The distance of this pixel from the CORNER of the DOE is:
    r = SQRT(i^2.0 + j^2.0)
    
       ; The angular coordinate of this pixel with respect to
       ; the bottom axis of the DOE is
       ; (given that the origin is in the CORNER):
    alpha = ATAN(j, i)
       ; The RELATIVE angle, then, with respect to the rotated coord system
       ; specified by dir is:
    beta = dir_rad - alpha
    
    
       ; Next, create an index showing how far this pixel is
       ; along the direction specified by dir:
    dist = r * COS(beta)
; It MUST be in this step that ROUNDING ERRORS creep in!
; TEST this by constructing an array that can be PLOTTED:
; prism[i,j] = dist
; The problem, for dir=45 on a square SLM,
; occurs exactly halfway to max_dist

       ; Now, create a phase blazing to shift everything away
       ; from the undiffracted spot:
;    prism[i,j] =  (2*!pi) * fullblaze * dist/(max_dist)
    prism[i,j] =  (2^nbit-1) * fullblaze * dist/(max_dist)

  ENDFOR
 ENDFOR

; print, 'Pre-Modulo phase range = [', min(prism), ',', max(prism),']'
WHILE ( min(prism) LT 0 ) DO prism = prism + (2^nbit-1)
; print, 'Prism Phase Range = [',min(prism), max(prism),']'
;modulo = 2*!pi
modulo = (2^nbit-1.)
prism = prism - (modulo * long(prism/modulo))
; print, 'Post-Modulo phase range = [', min(prism), ',', max(prism),']'

   ; Finally, scale everything from units of phase to 
   ; units of 'image level'
;prism = prism*((2^nbit-1)/(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.
; imgui, prism
return,prism
END