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

/*
Written by Kenneth E. Moore
May 8, 2003

Demonstrates use of FIXED_DATA2 to modify the polarization. See the real ray trace code (case 5) below.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface2(USER_DATA *UD, FIXED_DATA2 *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 */

int __declspec(dllexport) APIENTRY UserDefinedSurface2(USER_DATA *UD, FIXED_DATA2 *FD)
	{
   int i;
   double p2, alpha, power, a, b, c, rad, casp, t, zc;
   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,"Polarization");
               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, but display some dummy names to show how it works */
         /* returning a null string indicates that the parameter is unused. */
			UD->string[0] = '\0';
			if ((int) FD->param[0] == 0)
				{
				switch(FD->numb)
         		{
					case 0: strcpy(UD->string, "Method"); break;
					case 1: strcpy(UD->string, "s_rr"); break;
					case 2: strcpy(UD->string, "s_ri"); break;
					case 3: strcpy(UD->string, "s_tr"); break;
					case 4: strcpy(UD->string, "s_ti"); break;
					case 5: strcpy(UD->string, "p_rr"); break;
					case 6: strcpy(UD->string, "p_ri"); break;
					case 7: strcpy(UD->string, "p_tr"); break;
					case 8: strcpy(UD->string, "p_ti"); break;
					}
            }
			if ((int) FD->param[0] == 1)
				{
				switch(FD->numb)
         		{
					case 0: strcpy(UD->string, "Method"); break;
					case 1: strcpy(UD->string, "Exr"); break;
					case 2: strcpy(UD->string, "Exi"); break;
					case 3: strcpy(UD->string, "Eyr"); break;
					case 4: strcpy(UD->string, "Eyi"); break;
					case 5: strcpy(UD->string, "Ezr"); break;
					case 6: strcpy(UD->string, "Ezi"); 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);
            }
			else
				{
				/* okay, not a plane. */
				a = (UD->n) * (UD->n) * FD->k + 1;
				b = ((UD->n)/FD->cv) - (UD->x) * (UD->l) - (UD->y) * (UD->m);
				c = (UD->x) * (UD->x) + (UD->y) * (UD->y);
				rad = b * b - a * c;
				if (rad < 0) return(FD->surf);  /* ray missed this surface */
				if (FD->cv > 0) t = c / (b + sqrt(rad));
				else           t = c / (b - sqrt(rad));
				(UD->x) = (UD->l) * t + (UD->x);
				(UD->y) = (UD->m) * t + (UD->y);
				(UD->z) = (UD->n) * t + (UD->z);
				UD->path = t;
				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);
				}

			/* right here is where we optionally compute and return polarization data */
			if (FD->dbreserved[10] == -1.0)
				{
				int method;
				/*
				We are asked to compute the polarization data.
				There are two methods to do this:
				1) Compute the s and p coefficients for t and r (method 0 in this example).
				2) Directly modify the Ex, Ey, and Ez (method 1 in this example).

				Note if the coating data is provided by the DLL, the DLL must set FD->dbreserved[10] to 1.0.
				If the electric field data is modifed by the DLL, the DLL must set FD->dbreserved[10] to 2.0.
				*/

				method = (int) FD->param[0];
				
				if (method == 0)
					{
					/* compute s and p coefficients for t and r */
					double s_rr, s_ri, s_tr, s_ti;
					double p_rr, p_ri, p_tr, p_ti;
					/*
					Note there are eight values; 4 for s and 4 for p.
					There is r for reflection, t for transmission, and
					real and imagainary for each.
					*/

					/* note there is no default incoming data, use the values from the parameter spreadsheet */
					s_rr = FD->param[1]; // s pol reflect real
					s_ri = FD->param[2]; // s pol reflect imaginary
					s_tr = FD->param[3]; // s pol transmit real
					s_ti = FD->param[4]; // s pol transmit imaginary
					p_rr = FD->param[5]; // p pol reflect real
					p_ri = FD->param[6]; // p pol reflect imaginary
					p_tr = FD->param[7]; // p pol transmit real
					p_ti = FD->param[8]; // p pol transmit imaginary

					/* put the coating data in the array to pass back */
					FD->dbreserved[11] = s_rr;
					FD->dbreserved[12] = s_ri;
					FD->dbreserved[13] = s_tr;
					FD->dbreserved[14] = s_ti;
					FD->dbreserved[15] = p_rr;
					FD->dbreserved[16] = p_ri;
					FD->dbreserved[17] = p_tr;
					FD->dbreserved[18] = p_ti;
					/* now tell ZEMAX we are providing the coating data! */
					/* the "1" means we are specifying the coating data */
					FD->dbreserved[10] = 1.0;
					}

				/* use this code for method 1 */
				if (method == 1)
					{
					/* Specify E directly.. this is dangerous since energy might not be conserved... but this is how it's done... */
					double exr, exi, eyr, eyi, ezr, ezi;
					/* this is the incoming data, shown for illustration only */
					exr = FD->dbreserved[11];
					exi = FD->dbreserved[12];
					eyr = FD->dbreserved[13];
					eyi = FD->dbreserved[14];
					ezr = FD->dbreserved[15];
					ezi = FD->dbreserved[16];

					/* replace with the data in the spreadsheet */
					exr = FD->param[1];
					exi = FD->param[2];
					eyr = FD->param[3];
					eyi = FD->param[4];
					ezr = FD->param[5];
					ezi = FD->param[6];

					/* put the specified data back */
					FD->dbreserved[11] = exr;
					FD->dbreserved[12] = exi;
					FD->dbreserved[13] = eyr;
					FD->dbreserved[14] = eyi;
					FD->dbreserved[15] = ezr;
					FD->dbreserved[16] = ezi;
					/* now tell ZEMAX we modified the E field directly! */
					/* the "2" means we are specifying the E field data */
					FD->dbreserved[10] = 2.0;
					}

				}

         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 = 0; i <= 12; i++) FD->param[i] = 0.0;
         for (i = 1; i <= 200; i++) FD->xdata[i] = 0.0;
			FD->param[3] = 1.0;
			FD->param[7] = 1.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;
}

