;+
; NAME:
;	BESSEL
;
; PURPOSE:
;	To create a the DOE equivalent of an axicon (a conical lens)
;		-- to create approximated Bessel beams 
;			(or circle beams in far field)
;		out of a fundamental Gaussian or an LG mode input beam --
;	If this is the final DOE, it should be rounded.
;
; CATEGORY:
;	Hologram/DOE Generator
;
; CALLING SEQUENCE:
;	BessMode = bessel(ph_max)
;
; INPUTS:
;	ph_max: The maximum phase shift (in units of 2pi)
;	**WHAT IS THE OPTIMAL SETTING FOR ph_max? (in terms of alignment 
;	and effect)** Obviously, as ph_max goes to zero, the axicon effect 
;	disappears and we recover the incident beam, and yet ...
;	as ph_max decreases, the propagation distance for the Bessel-
;	like pattern goes up and less of the incident intensity goes 
;	into the rings (which can be good if what you want is to guide 
;	something along the core). However, as ph_max decreases, the width
;	of the core goes up, and so whatever is being guided is less-well
;	localized in the direction transverse to the beam. Note, however,
;	that Bessel beam tweezers do not make the axicon conjugate to the back
;	aperture of the objective (or final) lens; rather, the final lens is 
;	merely one of two lenses involved in a telescope that allows one to
;	independently scale down the width of the core.
;
;	Setting ph_max to, say, 6 gives a pretty image, but in order to
;	determine an appropriate setting, one must consider the wavelength,
;	the beam waist at the SLM, and the properties of the telescope used
;	to introduce the Bessel Beam into the sample cell. In order to make 
;	comparisons to Bessel Beams produced by glass axicons, I made up a
;	formula for converting from ph_max to an effective axicon angle:
;	gamma = ATAN(2*!pi*lambda*ph_max*(1-n)/(n*beamwaistataxicon))
;
; 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)
;
;	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
;
;	irad: a radius, inside of which everything will be set to zero
;		(creating an axicon with a hole in the middle = Loony Toons logo)
;		Units of irad are not in pixels, but in increments that
;		ensure no sharp discontinuities in phase are introduced 
;		(beyond those associated with our lack of calibration
;		and the fact that our '2pi' is therefore not identical to
;		'zero')
;
;	orad: a radius, outside of which everything will be set to zero
;		(creating an axicon circular boundary conditions)
;		Units of orad are the same as those of irad.
;
;	nbit: provides 2^nbit levels of phase, 
;		default setting is 8 (good for most SLMs)
;
; OUTPUTS:
;	BessMode: a npix x npix real array giving a conical phase plate
;
; SIDE EFFECTS:
;
; REQUIRES:
;
;	To add a blazing, run either BLAZE.pro or LG.pro (with l=0, p=0).
;
;	LG.pro (if the input is to be an LG beam)
;	The current 'workaround' is to make a DOE image using 
;	bessel(ph_max), and another using lg(l, p), and then to create a 
;	linear superposition:
;		DOE = (DOE1 + DOE2)MOD(2^nbit-1)
;
;	Also, blazing should be added if the beam is to be deflected away 
;	from the central, undiffracted spot.
;
;	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:
;	ph_max is limited since 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 by GS Nov, 2002 to avoid using mod function on 
;		real arguments (See Gumley, Practical IDL Programming, p. 42)
;	Modified by GS Feb, 2005 to improve commentary on optimal choice of 
;	ph_max
;	Modified May 2023 for use with SLMs with rectangular active areas
;
; KNOWN ISSUES:
;	At ph_max = 1550 there is a cutoff I don't understand, as if orad becomes
;	suddenly too small. When nbit and npix are set for the Boulder 512x512, this
;       cutoff shows up at around 1050.
;-


FUNCTION bessel,ph_max, irad=irad, orad=orad, nbit=nbit, $
	xpix=xpix, ypix=ypix, maxGSL=maxGSL, circ=circ, noscl=noscl

 IF NOT keyword_set(nbit) THEN nbit = 8
 IF NOT keyword_set(npix) THEN npix = 768
 IF NOT keyword_set(maxGSL) THEN maxGSL = 2^nbit
 IF NOT keyword_set(xpix) THEN xpix = 1920
 IF NOT keyword_set(ypix) THEN ypix = 1080
 IF NOT keyword_set(noscl) THEN noscl = 0
 IF NOT keyword_set(irad) THEN irad = 0
 IF NOT keyword_set(orad) THEN orad = 2*npix

 BessMode=fltarr(xpix,ypix)
 xoffset = long(0.5 * xpix)
 yoffset = long(0.5 * ypix)
 rho_max = SQRT((xoffset + 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)

	; The conical phase plate is created as follows
     BessMode(i,j) = (2*!pi)* ph_max * (rho_max - rho)/rho_max

     IF rho LT irad*(rho_max/ph_max) THEN BessMode(i,j) = 0
     IF rho GT orad*(rho_max/ph_max) THEN BessMode(i,j) = 0

  ENDFOR
 ENDFOR

modulo = 2*!pi
BessMode = BessMode - (modulo * long(BessMode/modulo))

   ; Optionally, we scale everything from units of phase to units of 
   ; 'greyscale level'
IF keyword_set(noscl) THEN BEGIN
  print, 'Output not scaled'
ENDIF ELSE BEGIN  
  BessMode = long(BessMode*((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(BessMode), max(BessMode)
ENDELSE

return,BessMode
END