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

/*
Written by Kenneth E. Moore
July 22, 1998

Modified Oct 27, 1998 to test for huge ray coordinates; returns error if rays get too big; see BIG_RAYS below.

Modified April 11, 2002 to add "slow factor". If 1, ray traces at full speed, but may not converge of steep parts/rays. Set to 0.5 or less to improve convergence

*/

#define BIG_RAY 1.0E+010

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 an ogive surface shape, which is just like a standard
surface except the axis of rotation is offset by a distance ro.

This DLL can be used as a template for any surface form that requires iteration
to determine the surface intercept.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, loop;
   double p2, alpha1, alpha2, power, r, r0, rpr0, t, tp, x, y, z, dz, sag, slowfactor;

   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,"Ogive");
               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 except parameter 1 */
         /* returning a null string indicates that the parameter is unused. */
         switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "r 0");
            	break;
            case 2:
            	strcpy(UD->string, "slowfactor");
            	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)
         	{
            default:
            	UD->string[0] = '\0';
            	break;
            }
      	break;
      case 3:
      	/* ZEMAX wants to know the sag of the surface */
         /* 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);

         r0 = FD->param[1];
         p2 = UD->x * UD->x + UD->y * UD->y;
         r = sqrt(p2);
         rpr0 = r + r0;

         alpha1 = 1 - (1+FD->k)*FD->cv*FD->cv*rpr0*rpr0;
         if (alpha1 < 0) return(-1);

         alpha2 = 1 - (1+FD->k)*FD->cv*FD->cv*r0*r0;
         if (alpha2 < 0) return(-1);

         UD->sag1 = FD->cv * ( rpr0*rpr0/(1 + sqrt(alpha1)) - r0*r0/(1 + sqrt(alpha2)) );
         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 the ogive surface, we will assume ro = 0 for this computation */

         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;

         r0 = FD->param[1];

			/* slow the iteration down if the ray is steep */
			slowfactor = FD->param[2];
			if (slowfactor <= 0 || slowfactor > 1.0) slowfactor = 1.0;

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

				/* make sure we don't iterate out of control */
            if (p2 >= BIG_RAY) return(-1);

	         r = sqrt(p2);
	         rpr0 = r + r0;

	         alpha1 = 1 - (1+FD->k)*FD->cv*FD->cv*rpr0*rpr0;
	         if (alpha1 < 0)
            	{
               /* we may be too high up on the surface. Try a value near the edge */
               alpha1 = 0.0001;
               }

	         alpha2 = 1 - (1+FD->k)*FD->cv*FD->cv*r0*r0;
	         if (alpha2 < 0)
            	{
               /* we may be too high up on the surface. Try a value near the edge */
               alpha2 = 0.0001;
               }

   	      sag = FD->cv * ( rpr0*rpr0/(1 + sqrt(alpha1)) - r0*r0/(1 + sqrt(alpha2)) );

            /* 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);
            t = dz * slowfactor;

            /* 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 at 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, this is the tricky part  */

         p2 = x * x + y * y;
         r = sqrt(p2);
         rpr0 = r + r0;

			if (r == 0.0)
         	{
				UD->ln = 0.0;
				UD->mn = 0.0;
				UD->nn = -1.0;
            }
         else
         	{
            double mx, my, temp;
	         alpha1 = 1 - (1+FD->k)*FD->cv*FD->cv*rpr0*rpr0;
	         if (alpha1 < 0) return -1;
            alpha1 = sqrt(alpha1);
            temp  = (FD->cv/r)*rpr0/(1 + alpha1);
            temp *= (2.0 + FD->cv*FD->cv*rpr0*rpr0*(1+FD->k)/(alpha1*(1.0+alpha1)));
	         mx = x*temp;
	         my = y*temp;
				UD->nn = -sqrt(1.0/(1.0 + (mx*mx) + (my*my)));
				UD->ln = -mx*UD->nn;
				UD->mn = -my*UD->nn;
            }

         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;
}

