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.