evenarray.c

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

         x -= cx;
         y -= cy;
         z = 0.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.
            */

	         r2 = x * x + y * y;

	         alpha = 1.0 - (1.0 + FD->k)*FD->cv*FD->cv*r2;
	         if (alpha < 0.0) return(-1);

     			sag = (FD->cv*r2)/(1 + sqrt(alpha));

	         /* now the aspheric terms */
	         /* this is numerically inefficient, but easy to read */
	         if (a[1] != 0.0)
	         	{
					sag += a[1] * r2;
	            }
   	      if (a[2] != 0.0)
   	      	{
					sag += a[2] * r2 * r2;
	            }
	         if (a[3] != 0.0)
   	      	{
					sag += a[3] * r2 * r2 * r2;
	            }
	         if (a[4] != 0.0)
   	      	{
					sag += a[4] * r2 * r2 * r2 * r2;
	            }
	         if (a[5] != 0.0)
   	      	{
					sag += a[5] * r2 * r2 * r2 * r2 * r2;
	            }
   	      if (a[6] != 0.0)
   	      	{
					sag += a[6] * r2 * r2 * r2 * r2 * r2 * r2;
	            }
	         if (a[7] != 0.0)
	         	{
					sag += a[7] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
	            }
	         if (a[8] != 0.0)
	         	{
					sag += a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2;
	            }

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

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

			UD->path = tp;

         /* okay, if the ray makes a steep angle of intercept,
         we may actually have hit the wrong element. Check it again! */
         if (miss_flag == 0)
         	{
            double new_cx, new_cy;

				/* avoid infinite loop */
            miss_flag = 1;

      		/* restore global coordinates prior to test */
	         UD->x = x + cx;
   	      UD->y = y + cy;
	         error = GetCellCenter(nx, ny, wx, wy, UD->x, UD->y, &new_cx, &new_cy);
	         if (error) goto outofbounds;
            if (new_cx != cx || new_cy != cy)
            	{
               /* we hit the wrong one! */
               /* go back to the tangent plane */
               t = -tp;
					(x) = (UD->l) * t + (x);
					(y) = (UD->m) * t + (y);

		         x = x + cx;
	   	      y = y + cy;

               cx = new_cx;
               cy = new_cy;

               goto try_again;
               }
            }

         r2 = x * x + y * y;

			/* now do the normals */
         if (r2==0)
         	{
            UD->ln =  0;
            UD->mn =  0;
            UD->nn = -1;
            }
         else
         	{
				alpha = 1.0 - (1.0+FD->k)*FD->cv*FD->cv*r2;
	         if (alpha < 0) return(-1); /* ray misses */
	         alpha = sqrt(alpha);

	 	      mm = (FD->cv/(1.0+alpha))*(2.0 + (FD->cv*FD->cv*r2*(1.0+FD->k))/(alpha*(1.0+alpha)));

		      /* now the aspheric terms */
		      /* this is numerically inefficient, but easy to read */
		      if (a[1] != 0.0)
		         {
					mm += 2.0 * a[1];
		         }
		      if (a[2] != 0.0)
		         {
					mm += 4.0 * a[2] * r2;
		         }
		      if (a[3] != 0.0)
		         {
					mm += 6.0 * a[3] * r2 * r2;
	   	      }
	   	   if (a[4] != 0.0)
	      	   {
					mm += 8.0 * a[4] * r2 * r2 * r2;
	         	}
	      	if (a[5] != 0.0)
		         {
					mm += 10.0 * a[5] * r2 * r2 * r2 * r2;
		         }
		      if (a[6] != 0.0)
	   	      {
					mm += 12.0 * a[6] * r2 * r2 * r2 * r2 * r2;
  	      	   }
		      if (a[7] != 0.0)
		         {
					mm += 14.0 * a[7] * r2 * r2 * r2 * r2 * r2 * r2;
		         }
	   	   if (a[8] != 0.0)
	   	      {
					mm += 16.0 * a[8] * r2 * r2 * r2 * r2 * r2 * r2 * r2;
		         }

	         mx = x * mm;
	         my = y *	mm;

		      UD->nn = -sqrt(1/(1 + (mx*mx) + (my*my)));
		      UD->ln = -mx*UD->nn;
		      UD->mn = -my*UD->nn;
				}

         /* restore coordinates */
         UD->x = x + cx;
         UD->y = y + cy;
         UD->z = z;

         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] = 5;
         FD->param[2] = 5;
         if (FD->cv == 0.0) FD->param[3] = 100.0;
         else FD->param[3] = 0.5/FD->cv;
         FD->param[4] = FD->param[3];
         for (i = 5; 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;
}

int GetCellCenter(int nx, int ny, double wx, double wy, double x, double y, double *cx, double *cy)
{
double hwx, hwy, tx, ty;
int tnx, tny;

*cx = 0.0;
*cy = 0.0;

tnx = 0;
tny = 0;

hwx = 0.5*wx;
hwy = 0.5*wy;

/* do cx */
if (fabs(x) > hwx)
	{
   tx = x;
   tnx = 0;
   if (x > 0)
   	{
      while(tx > hwx)
      	{
         tnx++;
         tx -= wx;
         }
      }
   else
   	{
      while(tx < -hwx)
      	{
         tnx--;
         tx += wx;
         }
      }
   *cx = tnx*wx;
   }

/* do cy */
if (fabs(y) > hwy)
	{
   ty = y;
   tny = 0;
   if (y > 0)
   	{
      while(ty > hwy)
      	{
         tny++;
         ty -= wy;
         }
      }
   else
   	{
      while(ty < -hwy)
      	{
         tny--;
         ty += wy;
         }
      }
   *cy = tny*wy;
   }

/* are there this many cells? */
if (fabs(tnx) > (nx-1)/2 || fabs(tny) > (ny-1)/2) return -1;
else return 0;
}