How to Write User-Defined Sources and Scatter Functions in Fortran

User-defined objects, scatter functions, sources etc are dll programs that communicate with Zemax. In this article I show how to write such programs in Fortran. Sample code is provided.
 
Jeff Casey
12/10/2007
Extensions

Introduction

Zemax has powerful abilities to let you write your own user-defined surfaces, objects, scattering functions, diffraction functions etc as a compiled function -known as a dll- and to link it into Zemax. Such user-defined functions work just like built-in functions. Zemax is supplied with many examples written in C.

In this note, I’ll show some specific techniques which were successful in getting Zemax to accept user-defined objects written in the Lahey/Fujitsu Fortran 95 Compiler. None of these UDOs have been widely tested in an exhaustive variety of platforms, so there may yet be bugs that haven’t been found – nevertheless, they have proven useful, and may be so to others. (I’d appreciate feedback and/or bugs reports to casey@rockfieldresearch.com).

The Lahey/Fujitsu Compiler is a popular and powerful Fortran compiler, with full 32 bit flat memory and many structured enhancements that allow easy mixing of old legacy code with more readable functionality. For those beyond a certain age, Fortran remains an easier language to think and generate code in. This article is aimed at  those people. The examples were compiled with Lahey/Fujitsu release 7.10.02.  All of these user-defined functions were written for non-sequential optics applications.

Compiling Fortran Code to a dll

User-defined objects are compiled code and have the file extension .dll. This extension shows that the code is a 'dynamic link library', which an executable program (.exe) or other dll can call and execute.

Once such code is written, it must be compiled with the proper switches. An appropriate command  line, executed from a command shell, is:

$ LF95 <sourcecode>.f90 –win –ml msvc –dll

This will generate the “source.dll” file, which can then be moved to the appropriate directory, generally under  {Zemaxroot}\ Objects \ DLL \ .

Multi-Threaded Code

Zemax takes full advantage of multiple CPUs if they are available. The extension software must be written to  be compatible with multiple threads. For example, some of my early attempts used the compiler-supplied function for random number generation. This works fine if Zemax is forced to only use a single CPU, but  generates erratic failures if multiple CPUs are used in a multi-core machine. (These are not frequent failures,  but will occur at low probability, depending on the likelihood of collision by multiple processes – as an example, this glitch usually only occurred if  > 106 rays were traced on a dual-core system).

The Lahey/Fujitsu compiler does not provide switches for multi-threaded compatibility, but it does provide such switches for generating re-entrant code. Since parameters are kept on the stack with re-entrant calls, this similarly solves the compatibility problem. Your code, therefore, should be completely self-contained, and internal routines should be specified with the RECURSIVE flag.

User-Defined Sources

Here is an example of the “SOURCE DLL”. This allows you to generate an initial cloud of rays which  match the peculiarities of your particular problem, and is particularly handy when you are mating Zemax to modeling codes for unusual radiating objects. These compiled DLL files will go into:  {Zemaxroot}\ Objects \ DLL \ Sources for Zemax to find and use.

Below, I include and annotate a complete source dll , which generates a bundle of rays from a single point in a planar fan, constrained between two angles. This is a very nice debugging tool.

function UserSourceDefinition (data)

implicit none

integer*4, dll_export :: UserSourceDefinition


! for the call to UserSourceDefinition, the data is stored as follows:

! data(0) = total number of values in the passed data array

! data(1) = x (to be returned)

! data(2) = y (to be returned)

! data(3) = z (to be returned)

! data(4) = l (to be returned)

! data(5) = m (to be returned)

! data(6) = n (to be returned)

! data(7) = relative intensity (to be returned)

! data(20) = wavelength in um

! data(21) = mm per unit length (1.0 for mm, 25.4 for inches, etc.)

! data(22) = random number seed

! data(30) = parameter 1 from user input

! data(31) = parameter 2 from user input

! etc.... up to data(maxdata), where maxdata = int(data(0))

! This DLL will compute and return data 1 through 7 given other data.

! Return value is normally 0 on success, -1 on failure.


real*8 data(0:*), pi, dr, theta, thetam, thetap, rx

integer*4 seed


parameter (pi = 3.14159265358979323846d0)

parameter (dr = pi/180.d0)


seed = int(data(22))

! start random sequence with Zemax's 4-byte random integer


thetam = data(30) ! one limit to fan

thetap = data(31) ! other limit to fan


data(1) = 0.

data(2) = 0.

data(3) = 0. ! all rays originate at a single point


if (thetam .eq. thetap) then

theta = thetam

else

theta = thetam + rx(seed) * (thetap-thetam)

end if

theta = theta * dr ! convert to radians


data(4) = sin(theta)

data(5) = 0.

data(6) = cos(theta)

! spread rays out in x-z plane (y-direction cosine always zero)

data(7) = 1.0 ! relative intensity is constant

UserSourceDefinition = 0 ! return with null error code

return

end function


The first function, “UserSourceDefinition” is the main function called by Zemax to generate rays. Since this is  declared as type “dll_export”, the parameter passing and name mangling are shuffled to be consistent with the  Zemax constraints (as programmed in C). (Note – UserSourceDefinition is case sensitive, unlike most FORTRAN variables, since it is an exported name.)

The comment header in UserSourceDefinition shows exactly what data is available on code entry, and what data must be generated on exit. Note that the xyz position (data(1) through data(3)) is relative to the placement of the source in Zemax. Since this is a point source, all ray origins are at (0,0,0)....but this would not be so for a distributed source. The ray direction (direction cosines data(4) through data(6)) is also relative to the orientation of the source in Zemax. Here, we have defined a fan in the x-z plane – any other orientation is accommodated by rotations in the placement of the source.

Next, we found that we needed to call an internal routine to generate random numbers. Since we wish to avoid  system calls that are not multi-thread compatible, we use our own. In this case, our needs are not very demanding, so a very simple linear congruential generator is applicable. This one is from Knuth, and taken out  of “Numerical Recipes”. For copyright purposes, I have suppressed the initial values for I1 and I2. (There are  numerous references to more sophisticated random number generators that make fascinating reading – but we will not delve into those here.) 

real*8 recursive function rx(seed)

implicit none

integer*4 i1,i2,seed

parameter (i1 = <suppressed for copyright reasons>)

parameter (i2 = <suppressed fro copyright raesons>)

seed = i1*seed + i2

rx = dble(seed)/2.**32 + .5

return

end function


Last, we have another exported call/entry, UserParamNames. This is called whenever Zemax wants to refresh the headers for the parameter tables, and describes what the various parameters are in short text. (This sort of  coding is where FORTRAN does not cope gracefully with text handling, but can nevertheless be overcome). The data array element 0 contains the pointer to the parameter number for which an ASCIIZ string should be returned in data (1:*).

function UserParamNames (data)

implicit none

integer*4, dll_export :: UserParamNames

byte data(0:*)

integer i

i = int(data(0))

if (i .lt. 0) i = i + 256 ! convert from unsigned

if (i .eq. 1) then

data(0) = ichar("t")

data(1) = ichar("h")

data(2) = ichar("e")

data(3) = ichar("t")

data(4) = ichar("a")

data(5) = ichar("-")

data(6) = 0

else if (i .eq. 2) then

data(0) = ichar("t")

data(1) = ichar("h")

data(2) = ichar("e")

data(3) = ichar("t")

data(4) = ichar("a")

data(5) = ichar("+")

data(6) = 0

else

data(0) = ichar("-")

data(1) = 0

end if

UserParamNames = 0

return

end function


Beyond this, the specific algorithmic needs of your source will determine the remainder of your code, which will replace the simple source calculations in UserSourceDefinition above. Remember that all internal routines must be marked recursive (as the function “rx” above was), and also that error checking within Zemax is not possible for user code – the most likely result of a programming error is hanging Zemax.

For many sources, especially distributed sources with a particular profile, you will need to properly invert the probability distribution to give you the correct distribution of origins (in space) or directions (in solid angle) from a simple uniform deviate (i.e. random number between 0 and 1). You will often need special function solutions, numerical root finders, or other complex calculations for this purpose. This is an algorithmic problem, however, and not specific to FORTRAN, and will not be amplified here.

User-Defined Objects

Another useful user element is the “USER OBJECT”, which allows you to generate an arbitrary surface or solid object. The  compiled DLL file must be placed in the following directory on the root disk for Zemax to find: {Zemaxroot}\ Objects \ DLL \ UserObjects.

Generation of a User Object is very similar to generation of a source. The main routine is necessarily more complicated, as it must serve several key functions:
  • it must generate the rendering and solid model information (as a surface of adjacent triangles),
  • it must annotate these triangular areas by surface/solid, grating, and CSG flags,
  • it must calculate the details of an intersection when passed the information for a ray,
  • it must calculate polarization data, when asked.

The User Object code has two Zemax calls (exported, case sensitive): UserObjectDefinition and UserParamNames. The header and comments for UserObjectDefinition are:

function UserObjectDefinition (data, tri_list)

implicit none

integer*4, dll_export :: UserObjectDefinition

real*8 data(0,*), tri_list (0:*)


! for the call to UserObjectDefinition, the data is stored as follows:

!      data(1) = calling mode (see below)

!      data(2-4) = x,y,z position for origin of ray traced

!      data(5-7) = cx,cy,cz direction cosines of ray traced

!      data(8) = exact code (see below)

!      data(10...) = returned values (depends on mode)

!      data(100) = number of parameters passed to DLL

!      data(101...) = first parameter in list, next, ...


! Return value is normally 0 on success, -1 on failure.


! if called with data(1) = 0: compute rendering and solid model info

!      return total # of triangular faces for rendering in data(10)

!      return solid/shell flag in data(11) (solid=1, surface=0)

!      return grating info in data(12) (not=0, yes, return CSG #)

!      return coating code in data(20) (default=0, compute coating info locally=1)


! if called with data(1) = 1: generate rendering mesh

!      return filled array of tri_list, real*8 array 10*N long

!          (10 values per triangle)

!      sequence for each triangle: (x1,y1,z1,x2,y2,z2,x3,y3,z3,code)

!      code is 6 digit integer EECCRV, decimally parsed:

!      V = visibility: 0=all legs, 1=1-2 inv, 2=2-3 inv,

!          4=3-1 inv...sum for multiples

!      for rectangles, make common leg of two triangles invisible

!      R = 0 for refractive, 1 for reflective

!      CC = coat scatter group...only 00 through 03 supported at this time

!      EE = exact code: 0 = triangle is exact solution, other values are used to

!      tag which internal algorithm may be used to calculate intersection

!      ALSO, must return # of triangular faces again in data(10)


! if called with data(1) = 2: find intersection and normal vector to surface

!      data(2,3,4) = ray origin

!      data(5,6,7) = ray direction cosines

!      data(8) = exact code, tells which algorithm to use

!      find the intersection by determining path length of ray, S

!      intersection at (x+Scx, y+Scy, z+Scz)

!      return S in data(10)

!      return normal vector to surface at intersection point in data(11,12,13)

!      if ray misses, return errorcode -1


! if called with data(1) = 3: do coating calculation

!      data(2,3,4) = ray origin

!      data(5,6,7) = ray direction cosines

!      data(8) = CSG (coating/scatter group)

!      data(9,10,11) = normal vector at surface

!      data(12,13) = current and next indices of refraction

!      data(14) = cosine of norm incidence angle

!      data(15) = wavelength in um

!      data(20) = zero flag (set to 1 if routine is to calculate coating,

!      otherwise uses object settings)

!      return data(20) = 1


! S-polarization amplitudes in data(21,22,23,24) (refl_r, refl_i, xmit_r, xmit_i)

! P-polarization amplitudes in data(31,32,33,34) (refl_r, refl_i, xmit_r, xmit_i)

! (don't forget that amplitude factor is sqrt of intensity factor)

! return error code -1 if this DLL doesn't do coating calcs

! if called with data(1) = 4: do initialization with safe data

! return values into data(101....) as default values


The coding of UserObjectDefinition is dominated by algorithmic concerns. Most of the above comments and specifications were taken from the Zemax manual and from other examples on the Zemax knowledge base. Not all options have been tested to date, but a variety of surfaces have been generated and found to work quite well.

It is most certain that your coding of a surface will require a root finder, a normal vector generator, and perhaps several other functions. Remember to code all of these as re-entrant (recursive), and make no system calls from within your function.

The coding of UserParamNames is very similar to that used for Sources. One example is shown mostly complete below:

function UserParamNames (data)

implicit none

integer*4, dll_export :: UserParamNames

character*10 data

integer i1, i2, i3

i1 = ichar(data(1:1))

i2 = ichar(data(2:2))

i3 = ichar(data(3:3))


data = char(0)

if (i2 .eq. 0) then 

if (i1 .eq. 48) then ! 0 

data = "All Faces" // char(0) 

else if (i1 .eq. 49) then ! 1 

data = "param 1" // char(0) 

else if (i1 .eq. 50) then ! 2 

data = "param 2" // char(0)

.......

end if


else if ((i3 .eq. 0) .and. (i1 .eq. 49)) then 

if (i2 .eq. 48) then ! 10 

data = "param 10" // char(0) 

else if (i2 .eq. 49) then ! 11 

data = "param 11" // char(0)

........

end if

end if

UserParamNames = 0

return


end function


For UserParamNames, we have parsed the input data a bit differently this time (also clumsily). In this case, we  can parse out more than 10 parameters, plus we must handle a pointer less than or equal to zero. For the USER  OBJECT types, Zemax will call UserParamNames with a zero or negative pointer to get the names of the CSG  (coating scatter group). For this example, we have are using a surface in 3D, which will be described by a single  CSG, thus we pass “All Faces” back for a pointer of 0.

User-Defined Surface Scattering Functions

One last example of FORTRAN code callable by Zemax is the Scattering type object. These compiled DLL files go into the directory: {Zemaxroot}\ Objects \ DLL \ SurfaceScatter.

The scattering objects allow a custom calculation for the scattering properties of a surface. Like the other examples, this also has two callable entries: UserScatterDefinition and UserParamNames. The latter, UserParamNames is identical to the call for sources, and serves only to name the user definable parameters. The entry for UserScatterDefinition is similar to the other examples. The header and comments are given below:

function UserScatterDefinition (data)

! for the call to UserScatterDefinition, the data is stored as follows:

!      data(0) = total number of values in the passed data array

!      data(1) = x position of specular ray

!      data(2) = y position of specular ray

!      data(3) = z position of specular ray

!      data(4) = l (x direction cosine of specular ray in, of scattered ray out)

!      data(5) = m (y direction cosine of specular ray in, of scattered ray out)

!      data(6) = n (z direction cosine of specular ray in, of scattered ray out)

!      data(7) = lxn (x direction cosine of normal vector to surface)

!      data(8) = myn (y direction cosine of normal vector to surface)

!      data(9) = nzn (z direction cosine of normal vector to surface)

!      data(10) = 0 initially, return 1 if DLL scatters ray,

!          2 if full polarization given

!      data(11) = mm per unit length (1.0 for mm, 25.4 for inches,

!          10.0 for cm, 1000 for m)

!      data(12) = relative energy (to be computed by dll and returned)

!      data(16) = random seed

!      data(17) = wavelength in um

!      data(20) = incident Ex real

!      data(21) = incident Ex imag

!      data(22) = incident Ey real

!      data(23) = incident Ey imag

!      data(24) = incident Ez real

!      data(25) = incident Ez imag

!      data(40) = output Ex real (to be returned if data(10) set to 2)

!      data(41) = output Ex imag (to be returned if data(10) set to 2)

!      data(42) = output Ey real (to be returned if data(10) set to 2)

!      data(43) = output Ey imag (to be returned if data(10) set to 2)

!      data(44) = output Ez real (to be returned if data(10) set to 2)

!      data(45) = output Ez imag (to be returned if data(10) set to 2)

!      data(51) = input parameter 1 from user

!      data(52) = input parameter 2 from user

!      etc... up to data(maxdata), where maxdata = int(data(0))

! This DLL will compute and return data(i) for i=4,5,6,10,12,40-45 given other data.

! Return value is normally 0 on success, -1 on failure.


implicit none

integer*4, dll_export :: UserScatterDefinition

real*8 data(0:*)


In general, the input data will far exceed the needs of the calculation – any particular scattering algorithm will likely be fairly restrictive, and hence use a small subset of this data. The case we have used is to model the grazing incidence reflection coefficients for far UV / soft X-ray optics using simple polynomial fits. In that particular case, the rest of the code consists of:

real*8 data(0:*), dot, pi, p1, p2, theta, thetamax

parameter (pi = 3.14159265358979323846)

parameter (thetamax = <constant>)

dot = data(4)*data(7) + data(5)*data(8) + data(6)*data(9)

theta = abs(asin(dot)) ! this is grazing angle

if (theta .gt. thetamax) then

data(12) = 0.0 ! cutoff

else

p1 = <polyfit numerator >

p2 = <polyfit denominator >

data(12) = p1/p2

end if

data(10) = 1.0 ! flag to show that we scattered

UserScatterDefinition = 0 ! return with null error code

return

end function


In short, we have only calculated the grazing angle (from the dot product of the specular ray and the normal ray), and then reduced the energy of the ray according to the reflection coefficient model.

Acknowledgements

These techniques were developed while consulting for Energetiq Inc., of Woburn, MA. I thank my colleagues at Energetiq for helpful and enjoyable brainstorming and support.

I have greatly appreciated email conversations with Tony Richards and Mark Nicholson, who helped tremendously in getting past the first nest of problems.