us_anamr.c

zemax源码： This DLL models an anamorphic aspheric surface. This surface is essentially an even asp

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

/*
May 28, 1998 Written by Kenneth E. Moore
Sep 15, 2000 decreased iteration step size, and added sphere intercept. JFS
May 15, 2001 corrected surface derivative. KEM
*/

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 anamorphic aspheric surface.

This surface is essentially an even aspheric surface with different terms for
the x and y directions.

The sag is given by:

Z = ((CX*x*x)+(CY*y*y)) / (1 + sqrt(1-((1+KX)*CX*CX*x*x)-((1+KY)*CY*CY*y*y)))

+ AR*( (1 - AP)*x*x + (1 + AP)*y*y )^2
+ BR*( (1 - BP)*x*x + (1 + BP)*y*y )^3
+ CR*( (1 - CP)*x*x + (1 + CP)*y*y )^4
+ DR*( (1 - DP)*x*x + (1 + DP)*y*y )^5

Note the terms AR, BR, CR, and DR ... have units of length to the -3, -5, -7, and -9 power.

The terms AP, BP, CP, and DP are dimensionless.

The surface is rotationally symmetric only if AP = BP = CP = DP == 0 and CX = CY and KX = KY.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, loop;
   double alpha, xpower, ypower, t, tp, x, y, z, dz, sag, mx, my;
   double AR, BR, CR, DR, AP, BP, CP, DP, CX, CY, KX, KY, X2, Y2, temp;


   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,"Anamorphic");
               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. */
         /* Only "q" in parameter 1 is used for this surface type */
         /* returning a null string indicates that the parameter is unused. */
         switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "Cx");
               break;
            case 2:
            	strcpy(UD->string, "Cy");
               break;
            case 3:
            	strcpy(UD->string, "Kx");
               break;
            case 4:
            	strcpy(UD->string, "Ky");
               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. */
         /* returning a null string indicates that the extradata value is unused. */
         switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "AR");
               break;
            case 2:
            	strcpy(UD->string, "BR");
               break;
            case 3:
            	strcpy(UD->string, "CR");
               break;
            case 4:
            	strcpy(UD->string, "DR");
               break;
            case 5:
            	strcpy(UD->string, "AP");
               break;
            case 6:
            	strcpy(UD->string, "BP");
               break;
            case 7:
            	strcpy(UD->string, "CP");
               break;
            case 8:
            	strcpy(UD->string, "DP");
               break;
            default:
            	UD->string[0] = '\0';
            	break;
            }
      	break;
      case 3:
      	/* ZEMAX wants to know the sag of the surface */

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

         CX = FD->param[1];
         CY = FD->param[2];
         KX = FD->param[3];
         KY = FD->param[4];
         AR = FD->xdata[1];
         BR = FD->xdata[2];
         CR = FD->xdata[3];
         DR = FD->xdata[4];
         AP = FD->xdata[5];
         BP = FD->xdata[6];
         CP = FD->xdata[7];
         DP = FD->xdata[8];

         x = UD->x;
         y = UD->y;

         X2 = x*x;
         Y2 = y*y;

         alpha = 1.0 - (1.0 + KX)*CX*CX*X2 - (1.0 + KY)*CY*CY*Y2;
         if (alpha < 0) return(-1);

			UD->sag1 = (CX*X2 + CY*Y2)/(1.0 + sqrt(alpha));

         /* now the aspheric terms */
         if (AR != 0.0)
         	{
	         temp = (1.0 - AP)*X2 + (1.0 + AP)*Y2;
				UD->sag1 += AR*temp*temp;
            }

         if (BR != 0.0)
         	{
	         temp = (1.0 - BP)*X2 + (1.0 + BP)*Y2;
				UD->sag1 += BR*temp*temp*temp;
            }

         if (CR != 0.0)
         	{
	         temp = (1.0 - CP)*X2 + (1.0 + CP)*Y2;
	         temp *= temp;
	         temp *= temp;
				UD->sag1 += CR*temp;
            }

         if (DR != 0.0)
         	{
	         temp = (1.0 - DP)*X2 + (1.0 + DP)*Y2;
				UD->sag1 += DR*temp*temp*temp*temp*temp;
            }

         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. */
         /* we will ignore the aspheric terms. */

         UD->ln =  0.0;
         UD->mn =  0.0;
         UD->nn = -1.0;
         CX = FD->param[1];
         CY = FD->param[2];
         xpower = (FD->n2 - FD->n1)*CX;
         ypower = (FD->n2 - FD->n1)*CY;
         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)*xpower)/(FD->n2);
            (UD->m) = (FD->n1*(UD->m) - (UD->y)*ypower)/(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 */
         /* do not allow n == 0 */
         if (UD->n == 0.0) return -1;

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

         CX = FD->param[1];
         CY = FD->param[2];
         KX = FD->param[3];
         KY = FD->param[4];
         AR = FD->xdata[1];
         BR = FD->xdata[2];
         CR = FD->xdata[3];
         DR = FD->xdata[4];
         AP = FD->xdata[5];
         BP = FD->xdata[6];
         CP = FD->xdata[7];
         DP = FD->xdata[8];

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

			t = 100.0;
         tp = 0.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.
            */

	         X2 = x*x;
   	      Y2 = y*y;

	         alpha = 1.0 - (1.0 + KX)*CX*CX*X2 - (1.0 + KY)*CY*CY*Y2;

				if (alpha < 0 && loop == 0 && (CX != 0 || CY != 0)) 
					{
					/* try tracing from tangent plane to sphere */
					/* to find new starting x,y if first loop fales.*/
					double tempcv, tempk, a, b, c, rad;

					/* Start with curvature in y */						
					if (CY != 0)
						{
						tempcv = CY;
						tempk = KY;
						}
					else
						{
						tempcv = CX;
						tempk = KX;
						}
					a = (UD->n) * (UD->n) * tempk + 1;
					b = ((UD->n)/tempcv) - (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 && CY != CX && CX != 0)
						{
						/* Try curvature in x */
						tempcv = CX;
						tempk = KX;
						a = (UD->n) * (UD->n) * tempk + 1;
						b = ((UD->n)/tempcv) - (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 (tempcv > 0) t = c / (b + sqrt(rad));
					else            t = c / (b - sqrt(rad));
					x = (UD->l) * t + (UD->x);
					y = (UD->m) * t + (UD->y);
					z = (UD->n) * t + (UD->z);
					tp += t;

					/* test again */
					X2 = x*x;
					Y2 = y*y;
					alpha = 1.0 - (1.0 + KX)*CX*CX*X2 - (1.0 + KY)*CY*CY*Y2;
					}
				
	         if (alpha < 0) return(FD->surf);

				sag = (CX*X2 + CY*Y2)/(1.0 + sqrt(alpha));

	         /* now the aspheric terms */
	         if (AR != 0.0)
   	      	{
	   	      temp = (1.0 - AP)*X2 + (1.0 + AP)*Y2;
					sag += AR*temp*temp;
               }

	         if (BR != 0.0)
   	      	{
	   	      temp = (1.0 - BP)*X2 + (1.0 + BP)*Y2;
					sag += BR*temp*temp*temp;
               }

	         if (CR != 0.0)
   	      	{
		         temp = (1.0 - CP)*X2 + (1.0 + CP)*Y2;
		         temp *= temp;
	   	      temp *= temp;
					sag += CR*temp;
               }

	         if (DR != 0.0)
   	      	{
		         temp = (1.0 - DP)*X2 + (1.0 + DP)*Y2;
					sag += DR*temp*temp*temp*temp*temp;
               }

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

            /* 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(FD->surf-1);
            }

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

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


			/* now do the normals */
         alpha = 1.0 - (1.0 + KX)*CX*CX*X2 - (1.0 + KY)*CY*CY*Y2;
         if (alpha < 0) return(FD->surf);
         alpha = sqrt(alpha);

			/* old incorrect expression */
			/*
			temp = CX*CX*X2+CY*CY*Y2;
         mx = (CX*x/(1.0+alpha))*(2.0 + temp*(1.0+KX)/(alpha*(1.0+alpha)));
         my = (CY*y/(1.0+alpha))*(2.0 + temp*(1.0+KY)/(alpha*(1.0+alpha)));
			*/

			/* corrected expression 5-15-01 */
			temp = CX*X2+CY*Y2;
         mx = (CX*x/(1.0+alpha))*(2.0 + CX*temp*(1.0+KX)/(alpha*(1.0+alpha)));
         my = (CY*y/(1.0+alpha))*(2.0 + CY*temp*(1.0+KY)/(alpha*(1.0+alpha)));

         if (AR != 0.0)
         	{
	         temp = (1.0 - AP)*X2 + (1.0 + AP)*Y2;
	         temp = 2.0*AR*temp;
	         mx += temp*2.0*(1.0 - AP)*x;
	         my += temp*2.0*(1.0 + AP)*y;
            }

         if (BR != 0.0)
         	{
	         temp = (1.0 - BP)*X2 + (1.0 + BP)*Y2;
   	      temp = 3.0*BR*temp*temp;
      	   mx += temp*2.0*(1.0 - BP)*x;
         	my += temp*2.0*(1.0 + BP)*y;
            }

         if (CR != 0.0)
         	{
	         temp = (1.0 - CP)*X2 + (1.0 + CP)*Y2;
   	      temp = 4.0*CR*temp*temp*temp;
      	   mx += temp*2.0*(1.0 - CP)*x;
         	my += temp*2.0*(1.0 + CP)*y;
            }

         if (DR != 0.0)
         	{
	         temp = (1.0 - DP)*X2 + (1.0 + DP)*Y2;
   	      temp = 5.0*DR*temp*temp*temp*temp;
      	   mx += temp*2.0*(1.0 - DP)*x;
         	my += temp*2.0*(1.0 + DP)*y;
            }

         UD->nn = -sqrt(1/(1 + (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;
}