Zemax POP macro: normalized plot and calculation of beam size

Zemax Physical Optics Propagation (POP) generates cross-sectional beam profile.  Two functions are lacking that are useful to me:
1. An irradiance-normalized plot. (The absolute "watts-per-sq-millimeter" value does not matter to me almost all the time).
2. Zemax gives "second-moment" beam diameter but again, almost all the time I just need a simple 1/e^2 diameter regardless of perfect Gaussian beam or not. (The second-moment beam diameter definition heavily weighs the tails or outer wings of the intensity profile. Therefore beams with side-lobes will have second-moment width substantially larger than their central lobe widths. [1])

So I wrote a macro file to do these two tasks:

!-----------------------------------------------------------
!FIRST, SAVE SETTINGS IN POP WINDOW. 
!The "Show As" in "Display" tab must be set as "Cross X (or Y)".
!----------------------------------------------------------- 
!Input POP parameters saved in POP window:
!-----------------------------------------------------------
sampling = 512
zoom = 8

AA$ = $TEMPFILENAME()
GETTEXTFILE AA$, POP    # Import POP text data file
OPEN AA$

!-----------------------------------------------------------
! Create a loop to read POP data
!-----------------------------------------------------------
invalid_row = 13        # The first 13 rows are text, not data
For i,1,invalid_row,1
    READ a,b
Next i

Grid_size = sampling/zoom        # POP grid size (=Sampling/zoom)
DECLARE xx, DOUBLE, 1, Grid_size+1                   
DECLARE Irrad, DOUBLE, 1, Grid_size+1
For i,1,Grid_size+1,1
    READ a,b
    xx(i) = a
    Irrad(i) = b
    print xx(i), ", ", Irrad(i)
Next i
CLOSE

!-----------------------------------------------------------
! Find peak
!-----------------------------------------------------------
Max_i = 1
For j,2,Grid_size+1,1
!    PRINT "j=",j,"xx=",xx(j),"Irrad=",Irrad(j)
    IF ( Irrad(j)>Irrad(j-1) ) & ( Irrad(j)>Irrad(Max_i) ) THEN Max_i = j
Next j

!-----------------------------------------------------------
! Normalize lineshape
!-----------------------------------------------------------
DECLARE Irrad_norm, DOUBLE, 1, Grid_size+1
For j,1,Grid_size+1,1
    Irrad_norm(j) = Irrad(j)/Irrad(Max_i)
    PRINT "j=",j," xx=",xx(j)," Irrad_norm=",Irrad_norm(j)
Next j
!-----------------------------------------------------------
! Find 1/e^2 radius at extreme left (linear interpolation)
!-----------------------------------------------------------
Radius_Li = 1
For j,1,Max_i,1
    IF ( Irrad_norm(j)<=0.135 ) & ( Irrad_norm(j+1)>=0.135 ) THEN GOTO 1    #once the condition is met, jump out loop.
Next j
LABEL 1
Radius_Li = j
R_L1 = xx(Radius_Li)
R_L2 = xx(Radius_Li+1)
Irrad_L1 = Irrad_norm(Radius_Li)
Irrad_L2 = Irrad_norm(Radius_Li+1)
R_L = (0.135-Irrad_L1)*(R_L1-R_L2)/(Irrad_L1-Irrad_L2) + R_L1
PRINT "At ", R_L, " Irradiance = 0.135"
!-----------------------------------------------------------
! Find 1/e^2 radius at extreme right (linear interpolation)
!-----------------------------------------------------------
Radius_Ri = 1
For j,Max_i,Grid_size,1
    IF ( Irrad_norm(j)>=0.135 ) & ( Irrad_norm(j+1)<=0.135 ) THEN Radius_Ri = j+1
Next j
R_R1 = xx(Radius_Ri)
R_R2 = xx(Radius_Ri-1)
Irrad_R1 = Irrad_norm(Radius_Ri)
Irrad_R2 = Irrad_norm(Radius_Ri-1)
R_R = (0.135-Irrad_R1)*(R_R1-R_R2)/(Irrad_R1-Irrad_R2) + R_R1
PRINT "At ", R_R, " Irradiance = 0.135"
!-----------------------------------------------------------
! Calculate 1/e^2 beam diameter
!-----------------------------------------------------------
dia = R_R - R_L
PRINT "Beam Diameter = ", dia
!-----------------------------------------------------------
! Plot normalized lineshape
!-----------------------------------------------------------
PLOT NEW
PLOT DATA, xx,Irrad_norm,Grid_size+1,0,0,0
PLOT BANNER, "Normalized POP beam profile"
PLOT RANGEX, xx(1), xx(Grid_size+1)
PLOT RANGEY, 0, 1
PLOT COMM1, $DATE()
PLOT COMM2, $FILEPATH()
PLOT COMM3, "1/e^2 Gaussian Diameter = ", $STR(dia)
PLOT GO

Save the above text as a "xxx.zpl" file in Zemax's ZPL folder (the folder path can be viewed/changed in File -> Preferences). To run this macro, press F9, choose this file and click "Execute".

Here is an example. The POP window generated by Zemax shows an aberrated beam profile:


My macro generates a plot shown below.  First, it is normalized.  Second, the 1/e^2 beam diameter is given.



Compare three beam sizes:
1. POP pilot beam radius (paraxial 1/e^2) = 0.0070
2. POP "second moment" beam radius = 0.0071
3. Simple 1/e^2 radius calculated by my POP macro = 0.0024
The three sizes are different. Most often I just need know the 3rd value.

Next: It would be nice to also generate a Gaussian-fit curve.  I will perhaps write a MATLab-linked code to do this task.

Reference
[1] A.E.Siegman, "How to measure laser beam quality". PDF available on web by googling it.

Zemax as Gaussian beam calculator

This is version 2 of the last note on this subject.

Example 1: Given an input Gaussian beam of waist radius W01, a lens of focal length f, what is the output beam's waist radius W02 and waist location z2?

To quickly calculate the paraxial result, a simple Zemax model can be created as follows. In this example, W01 = 0.5 mm, z1 = 10 mm, f = 50 mm and wavelength = 640 nm.

The aperture and field settings are irrelevant to Gaussian calculation; they are set for making the optical layout look good:

Pressing [Ctrl]+[b] opens the Paraxial Gaussian Beam Data window. Right-click mouse opens the setting window:

So the waist radius (1/e^2 value) is set to 0.5 mm at the surface 1. Press OK the paraxial Gaussian beam results are given:

Notice that the image plane, the geometrical focus, is not exactly at the Gaussian beam waist. So a simple Merit Function can be used to move the image plane to the waist:

Optimize it and the image plane is adjusted to be at the waist. As the result, the Gaussian beam focus is about 66 um in front of the geometrical focus:

Example 2: We still have the same input Gaussian beam W01 = 0.5 mm and the lens f = 50 mm, but we want the focused beam waist radius to be 5 um. To do this, another lens is needed to expand the input beam first; say we have a f = -10 mm lens in hand. Where does this lens need to be? (d = ?)

To solve this, simply insert the lens f1 = -10 mm into the Zemax model:

Add in a GBPW operand into the Merit Function to calculate Gaussian waist. The target value is 5 um:

After optimization, Zemax gives the correct thickness of lens 1 to be 80.74 mm and the Paraxial Gaussian result shows that the beam is focused to 5 um radius:

Note that these are the first-order, paraxial calculations. Real lenses have aberrations and usually cannot give exactly the same results and will usually yield strong side-lobes. The next step should be to replace with the real lens models and use POP for actual beam evaluation.