One of the example scattering DLLs provided with the ZEMAX installation is Gaussian_XY. The scattering distribution modeled by this DLL is a Gaussian about the projected specular ray. If p and q represent the coordinate axes on the scattering surface (i.e. perpendicular to the surface normal), then this distribution is represented analytically by:

P(p,q) = (4/(p*s_p*s_q))*exp[-((p/s_p)² + (q/s_q)²)]

As the probability is represented in Cartesian coordinates in this case, the corresponding integral is slightly simplified:

I = ∫∫ P(p,q) dp dq = Erf(p/s_p)*Erf(q/s_q)

where Erf is the Error Function. Since P is separable, it is straightforward to assign values for the scatter direction based on random numbers:

RN1 = Erf(p/s_p); p = s_p*Erf^-1(RN1)
RN2 = Erf(q/s_q); q = s_q*Erf^-1(RN2)

Unfortunately the Error Function and its inverse are not built-in numeric functions in C. However, an alternate solution for this case exists:

p = s_p*√(-ln(RV1))*cos(2p*RV2)
q = s_q*√(-ln(RV1))*sin(2p*RV2)

This solution is based on the Box-Muller method for generating random numbers with a Gaussian distribution, as described in Chapter 7.2 of Numerical Recipes in C, Second Edition¹.

The source code for this DLL (Gaussian_XY.c) is located in the \ZEMAX\Objects\DLL\SurfaceScatter\ directory. The main code starts with the Windows APIENTRY UserScatterDefinition. In the first part of this code, variables are defined, the random number generator is initialized, and the DLL reports to ZEMAX that a ray has been scattered and what the energy of that ray is (relative to the input energy available for scattering):

/* this DLL models a surface that has a Gaussian scattering
   distribution about the projected specular ray */

Next, the coordinate axes p and q are determined such that they are perpendicular to the surface normal n, and such that the projected specular ray lies in the p-n plane:

     /* find vectors p,q perpendicular to n, such that the
        specular ray lies along the p-n plane */
     px = sx; // initialize p with values from the specular ray
     py = sy;
     pz = sz;

     MAG = px*nx + py*ny + pz*nz;

     if (MAG > 0.9999999)
           {
           /* specular is very close to the normal, 
              any p will do */
           if (fabs(nz) < 0.9)
                   {
                   px = 0.0;
                   py = 0.0;
                   pz = 1.0;
                   }
           else
                   {
                   px = 1.0;
                   py = 0.0;
                   pz = 0.0;
                   }
           }

     /* this creates q normal to n and p */
     CrossProduct(nx, ny, nz, px, py, pz, &qx, &qy, &qz);
     Normalize(&qx, &qy, &qz);
     /* this creates p normal to both q and n */
     CrossProduct(nx, ny, nz, qx, qy, qz, &px, &py, &pz);
     //Normalize(&px, &py, &pz); // not needed since n and q
                                 // are orthonormal already

     /* the specular ray, when projected onto the surface,
        lies along p */
     bo = px*sx + py*sy + pz*sz;

After the data for s_p and s_q are read in from the user input array, the scattered ray direction is determined using two random numbers:

     get_another_scattered_ray:;

     MAG = 2.0;
     while (MAG > 1.0)
           {
           random_number = (double)rand()/(double)RAND_MAX;
           if (random_number == 0.0) goto get_another_scattered_ray;
           random_val1 = random_number;
           random_number = (double)rand()/(double)RAND_MAX;
           random_val2 = random_number;
           bp = sig_x*sqrt(-1.0*log(random_val1))*cos(2.0*PI*random_val2);  // Gaussian distribution along p
           bq = sig_y*sqrt(-1.0*log(random_val1))*sin(2.0*PI*random_val2);  // Gaussian distribution along q
           tx = bo*px + bp*px + bq*qx;  // t is the scattered ray
                                        // vector projected onto
                                        // the surface
           ty = bo*py + bp*py + bq*qy;
           tz = bo*pz + bp*pz + bq*qz;
           MAG = tx*tx + ty*ty + tz*tz; // validate that scattered
                                        // ray vector is within
                                        // unit circle
           }
     MAG = sqrt(1.0 - MAG);

     /* this is the scattered ray */
     data[4] = tx + MAG*nx;
     data[5] = ty + MAG*ny;
     data[6] = tz + MAG*nz;

The random numbers are used to determine the scattering direction with respect to the projected specular ray according to the Box-Muller equations. The direction of the scattered vector in the p-q plane is then calculated using those values and the known direction of the projected specular ray in this plane. A check is performed to make sure that the scattered vector in the projected plane has a magnitude less than unity; if it does not, two new random numbers are drawn and the process is repeated. The actual trajectory of scattered ray in (x,y,z) space is finally determined by accounting for the component of scattering along the surface normal, which is calculated from the constraint that the magnitude of the scatter vector is unity. Once this calculation is complete, a check is performed to ensure that the ray did not scatter “beyond the horizon”, i.e. more than 90 degrees from the surface normal:

     /*
     NOW, make sure we didn't scatter in such a way that the ray 
     is no more than 90 degrees from normal!
     */

     if (nx*data[4] + ny*data[5] + nz*data[6] < 0.0)
           {
           goto get_another_scattered_ray;
           }

     return 0;
     }

If this condition is violated, the process is repeated from the beginning by drawing two new random numbers; otherwise the scattering calculation is complete for the given input ray.

A simple example (Gaussian_XY.zmx) which illustrates the use of this DLL is provided in the .ZIP file located at the end of this article. For s_p = s_q = 0.5, the radiant intensity distribution looks like: