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

/*
Written by Kenneth E. Moore
August 10, 2000
*/

#define DEGTORAD 0.0174532925199433

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 plane surface that refracts light as though it were a DMD */

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, k, p, nx, ny, ixp, iyp, ip, is;
	double A, B, wx, wy, dx, dy, tnx, tny, tnz, tnr;

   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,"DMD");
               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. */
         switch(FD->numb)
         	{
            default:
            	UD->string[0] = '\0';
            	break;
            case 1:
            	strcpy(UD->string,"# X pixels");
               break;
            case 2:
            	strcpy(UD->string,"# Y pixels");
               break;
            case 3:
            	strcpy(UD->string,"X-Width");
               break;
            case 4:
            	strcpy(UD->string,"Y-Width");
               break;
            case 5:
            	strcpy(UD->string,"Angle 0");
               break;
            case 6:
            	strcpy(UD->string,"Angle 1");
               break;
            case 7:
            	strcpy(UD->string,"Angle 2");
               break;
            case 8:
            	strcpy(UD->string,"Rot Angle");
               break;
            }
      	break;
      case 2:
      	/* ZEMAX is requesting the names of the extra data columns */
         /* the value FD->numb will indicate which value ZEMAX wants. */
         switch(FD->numb)
         	{
            default:
					nx = (int)floor(0.5+FD->param[1]);
					ny = (int)floor(0.5+FD->param[2]);
					ip = nx*ny;
					if (FD->xdata[1] == 0)
						{
						/* by row */
						ip = ny;
						}
					if (FD->xdata[1] == 1)
						{
						/* by column */
						ip = nx;
						}
					k = 2 + (ip / 15);
					if (FD->numb <= k)
						{
						sprintf(UD->string, "Pixels %i - %i", 1+15*(FD->numb-2), 15*(FD->numb-1));
						if (FD->xdata[1] == 0)
							{
							/* by row */
							sprintf(UD->string, "Rows %i - %i", 1+15*(FD->numb-2), 15*(FD->numb-1));
							}
						if (FD->xdata[1] == 1)
							{
							/* by column */
							sprintf(UD->string, "Columns %i - %i", 1+15*(FD->numb-2), 15*(FD->numb-1));
							}
						}
					else
						{
	            	UD->string[0] = '\0';
						}			

            	break;
            case 1:
            	strcpy(UD->string,"P flag");
               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 */
			return(0);
      	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;
         if ((UD->n) != 0.0)
         	{
            (UD->l) = (UD->l)/(UD->n);
            (UD->m) = (UD->m)/(UD->n);

            (UD->l) = (FD->n1*(UD->l))/(FD->n2);
            (UD->m) = (FD->n1*(UD->m))/(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 */
	      UD->ln =  0.0;
   	   UD->mn =  0.0;
      	UD->nn = -1.0;

			/* what pixel did we hit ? */
			nx = (int)floor(0.5+FD->param[1]);
			ny = (int)floor(0.5+FD->param[2]);
			if (nx < 1 || ny < 1)
         	{
			   if (Refract(FD->n1, FD->n2, &UD->l, &UD->m, &UD->n, UD->ln, UD->mn, UD->nn)) return(-FD->surf);
            return(0);
            }
			wx = FD->param[3];
			wy = FD->param[4];
			dx = wx / nx;
			dy = wy / ny;

			ixp = 1 + (int) floor( (UD->x + 0.5*wx) / dx );
			iyp = 1 + (int) floor( (UD->y + 0.5*wy) / dy );
			if (ixp < 1 || ixp > nx || iyp < 1 || iyp > ny)
         	{
			   if (Refract(FD->n1, FD->n2, &UD->l, &UD->m, &UD->n, UD->ln, UD->mn, UD->nn)) return(-FD->surf);
            return(0);
            }

			ip = ixp + (iyp-1)*nx;
			p = 0;
			/* ok, we hit a pixel */
			/* how are we addressing the pixels? */
			if (FD->xdata[1] == 0)
				{
				/* by row */
				ip = iyp;
				}
			if (FD->xdata[1] == 1)
				{
				/* by column */
				ip = ixp;
				}
			if (FD->xdata[1] == 2)
				{
				/* by individual pixel */
				ip = ixp + (iyp-1)*nx;
				}
			/* where is the data for ip stored? 15 pixels per data value... */
			k = 2;
			while(ip > 15)
				{
				ip -= 15;
				k++;
				}

			if (k < 2 || k > 200)
				{
			   if (Refract(FD->n1, FD->n2, &UD->l, &UD->m, &UD->n, UD->ln, UD->mn, UD->nn)) return(-FD->surf);
            return(0);
				}

			is = (int) FD->xdata[k];

			k = 1;
			for (i = 1; i <= ip; i++) k *= 3;
			while (is >= k) is -= k;
			k /= 3;
			p = is / k;

			/* p is now either 0, 1, or 2 */
			B = DEGTORAD * FD->param[8];
			A = DEGTORAD * FD->param[5];
			if (p == 1) A = DEGTORAD * FD->param[6];
			if (p == 2) A = DEGTORAD * FD->param[7];

      	tnz = -cos(A);
			tnr =  sin(A);
			tnx = tnr*sin(B);
			tny = tnr*cos(B);
			UD->ln = tnx;
   		UD->mn = tny;
      	UD->nn = tnz;

         /* now refract */
			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] = 10;
			FD->param[2] = 10;
			FD->param[3] = 10;
			FD->param[4] = 10;
			FD->param[5] = -10;
			FD->param[6] = 0;
			FD->param[7] = +10;
			FD->param[8] = 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;
}

