evenarray.c

来自「光学设计分析软件zemax的非球面旋转对称透镜的源代码」· C语言 代码 · 共 581 行 · 第 1/2 页

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

/*
Written by James Sutter
Sept 21, 2001

This DLL models an arbitrary number of lens elements in a lens array.
The individual lenses are even aspheres.

The user provides the number of elements in x and y, the size in x and y of
each element, and the radius, conic, asphere coefficients, and glass.

The asphere sag in each element is given by:

Z = (c*r*r) / (1+(1-((1+k)*c*c*r*r))^ 1/2 ) + a2*(r)^2 + a4*(r)^4 + ... a8*(r)^16

Note the terms a2 ... a16 have units of length to the -1, -3 ... -15 power.

This surface breaks up the beam into numerous separate
beams, and so most ZEMAX features will fail to work with this surface.

However, the spot diagrams, image analysis, etc, all work okay.

*/

int __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD);
int GetCellCenter(int nx, int ny, double wx, double wy, double x, double y, double *cx, double *cy);

/* 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 lens array surface type */

int  __declspec(dllexport) APIENTRY UserDefinedSurface(USER_DATA *UD, FIXED_DATA *FD)
	{
   int i, nx, ny, error, miss_flag;
   double alpha, power, t;
   double wx, wy, cx, cy, x, y, z;
   int loop;
   double r2, tp, dz, sag, mm, mx, my;
   double a[9];
   
	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,"Even Array");
               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. */
         /* they are all "Unused" for this surface type */
         /* returning a null string indicates that the parameter is unused. */
         switch(FD->numb)
         	{
            case 1:
            	strcpy(UD->string, "Number X");
               break;
            case 2:
            	strcpy(UD->string, "Number Y");
               break;
            case 3:
            	strcpy(UD->string, "Width X");
               break;
            case 4:
            	strcpy(UD->string, "Height Y");
               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)
         	{
            case 1:
            	strcpy(UD->string, "2nd Order Term");
               break;
            case 2:
            	strcpy(UD->string, "4th Order Term");
               break;
            case 3:
            	strcpy(UD->string, "6th Order Term");
               break;
            case 4:
            	strcpy(UD->string, "8th Order Term");
               break;
            case 5:
            	strcpy(UD->string, "10th Order Term");
               break;
            case 6:
            	strcpy(UD->string, "12th Order Term");
               break;
            case 7:
            	strcpy(UD->string, "14th Order Term");
               break;
            case 8:
            	strcpy(UD->string, "16th Order Term");
               break;
            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. */

			/* aspheric terms*/
         for (i = 1; i<=8; i++)
	        	{
				a[i] = FD->xdata[i]; // even terms
	  	      }

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

			/* if a plane, just return */
			if (FD->cv == 0) return(0);

         /* figure out the center coordinates of which "cell" we are in */
         nx = (int) FD->param[1];
         ny = (int) FD->param[2];
         wx = FD->param[3];
         wy = FD->param[4];
         x = UD->x;
         y = UD->y;

         /* make sure nx and ny are both odd, otherwise the chief ray is a problem... */
         if (!nx&1)nx++;
         if (!ny&1)ny++;
         if (wx <= 0.0 || wy <= 0.0) return(-1);
         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &cx, &cy);

         if (error) return(0);

         /* offset the coordinates */
         x -= cx;
         y -= cy;

         r2 = x*x + y*y;
         alpha = 1 - (1+FD->k)*FD->cv*FD->cv*r2;
         if (alpha < 0)
         	{
            /* lens is not big enough to fill the cell, radius is too short */
            /* assume it is plane between lenses */
            alpha = 0.0;
            }
         UD->sag1 = (FD->cv*r2)/(1 + sqrt(alpha));

         /* now the aspheric terms */
         /* this is numerically inefficient, but easy to read */
         if (a[1] != 0.0)
         	{
				UD->sag1 += a[1] * r2;
            }
         if (a[2] != 0.0)
         	{
				UD->sag1 += a[2] * r2 * r2;
            }
         if (a[3] != 0.0)
         	{
				UD->sag1 += a[3] * r2 * r2 * r2;
            }
         if (a[4] != 0.0)
         	{
				UD->sag1 += a[4] * r2 * r2 * r2 * r2;
            }
         if (a[5] != 0.0)
         	{
				UD->sag1 += a[5] * r2 * r2 * r2 * r2 * r2;
            }
         if (a[6] != 0.0)
         	{
				UD->sag1 += a[6] * r2 * r2 * r2 * r2 * r2 * r2;
            }
         if (a[7] != 0.0)
         	{
				UD->sag1 += a[7] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
            }
         if (a[8] != 0.0)
         	{
				UD->sag1 += a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2;
            }
			/* forget supporting a hyper hemisphere! */
         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 a lens array, only consider the single lens on axis */
         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 */

			/* clear the multiple intercept test flag */
         miss_flag = 0;

         if (FD->cv == 0.0)
         	{
	         outofbounds:;
	         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. */
         nx = (int) FD->param[1];
         ny = (int) FD->param[2];
         wx = FD->param[3];
         wy = FD->param[4];
         x = UD->x;
         y = UD->y;

         /* make sure nx and ny are both odd, otherwise the chief ray is a problem... */
         if (!nx&1)nx++;
         if (!ny&1)ny++;
         if (wx <= 0.0 || wy <= 0.0) return(-1);
         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &cx, &cy);
         if (error) goto outofbounds;

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

			/* aspheric terms*/
         for (i = 1; i<=8; i++)
	        	{
				a[i] = FD->xdata[i]; // even terms
	  	      }

         try_again:;

			/* make sure we do at least 1 loop */
			t = 100.0;
         tp = 0.0;
         loop = 0;

         /* offset the coordinates */