#include <windows.h>
#include <math.h>
#include <string.h>
#include "usersurf.h"

/*
Written by James Sutter
Sept 21, 2001

This DLL models an arbitrary number of lens elements in a lens array.
The individual lenses are even aspheres.

The user provides the number of elements in x and y, the size in x and y of
each element, and the radius, conic, asphere coefficients, and glass.

The asphere sag in each element is given by:

Z = (c*r*r) / (1+(1-((1+k)*c*c*r*r))^ 1/2 ) + a2*(r)^2 + a4*(r)^4 + ... a8*(r)^16

Note the terms a2 ... a16 have units of length to the -1, -3 ... -15 power.

This surface breaks up the beam into numerous separate
beams, and so most ZEMAX features will fail to work with this surface.

However, the spot diagrams, image analysis, etc, all work okay.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD);
int GetCellCenter(int nx, int ny, double wx, double wy, double x, double y, double *cx, double *cy);

/* a generic Snells law refraction routine */
int Refract(double thisn, double nextn, double *l, double *m, double *n, double ln, double mn, double nn);

BOOL WINAPI DllMain (HANDLE hInst, ULONG ul_reason_for_call, LPVOID lpReserved)
	{
   return TRUE;
   }

/* this DLL models a lens array surface type */

int  __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, nx, ny, error, miss_flag;
   double alpha, power, t;
   double wx, wy, cx, cy, x, y, z;
   int loop;
   double r2, tp, dz, sag, mm, mx, my;
   double a[9];
   
	switch(FD->type)
   	{
      case 0:
      	/* ZEMAX is requesting general information about the surface */
         switch(FD->numb)
         	{
            case 0:
            	/* ZEMAX wants to know the name of the surface */
		         /* do not exceed 12 characters */
		         strcpy(UD->string,"Even Array");
               break;
            case 1:
            	/* ZEMAX wants to know if this surface is rotationally symmetric */
               /* it is not, so return a null string */
            	UD->string[0] = '\0';
               break;
            case 2:
            	/* ZEMAX wants to know if this surface is a gradient index media */
               /* it is not, so return a null string */
            	UD->string[0] = '\0';
            	break;
            }
         break;
      case 1:
      	/* ZEMAX is requesting the names of the parameter columns */
         /* the value FD->numb will indicate which value ZEMAX wants. */
         /* they are all "Unused" for this surface type */
         /* returning a null string indicates that the parameter is unused. */
         switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "Number X");
               break;
            case 2:
            	strcpy(UD->string, "Number Y");
               break;
            case 3:
            	strcpy(UD->string, "Width X");
               break;
            case 4:
            	strcpy(UD->string, "Height Y");
               break;
            default:
            	UD->string[0] = '\0';
            	break;
            }
      	break;
      case 2:
      	/* ZEMAX is requesting the names of the extra data columns */
         /* the value FD->numb will indicate which value ZEMAX wants. */
         /* they are all "Unused" for this surface type */
         /* returning a null string indicates that the extradata value is unused. */
			switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "2nd Order Term");
               break;
            case 2:
            	strcpy(UD->string, "4th Order Term");
               break;
            case 3:
            	strcpy(UD->string, "6th Order Term");
               break;
            case 4:
            	strcpy(UD->string, "8th Order Term");
               break;
            case 5:
            	strcpy(UD->string, "10th Order Term");
               break;
            case 6:
            	strcpy(UD->string, "12th Order Term");
               break;
            case 7:
            	strcpy(UD->string, "14th Order Term");
               break;
            case 8:
            	strcpy(UD->string, "16th Order Term");
               break;
            default:
	           	UD->string[0] = '\0';
            	break;
				}
			break;
      case 3:
      	/* ZEMAX wants to know the sag of the surface */
         /* if there is an alternate sag, return it as well */
         /* otherwise, set the alternate sag identical to the sag */
         /* The sag is sag1, alternate is sag2. */

			/* aspheric terms*/
         for (i = 1; i<=8; i++)
	        	{
				a[i] = FD->xdata[i]; // even terms
	  	      }

         UD->sag1 = 0.0;
         UD->sag2 = 0.0;

			/* if a plane, just return */
			if (FD->cv == 0) return(0);

         /* figure out the center coordinates of which "cell" we are in */
         nx = (int) FD->param[1];
         ny = (int) FD->param[2];
         wx = FD->param[3];
         wy = FD->param[4];
         x = UD->x;
         y = UD->y;

         /* make sure nx and ny are both odd, otherwise the chief ray is a problem... */
         if (!nx&1)nx++;
         if (!ny&1)ny++;
         if (wx <= 0.0 || wy <= 0.0) return(-1);
         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &cx, &cy);

         if (error) return(0);

         /* offset the coordinates */
         x -= cx;
         y -= cy;

         r2 = x*x + y*y;
         alpha = 1 - (1+FD->k)*FD->cv*FD->cv*r2;
         if (alpha < 0)
         	{
            /* lens is not big enough to fill the cell, radius is too short */
            /* assume it is plane between lenses */
            alpha = 0.0;
            }
         UD->sag1 = (FD->cv*r2)/(1 + sqrt(alpha));

         /* now the aspheric terms */
         /* this is numerically inefficient, but easy to read */
         if (a[1] != 0.0)
         	{
				UD->sag1 += a[1] * r2;
            }
         if (a[2] != 0.0)
         	{
				UD->sag1 += a[2] * r2 * r2;
            }
         if (a[3] != 0.0)
         	{
				UD->sag1 += a[3] * r2 * r2 * r2;
            }
         if (a[4] != 0.0)
         	{
				UD->sag1 += a[4] * r2 * r2 * r2 * r2;
            }
         if (a[5] != 0.0)
         	{
				UD->sag1 += a[5] * r2 * r2 * r2 * r2 * r2;
            }
         if (a[6] != 0.0)
         	{
				UD->sag1 += a[6] * r2 * r2 * r2 * r2 * r2 * r2;
            }
         if (a[7] != 0.0)
         	{
				UD->sag1 += a[7] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
            }
         if (a[8] != 0.0)
         	{
				UD->sag1 += a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2;
            }
			/* forget supporting a hyper hemisphere! */
         UD->sag2 = UD->sag1;
      	break;
      case 4:
      	/* ZEMAX wants a paraxial ray trace to this surface */
         /* x, y, z, and the optical path are unaffected, at least for this surface type */
         /* for paraxial ray tracing, the return z coordinate should always be zero. */
         /* paraxial surfaces are always planes with the following normals */

         /* for a lens array, only consider the single lens on axis */
         UD->ln =  0.0;
         UD->mn =  0.0;
         UD->nn = -1.0;
         power = (FD->n2 - FD->n1)*FD->cv;
         if ((UD->n) != 0.0)
         	{
            (UD->l) = (UD->l)/(UD->n);
            (UD->m) = (UD->m)/(UD->n);

            (UD->l) = (FD->n1*(UD->l) - (UD->x)*power)/(FD->n2);
            (UD->m) = (FD->n1*(UD->m) - (UD->y)*power)/(FD->n2);

            /* normalize */
            (UD->n) = sqrt(1/(1 + (UD->l)*(UD->l) + (UD->m)*(UD->m) ) );
            /* de-paraxialize */
            (UD->l) = (UD->l)*(UD->n);
            (UD->m) = (UD->m)*(UD->n);
            }
         break;
      case 5:
      	/* ZEMAX wants a real ray trace to this surface */

			/* clear the multiple intercept test flag */
         miss_flag = 0;

         if (FD->cv == 0.0)
         	{
	         outofbounds:;
	         UD->ln =  0.0;
   	      UD->mn =  0.0;
      	   UD->nn = -1.0;
			   if (Refract(FD->n1, FD->n2, &UD->l, &UD->m, &UD->n, UD->ln, UD->mn, UD->nn)) return(-FD->surf);
            return(0);
            }
         /* okay, not a plane. */
         nx = (int) FD->param[1];
         ny = (int) FD->param[2];
         wx = FD->param[3];
         wy = FD->param[4];
         x = UD->x;
         y = UD->y;

         /* make sure nx and ny are both odd, otherwise the chief ray is a problem... */
         if (!nx&1)nx++;
         if (!ny&1)ny++;
         if (wx <= 0.0 || wy <= 0.0) return(-1);
         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &cx, &cy);
         if (error) goto outofbounds;

         /* Now, we illustrate an iterative method of finding
            the intercept for a general surface. */

			/* aspheric terms*/
         for (i = 1; i<=8; i++)
	        	{
				a[i] = FD->xdata[i]; // even terms
	  	      }

         try_again:;

			/* make sure we do at least 1 loop */
			t = 100.0;
         tp = 0.0;
         loop = 0;

         /* offset the coordinates */
         x -= cx;
         y -= cy;
         z = 0.0;

         while (fabs(t) > 1e-10)
         	{
				/*
            First, compute the sag using whatever the surface sag expression is.
            This is given the x and y starting points. The following block of code
            will change depending upon the surface shape, the rest of this iteration
            is typically common to all surface shapes.
            */

	         r2 = x * x + y * y;

	         alpha = 1.0 - (1.0 + FD->k)*FD->cv*FD->cv*r2;
	         if (alpha < 0.0) return(-1);

     			sag = (FD->cv*r2)/(1 + sqrt(alpha));

	         /* now the aspheric terms */
	         /* this is numerically inefficient, but easy to read */
	         if (a[1] != 0.0)
	         	{
					sag += a[1] * r2;
	            }
   	      if (a[2] != 0.0)
   	      	{
					sag += a[2] * r2 * r2;
	            }
	         if (a[3] != 0.0)
   	      	{
					sag += a[3] * r2 * r2 * r2;
	            }
	         if (a[4] != 0.0)
   	      	{
					sag += a[4] * r2 * r2 * r2 * r2;
	            }
	         if (a[5] != 0.0)
   	      	{
					sag += a[5] * r2 * r2 * r2 * r2 * r2;
	            }
   	      if (a[6] != 0.0)
   	      	{
					sag += a[6] * r2 * r2 * r2 * r2 * r2 * r2;
	            }
	         if (a[7] != 0.0)
	         	{
					sag += a[7] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
	            }
	         if (a[8] != 0.0)
	         	{
					sag += a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2;
	            }

            /* okay, now with sag in hand, how far are we away in z? */
            dz = sag - z;

            /* now compute how far along the z axis this is */
            /* note this will crash if n == 0!! */
            t = dz / (UD->n);

            /* propagate the additional "t" distance */
            x += UD->l*t;
            y += UD->m*t;
            z += UD->n*t;

            /* add in the optical path */
				tp += t;

            /* prevent infinte loop if no convergence */
            loop++;
            if (loop > 1000) return(-1);
            }

			UD->path = tp;

         /* okay, if the ray makes a steep angle of intercept,
         we may actually have hit the wrong element. Check it again! */
         if (miss_flag == 0)
         	{
            double new_cx, new_cy;

				/* avoid infinite loop */
            miss_flag = 1;

      		/* restore global coordinates prior to test */
	         UD->x = x + cx;
   	      UD->y = y + cy;
	         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &new_cx, &new_cy);
	         if (error) goto outofbounds;
            if (new_cx != cx || new_cy != cy)
            	{
               /* we hit the wrong one! */
               /* go back to the tangent plane */
               t = -tp;
					(x) = (UD->l) * t + (x);
					(y) = (UD->m) * t + (y);

		         x = x + cx;
	   	      y = y + cy;

               cx = new_cx;
               cy = new_cy;

               goto try_again;
               }
            }

         r2 = x * x + y * y;

			/* now do the normals */
         if (r2==0)
         	{
            UD->ln =  0;
            UD->mn =  0;
            UD->nn = -1;
            }
         else
         	{
				alpha = 1.0 - (1.0+FD->k)*FD->cv*FD->cv*r2;
	         if (alpha < 0) return(-1); /* ray misses */
	         alpha = sqrt(alpha);

	 	      mm = (FD->cv/(1.0+alpha))*(2.0 + (FD->cv*FD->cv*r2*(1.0+FD->k))/(alpha*(1.0+alpha)));

		      /* now the aspheric terms */
		      /* this is numerically inefficient, but easy to read */
		      if (a[1] != 0.0)
		         {
					mm += 2.0 * a[1];
		         }
		      if (a[2] != 0.0)
		         {
					mm += 4.0 * a[2] * r2;
		         }
		      if (a[3] != 0.0)
		         {
					mm += 6.0 * a[3] * r2 * r2;
	   	      }
	   	   if (a[4] != 0.0)
	      	   {
					mm += 8.0 * a[4] * r2 * r2 * r2;
	         	}
	      	if (a[5] != 0.0)
		         {
					mm += 10.0 * a[5] * r2 * r2 * r2 * r2;
		         }
		      if (a[6] != 0.0)
	   	      {
					mm += 12.0 * a[6] * r2 * r2 * r2 * r2 * r2;
  	      	   }
		      if (a[7] != 0.0)
		         {
					mm += 14.0 * a[7] * r2 * r2 * r2 * r2 * r2 * r2;
		         }
	   	   if (a[8] != 0.0)
	   	      {
					mm += 16.0 * a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
		         }

	         mx = x * mm;
	         my = y *	mm;

		      UD->nn = -sqrt(1/(1 + (mx*mx) + (my*my)));
		      UD->ln = -mx*UD->nn;
		      UD->mn = -my*UD->nn;
				}

         /* restore coordinates */
         UD->x = x + cx;
         UD->y = y + cy;
         UD->z = z;

         if (Refract(FD->n1, FD->n2, &UD->l, &UD->m, &UD->n, UD->ln, UD->mn, UD->nn)) return(-FD->surf);
         break;
      case 6:
      	/* ZEMAX wants the index, dn/dx, dn/dy, and dn/dz at the given x, y, z. */

         /* This is only required for gradient index surfaces, so return dummy values */
         UD->index = FD->n2;
         UD->dndx = 0.0;
         UD->dndy = 0.0;
         UD->dndz = 0.0;
      	break;
      case 7:
      	/* ZEMAX wants the "safe" data. */
         /* this is used by ZEMAX to set the initial values for all parameters and extra data */
         /* when the user first changes to this surface type. */
         /* this is the only time the DLL should modify the data in the FIXED_DATA FD structure */
         FD->param[1] = 5;
         FD->param[2] = 5;
         if (FD->cv == 0.0) FD->param[3] = 100.0;
         else FD->param[3] = 0.5/FD->cv;
         FD->param[4] = FD->param[3];
         for (i = 5; i <= 8; i++) FD->param[i] = 0.0;
         for (i = 1; i <= 200; i++) FD->xdata[i] = 0.0;
         break;
      }
   return 0;
   }

int Refract(double thisn, double nextn, double *l, double *m, double *n, double ln, double mn, double nn)
{
double nr, cosi, cosi2, rad, cosr, gamma;
if (thisn != nextn)
	{
	nr = thisn / nextn;
	cosi = fabs((*l) * ln + (*m) * mn + (*n) * nn);
	cosi2 = cosi * cosi;
	if (cosi2 > 1) cosi2 = 1;
	rad = 1 - ((1 - cosi2) * (nr * nr));
	if (rad < 0) return(-1);
	cosr = sqrt(rad);
	gamma = nr * cosi - cosr;
	(*l) = (nr * (*l)) + (gamma * ln);
	(*m) = (nr * (*m)) + (gamma * mn);
	(*n) = (nr * (*n)) + (gamma * nn);
	}
return 0;
}

int GetCellCenter(int nx, int ny, double wx, double wy, double x, double y, double *cx, double *cy)
{
double hwx, hwy, tx, ty;
int tnx, tny;

*cx = 0.0;
*cy = 0.0;

tnx = 0;
tny = 0;

hwx = 0.5*wx;
hwy = 0.5*wy;

/* do cx */
if (fabs(x) > hwx)
	{
   tx = x;
   tnx = 0;
   if (x > 0)
   	{
      while(tx > hwx)
      	{
         tnx++;
         tx -= wx;
         }
      }
   else
   	{
      while(tx < -hwx)
      	{
         tnx--;
         tx += wx;
         }
      }
   *cx = tnx*wx;
   }

/* do cy */
if (fabs(y) > hwy)
	{
   ty = y;
   tny = 0;
   if (y > 0)
   	{
      while(ty > hwy)
      	{
         tny++;
         ty -= wy;
         }
      }
   else
   	{
      while(ty < -hwy)
      	{
         tny--;
         ty += wy;
         }
      }
   *cy = tny*wy;
   }

/* are there this many cells? */
if (fabs(tnx) > (nx-1)/2 || fabs(tny) > (ny-1)/2) return -1;
else return 0;
}

