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.