FORT_structure_overview.hp42s

LBL "FORT"
@ Finds the minimum of the FORT potential and optionally calculates the
@ trap frequencies.
@
@
@ input:
@  MEH: Matrix of H-FORT unit column vectors 
@       (waist 1 direction, waist 2 direction, propagation direction)
@  MEV: Matrix of V-FORT unit colums vectors 
@       (waist 1 direction, waist 2 direction, propagation direction)
@
@ output:
@  TRAPF: 3x1 column vector with trap frequencies (Hz)
@  UMIN: Potential at trap minimum (in µK)
@
@ internally used (global storage):
@  PWRH: Power of horizontal FORT (used by IFORT)
@  PWRV: Power of vertical FORT (used by IFORT)
@  MASS: atomic mass in amu
@  VR: 3x1 column vector of position
@  VRMIN: 3x1 column vector of trap minimum position
@  VRMAX: 3x1 column vector of trap minimum position
@  H: 3x3 Hessian matrix at trap minimum
@  ATOM: name of species (e.g. 168Er)
@  FLAG 01: Calculate Jacobian or Hessian
@  FLAG 02: Search for maximum instead of minimum
@  FLAG 03: Error in searching extremum (no trap?)
@  FLAG 04: Recalculation of trap geometry needed if set
@  FLAG 05: Minimum of trap potential not known if set
@  FLAG 06: Plot to printer instead of to LCD if set

	LBL A
	LBL a
	@ Program menus
	@
	LBL 20
	@ Dummy loop for program menu
	LBL 11
	@ H-FORT power setting
	LBL 12
	@ V-FORT power setting
	
	LBL B
	@ Search for trap minimum and trap maximum.
	@ Calls LBL C two times, without and with FLAG 02 set.
	@ Also calculates final Hessian and potentials,
	@ and save all calculation results
	@
	@ output:
	@  H: 3x3 Hessian matrix at trap minimum
	@  UMIN: Potential at trap minimum (in µK)
	@  UMAX: Potential at trap maximum (in µK)
	@  VRMIN: 3x1 column vector of trap minimum position
	@  VRMAX: 3x1 column vector of trap maximum position
	LBL 21
	@ Finish trap calculation by calling function to calculate
	@ trap frequencies and print out result.

	LBL C
	@ Newton Algorithm.
	@ Determines either minimum (FLAG 02 false) or maximum (FLAG 02 true).
	@ Before returning, FLAG 03 will be set if extremum was not found.
	@
	@  N: counter for number of Newton iterations
	LBL 02
	@ Loop for initial minimum/maximum search.
	@ If FLAG 02 is false, searches for minimum.
	@ If FLAG 02 is true, searches for maximum.
	@ Starts at current position VR and goes down in steps of 2 µm
	@ until beyond minimum or maximum.
	@
	@ internally used (local storage):
	@  VAL: current potential value
	@  M: counter of search interations (30 iterations max)
	LBL 01
	@ Newton Algorithm main loop
	@
	@ internally used (local storage):
	@  J: Jacobian 3x1 column vector
	@  DELTA: 3x1 column vector, current Newton step size
	LBL 03
	@ Could not find trap extremum. Cancel Newton.
	@
	@ output:
	@  FLAG 03: set to true
	LBL 07
	@ Newton algorithm converged
	LBL 08
	@ Newton algorithm did not converge within N steps
	LBL 06
	@ Display trap depth and gravitational sag
	@
	@ input:
	@  UMIN: Potential at trap minimum (in µK)
	@  UMAX: Potential at trap maximum (in µK)
	@  VRMIN: 3x1 column vector of trap minimum position
	@
	@ output:
	@  ALPHA: informational message
	LBL 22
	@ helper label to jump directly to grav. sag printing in
	@ case the trap depth is unknown

	LBL D
	@ Calculates the eigenvectors of a 3x3 matrix.
	@ Follows the implementation of 
	@ https://en.wikipedia.org/wiki/Eigenvalue_algorithm
	@
	@ input:
	@  ST X: 3x3 matrix (in this case the Hessian matrix)
	@  MASS: atomic mass in amu
	@
	@ output:
	@  TRAPF: 3x1 column vector with trap frequencies (Hz)
	@
	@ internally used (local storage):
	@  A: copy of 3x3 matrix for which eigenvalues are computed
	@  Q: trace(H)/3
	@  P1, P2, P, PHI: intermediate values
	@  B: 3x3 intermediate matrix
	@  EIG1, EIG2, EIG3: eigenvalues
	
	LBL E
	@ Plot trap potential in vertical direction.
	@ If FLAG 04 is clear plots along the line connecting trap
	@ minimum and maximum, otherwise through zero along Z.
	@ If FLAG 06 is set plots to printer instead of to LCD.
	@
	@ internally used (local storage):
	@  VDIR: 3x1 unit column vector of plot direction
	LBL 30
	@ Plot direction is determined, start plot
	@ If FLAG 06 is set jump to LBL d to plot to printer instead.
	@
	@ input:
	@  VDIR: 3x1 unit column vector of plot direction
	@
	@ internally used:
	@  VR: 3x1 vector of current plot position
	LBL 31
	@ main plotting loop
	LBL 32
	@ Determine vertical position of pixel from ST X function value
	LBL 33
	@ Draw a vertical line
	@
	@ input:
	@  ST X: position (in spatial units)
	LBL 34
	@ Draw a horizontal line
	@
	@ input:
	@  ST X: position (in µK units)
	LBL c
	@ Calculate potential at given position
	@ The position will be translated into a 3D point along
	@ the plot axis.
	@
	@ input:
	@  ST X: position (in m)
	@
	@ output:
	@  ST X: potential in µK
	@
	@ internally used (global storage):
	@  XMIN, XMAX: minimal and maximal vertical value
	@  YMIN, YMAX: minimal and maximal vertical value
	
	LBL 41/42/43
	@ Species selection for Er/Yb/Li
	@
	@ output:
	@  POLARI: polarizability in atomic units (positive for attractive pot.)
	@  MASS: atomic mass in amu
	
	LBL 51
	@ Entry point for summary printout to HP82240A printer
	@ Starts by executing PRON and printing 
	LBL 52
	@ Loop to print out all three trap frequencies

	LBL d
	@ Start of actual function plotting to printer.
	@ This is called from LBL E (our standard plotting function)
	@ when FLAG 06 is set)
	@
	@ internally used (local storage):
	@  _GrMod: graphics mode before chaning to low-res for printing
	@  _DX: x-axis resolution in m/pixel
	@  _XPOS: current position of plot in m
	LBL 53
	@ Outer loop over the single screens, built and printed one-by-one
	@
	@ internally used (local storage):
	@  _YSCR: current vetical position on screen in pixel (loop variable)
	LBL 54
	@ Inner loop over vertical positions (_YSCR) on current screen
	LBL 55
	@ Print final y-axis (via XEQ 57 and XEQ 60 with FLAG 01 set)
	LBL 56
	@ Function drawing callback.
	@ Takes position _XPOS, evaluates the potential (via XEQ c), converts
	@ into screen coordinates and draws a PIXEL at that position
	LBL 57
	@ Draws a horizontal line on the LCD at position ST X.
	@ After the routine the value of ST X is unchanged.
	LBL 58
	@ Draws one pixel at the left and one pixel at the right LCD boundary.
	@ After the routine the value of ST X is unchanged.
	LBL 59
	@ Draws x-axis tic marks and tic labels.
	@ Does not try to restore the stack before returning
	LBL 60
	@ Draws y-axis tic marks at vertical position ST X. 
	@ With FLAG 01 cleared initial marks are drawn below ST X.
	@ With FLAG 01 set final marks are drawn above ST X.
	@ FLAG 01 is cleared before returning.
	@ Implemented with FUNC 00 so the stack is unchanged.
	LBL 61
	@ Prints lower, upper and central y-axis value labels.
	@ The labels are positioned at the bottom of the screen so that the
	@ actual plot with box and tics start at the top of the next screen.
	LBL 62
	@ Calculates estimated width of current ALPHA string in pixel.
	@ For now every character is assumed to take 4 pixels.
	@
	@ input:
	@  ALPHA: input string
	@
	@ output:
	@  ST X: width in pixel
	LBL 63
	@ Finishes plotting to printer by restoring initical graphics mode
	
	LBL 99
	@ Exits the program cleanly.
	@ Clears FLAG 04, FLAG 05 and the programmable menu.

@ Polarizabilities used:
@ Er: 166 a.u. at 1064.5 nm from Tab. I of
@     163 a.u. at 1570.0 nm from Tab. I of
@     J. H. Becher et al., Phys. Rev. A 97, 012509 (2018).
@ Yb: 162.5 a.u. at 1064 nm from Tab. B.2 of
@     150.4 at 1550.0 nm by repetition of calculation Eq. (B.2)
@     C. D. Herold, Dissertation, University of Maryland, 2014.
@ Li: 269 a.u. by interpolation to 1070 nm of data in Tab. I of
@     201.3 a.u. from data in Tab. I of
@     M. S. Safronova et al., Phys. Rev. A 73, 022505 (2006).
@
@ Note that the calculation of Herold (PhD) for Rb is in close agreement with
@ the results given by Safronoa (PRA), giving credibility to Herold's data.

LBL "MROT"
@ Creates 3x3 rotation matrix Mrot = Mrot_z * Mrot_y * Mrot_x, where Mrot_i
@ is a rotation about axis i.
@
@ input: 
@  ST X, ST Y, ST Z: rotation angles along X, Y, Z
@ output: 
@  MROT: rotation matrix
@  ST X: rotation matrix (copy of MROT)
@
@ interally used:
@  REG 01: initial ST X (rotation angle about x)
@  REG 02: initial ST Y (rotation angle about y)
@  REG 03: initial ST Z (rotation angle about z)


LBL "WSTZ"
@ Calculates beam radius of Gaussian beam at distance Z from waist position.
@
@ input:
@  ST x: distance Z from waist
@  LAMBDA: wavelength
@  WAIST: waist size
@
@ output:
@  ST X: beam radius at Z


LBL "IGAUSS"
@ Calculates the intensity of an elliptical Gauss beam oriented in an 
@ arbitrary direction at and arbitrary position in space. The beam is
@ assumed to go through the origin of the coordinate system and to also
@ have its waist there.
@
@ input:
@  ME: 3x3 matrix containing unit vectors with beam geometry.
@      column 1: direction of 1st waist
@      column 2: direction of 2nd waist
@      column 3: direction of beam propagation
@      (Effectively this is just MROT)
@  WST1: size of 1st waist
@  WST2: size of 2nd waist
@  LAMBDA: wavelength
@  VR: 3x1 vector with (x, y, z) laboratory position where intensity is wanted
@  POWER: Total power of beam
@
@ output:
@  ST X: intensity at position VR
@
@ internally used:
@  REG 04: position in 1st waist direction
@  REG 05: position in 2nd waist direction
@  REG 06: position along beam axis
@  REG 07: radius of beam in 1st waist direction
@  REG 08: radius of beam in 2nd waist direction
@  REG 10-12: used during plotting
	LBL A
	@ Get (rotated) unit vector that defines beam coordinate system
	@ input: 
	@  ST X: axis selection (1, 2 or 3)
	@  ME: 3x3 matrix containing unit vectors as column vectors
	@ return: 
	@  ST X: 3x1 vector that is the column vector of MROT at position ST X
	LBL "B"
	@ Get position of laboratory frame vector for one axis in the beam frame
	@ input:
	@  ST X: beam frame axis selection (1, 2 or 3)
	@  VR: 3x1 column vector with (x, y, z) laboratory frame position
	@  ME: 3x3 matrix containing unit vectors as column vectors
	@ output:
	@  ST X: position of vector along axis that was selected by ST X

LBL "IFORT"
@ Calculates the intensity of the crossed FORT at the given position and
@ with the given powers.
@
@ For the ErYbLi experiment we assume the following laboratory coordinate frame:
@  X: Zeeman slower direction
@  Y: H-FORT direction
@  Z: V-FORT direction
@
@ input:
@  VR: 3x1 vector with (x, y, z) laboratory position where intensity is wanted
@  PWRH: Power of horizontal FORT
@  PWRV: Power of vertical FORT
@  MEH: Maxtrix of H-FORT unit column vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@  MEV: Maxtrix of V-FORT unit colums vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@
@ output:
@  ST X: intensity at position VR
@
@ internally used:
@  REG 04-08: in the execution of IGAUSS
@  REG 09: intermediate result from intensity calculation
@  LAMBDA: wavelength of beams
@  POWER: power of current FORT beam
@  WST1, WST2: waists of beams

LBL "UFORT"
@ Calculates the potential of the crossed FORT including gravity.
@ Currently polarization and mass are hard-coded for Erbium.
@
@ input:
@  VR: 3x1 vector with (x, y, z) laboratory position where potential is wanted
@  POLARI: Polarizability in atomic units (positive for attractive potential)
@  MASS: Atomic mass in amu
@  PWRH: Power of horizontal FORT (used by IFORT)
@  PWRV: Power of vertical FORT (used by IFORT)
@  MEH: Matrix of H-FORT unit column vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@  MEV: Matrix of V-FORT unit column vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@
@ output:
@  ST X: potential at VF (in Joule)
@
@ internally used:
@  REG 04-08: in the execution of IGAUSS
@  REG 09: intermediate result from intensity calculation
@  LAMBDA, WST1, WST2: in the execution of IFORT

LBL "UDERIV"
@ Calculates either the 1st derivatives (Jacobian) of the potential (FLAG 01 is 
@ false) or the 2nd derivates (Hessian) (FLAG 01 is true) at position VR.
@ 
@ input:
@  VR: 3x1 vector with (x, y, z) laboratory position where potential is wanted
@  PWRH: Power of horizontal FORT (used by IFORT)
@  PWRV: Power of vertical FORT (used by IFORT)
@  MEH: Matrix of H-FORT unit column vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@  MEV: Matrix of V-FORT unit colums vectors
@       (waist 1 direction, waist 2 direction, propagation direction)
@  FLAG 01: 1st or 2nd derivative selection
@
@ output:
@  ST X: derivative (either 3x1 or 3x3 matrix)
@
@ internally used (global storage):
@  REG 04-08: in the execution of IGAUSS
@  REG 09: in the execution of IFORT
@  REG 10: step size in numerical derivative
@  LAMBDA, WST1, WST2: in the execution of IFORT
@ internally used (local storage):
@  COOR: coordinate for which to calculate the derivative (1, 2 or 3)
@  VDER1: 3x1 column vector that holds the 1st derivative
	LBL A
	@ Calculate a single 1st derivative and store result into the 
	@ correct position of VDER1.
	@
	@ input:
	@  COOR: coordinate for which to calculate the derivative (1, 2 or 3)
	@
	@ output:
	@  VDER1: 3x1 column vector of 1st derivatives with updated element
	@
	@ internally used (local storage):
	@  VAL1: function value at positive step
	@  VAL2: function value at negative step
	LBL B
	@ Calculate the complete Hessian matrix.
	@
	@ input:
	@  -- nothing additional --
	@
	@ output:
	@  ST X: 2nd derivative 3x3 matrix
	@
	@ internally used (local storage):
	@  MDER2: 3x3 matrix that hold the 2nd derivative
	@  C2: coordinate of the 2nd derivative
	LBL C
	@ Calculate a single element of the Hessian matrix and store
	@ the value into the correct position of MDER2.
	@
	@ input:
	@  COOR: coordinate of the 'inner' 1st derivate
	@  ST X: coordinate of the 2nd derivative
	@
	@ output:
	@  MDER2: 3x3 matrix of 2nd derivatives with updated element
	@
	@ internally used (local storage):
	@  VAL1: function value at positive step
	@  VAL2: function value at negative step