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

/*
Written by Kenneth E. Moore
April 22, 1997
*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD);

/* 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 standard ZEMAX surface type, either plane, sphere, or conic.
It is just like US_STAND, except it uses "dumb" iteration to find the intercept
point rather than a closed-form expression. This can be used as a template for
defining aspheric UDS DLL's which require iteration.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, loop;
   double p2, alpha, power, rad, casp, t, zc, tp, x, y, z, dz, sag;

   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,"Iter Std Dll");
               break;
            case 1:
            	/* ZEMAX wants to know if this surface is rotationally symmetric */
               /* it is, so return any character in the string; otherwise, return a null string */
            	strcpy(UD->string, "1");
               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)
         	{
            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)
         	{
            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. */

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

			/* if a plane, just return */
			if (FD->cv == 0) return(0);
         p2 = UD->x * UD->x + UD->y * UD->y;
         alpha = 1 - (1+FD->k)*FD->cv*FD->cv*p2;
         if (alpha < 0) return(-1);
         UD->sag1 = (FD->cv*p2)/(1 + sqrt(alpha));
         if (alpha != 1.0) UD->sag2 = (FD->cv*p2)/(1 - sqrt(alpha));
      	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 */

         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 */
         if (FD->cv == 0.0)
         	{
	         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. */
         /* do not allow n == 0 */
         if (UD->n == 0.0) return -1;



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

			/* make sure we do at least 1 loop */

   		tp = 0.0;
			t = 100.0;
         x = UD->x;
         y = UD->y;
         z = UD->z;
         loop = 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.

            If you have (and understand) a faster routine, such as a Newton-Raphson technique,
            or if you know a way of computing a close first guess, then use that instead!

            This basic method is not optimal speed wise, but it works fine for most simple surface
            shapes that have clearly defined intercept points. This method may fail for surfaces with
            multiple possible solutions, like toroids, or highly aspheric shapes.
            */

	         p2 = x * x + y * y;
	         alpha = 1.0 - (1.0 + FD->k)*FD->cv*FD->cv*p2;
	         if (alpha < 0.0) return(-1);
	         sag = (FD->cv*p2)/(1.0 + sqrt(alpha));

            /* 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);

            /* for some aspheres, it is safer to use dz directly, as it is a smaller number */
            /* the convergence will be slower if the ray angle is steep, fast if near parallel to the axis */


            /* 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);
            }

         /* okay, we should be a the intercept coordinates now */
         UD->x = x;
         UD->y = y;
         UD->z = z;

      	/* don't forget the path! */
			UD->path = tp;

			/* now do the normals in the usual way */
			zc = (UD->z) * FD->cv;
			rad = zc * FD->k * (zc * (FD->k + 1) - 2) + 1;
			casp = FD->cv / sqrt(rad);
			UD->ln = (UD->x) * casp;
			UD->mn = (UD->y) * casp;
			UD->nn = ((UD->z) - ((1/FD->cv) - (UD->z) * FD->k)) * casp;
         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 */
         for (i = 1; 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;
}

