source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/Ntrace_lenses.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 36.7 KB
Line 
1/* -------------------------------------------------------------------------
2 *   trace_lenses.c
3 *
4 *   --- main program of raytracing in the Fresnel lens system
5 *
6 * Copyright (c) 2000-2012 N.Sakaki, Y.Takizawa, Y.Kawasaki
7 * All rights reserved
8 * $Id$
9 * -------------------------------------------------------------------------
10 */
11#define NTRACE
12#include <stdio.h>
13#include <stdlib.h>
14#include <string.h>
15#include <math.h>
16#include "EConst.hh"
17#include "EsafRandom.hh"
18#include "Nrand_gen.hh"
19#include "Nspline.hh"
20#include "Nutils.hh"
21#include "Ntrace_lenses.hh"
22#include "Nphoton.hh"
23//#include "EConst.hh"
24
25#define RT_SEARCH_ITER_MAX        100
26#define RT_INTERACTION_MAX         20
27//#define DEBUG
28using namespace NTraceLens;
29
30static double ref(double n1, double n2, double ci, double ct);
31static int check_intersect_cylinder(double r1, double z1, double r2, double z2,
32                                    NPhoton *photon, double *dis);
33static int check_intersect_cone    (double r1, double z1, double r2, double z2,
34                                    NPhoton *photon, double *dis);
35static int CheckIris(tel_param *para, NPhoton *ph, int flag);
36
37static int fGen;
38static double fEffectiveEfficiencyFactor;
39
40void  NTraceLens::SetRunGen(int i){
41
42    fGen=i;
43
44}
45
46void NTraceLens::SetEffFactor(double ef){
47
48fEffectiveEfficiencyFactor=ef;
49
50}
51
52#define LARGECURV                1.0e10
53/*
54 * ========================================================================
55 * functions for diffractive optics start
56 * ========================================================================
57 */
58/*
59 * ------------------------------------------------------------------------
60 *  Calculate normal vector to the sphere approximating the lens surface
61 *  Input:
62 *    aX, aY, aZ: photon position on the lens surface
63 *    aCZ       : center of the approximated sphere
64 *
65 *  Output:
66 *    nnnDIFF[3]: normal unit vector
67 * ------------------------------------------------------------------------
68 */
69int NTraceLens::CalSpheVec(double aX, double aY, double aZ, double aCZ, double nnnDIFF[3])
70{
71
72  double N_;
73
74  nnnDIFF[0] = aX;
75  nnnDIFF[1] = aY;
76  nnnDIFF[2] = aZ - aCZ;
77
78  N_ = sqrt(vdot(nnnDIFF,nnnDIFF));
79
80  nnnDIFF[0] = nnnDIFF[0]/N_;
81  nnnDIFF[1] = nnnDIFF[1]/N_;
82  nnnDIFF[2] = nnnDIFF[2]/N_;
83
84  return 0;
85}
86
87
88/*
89 * ------------------------------------------------------------------------
90 *  read parameters of diffractive optics
91 *    Input:
92 *        dirname    the name of the directory name
93 *                        where data files are stored
94 *    Output:
95 *        (return value) status
96 * ------------------------------------------------------------------------
97 */
98int NTraceLens::ReadDiffractData(tel_param *param, FILE *fplog)
99{
100  FILE *fr;
101  char Infname[256];
102  double SD1_R,SD1_D,*tptr;
103
104#ifndef QUIET
105  if(fplog)
106    fprintf(fplog,"Reading Diffractive Surface Data (%s) ....  ",
107            param->surface[SURFDOE].filename);
108#endif
109
110  /* read diffractive optics data */
111  sprintf(Infname,"%.128s/%.127s",param->lens_dir,param->surface[SURFDOE].filename);
112
113  if(( fr = fopen(Infname,"r" )) == NULL)
114    {
115      if(fplog)
116        fprintf(fplog,"Cannot not OPEN %s\n",Infname );
117      return (-1);
118    }
119  param->surface[SURFDOE].size=NSURFPTS;
120  param->surface[SURFDOE].r=(double *)malloc(sizeof(double)*NSURFPTS);
121  param->surface[SURFDOE].z=(double *)malloc(sizeof(double)*NSURFPTS);
122  if(param->surface[SURFDOE].r==NULL || param->surface[SURFDOE].z==NULL)
123    {
124      if(fplog)
125        fprintf(fplog,"cannot allocate memory (DOE)\n");
126      fclose(fr);
127      return (-1);
128    }
129  param->surface[SURFDOE].ndata=0;
130
131  while(fscanf(fr,"%lf %lf",&SD1_R,&SD1_D)==2)
132    {
133      if(param->surface[SURFDOE].ndata >= param->surface[SURFDOE].size)
134        {
135          param->surface[SURFDOE].size+=NSURFPTS;
136          tptr=(double *)realloc(param->surface[SURFDOE].r,
137                                 sizeof(double)*param->surface[SURFDOE].size);
138          if(tptr==NULL)
139            {
140              if(fplog)
141                fprintf(fplog,"cannot allocate memory (DOE,r)\n");
142              fclose(fr);
143              free(param->surface[SURFDOE].r);
144              free(param->surface[SURFDOE].z);
145              param->surface[SURFDOE].r=NULL;
146              param->surface[SURFDOE].z=NULL;
147              param->surface[SURFDOE].size=0;
148              param->surface[SURFDOE].ndata=0;
149              return (-1);
150            }
151          else
152            {
153              param->surface[SURFDOE].r=tptr;
154              tptr=(double *)realloc(param->surface[SURFDOE].z,
155                                     sizeof(double)*param->surface[SURFDOE].size);
156              if(tptr==NULL)
157                {
158                  if(fplog)
159                    fprintf(fplog,"cannot allocate memory (DOE,z)\n");
160                  fclose(fr);
161                  free(param->surface[SURFDOE].r);
162                  free(param->surface[SURFDOE].z);
163                  param->surface[SURFDOE].r=NULL;
164                  param->surface[SURFDOE].z=NULL;
165                  param->surface[SURFDOE].size=0;
166                  param->surface[SURFDOE].ndata=0;
167                  return (-1);
168                }
169              else
170                param->surface[SURFDOE].z=tptr;
171            }
172        }
173      param->surface[SURFDOE].r[param->surface[SURFDOE].ndata]=SD1_R;
174      param->surface[SURFDOE].z[param->surface[SURFDOE].ndata]=SD1_D;
175      param->surface[SURFDOE].ndata++;
176    }
177  fclose(fr);
178
179#ifndef QUIET
180  if(fplog)
181    fprintf(fplog,"Done.\n");
182#endif
183
184  return (0);
185}
186
187/*
188 * ------------------------------------------------------------------------
189 *  calculate the diffractive structure
190 *   Input:
191 *      r   the distance in mm to the optical axis at the intersection point
192 *      s   surface ID (i.e. S1,S2,...)
193 *   Output:
194 *     (return value)
195 * ------------------------------------------------------------------------
196 */
197double NTraceLens::Cal_D(double r, tel_param *param)
198{
199  int i;
200
201  i=rt_bsearch(param->surface[SURFDOE].r,r,0,param->surface[SURFDOE].ndata-1);
202  return (param->surface[SURFDOE].z[i]);
203}
204
205/*
206 * ------------------------------------------------------------------------
207 *  calculate the diffracted ray
208 *  Input:
209 *        *photon       pointer to photon
210 *        nnn[3]        normal vector to the diffractive surface
211 *                         at the intersection
212 *        n1            refractive index of xxx
213 *        n2            refractive index of xxx
214 *        d
215 *        wl            wavelength in mm
216 *        m             index
217 *  Output:
218 *        *photon
219 * ------------------------------------------------------------------------
220 */
221int NTraceLens::diffract(int ID, tel_param *param, NPhoton *photon, double Zc,
222             double n1, double n2, double wl, double wl0, int m)
223{
224  int    i;
225  double CosIN, CosD, dir, nnnR[3], r3, ppp[3], nnn[3], d, norm;
226  double x3, y3, z3;
227
228  x3=photon->x[X];
229  y3=photon->x[Y];
230  z3=photon->x[Z];
231
232  if(fabs(param->surface[SURFDOE].Rc)>LARGECURV) /* plane */
233    {
234      nnn[0]= 0.0;
235      nnn[1]= 0.0;
236      nnn[2]= 1.0;
237    }
238  else /* sphere */
239    {
240      CalSpheVec(x3, y3, z3, Zc, nnn);
241    }
242
243  r3=sqrt(x3*x3+y3*y3);
244  if(r3<kTolerance){
245    return 0;
246  }
247  nnnR[0]= x3/r3;
248  nnnR[1]= y3/r3;
249  nnnR[2]=0.0;
250
251  d=Cal_D(r3,param);
252
253  if(photon->px[Z]>0.)
254    dir=1.0;
255  else
256    dir=-1.0;
257
258  CosIN=vdot(photon->px,nnn);
259  CosD=vdot(nnn,nnnR);
260
261  for(i=0;i<3;i++)
262    ppp[i]=(n1*(photon->px[i]-CosIN*nnn[i])
263            +(double)m*wl/wl0*d*(nnnR[i]-CosD*nnn[i])*dir)/n2;
264
265  norm=vdot(ppp,ppp);
266  if(norm>1.0)
267    return -1;
268
269  norm=sqrt(1.0-norm);
270
271#ifdef DEBUG
272  fprintf(stderr,"Diffractive: (%.3f , %.3f , %.3f)=>",
273          photon->px[X],photon->px[Y],photon->px[Z]);
274#endif
275
276  for(i=0;i<3;i++)
277    photon->px[i]=ppp[i]+norm*nnn[i]*(CosIN>0.0?(1.0):(-1.0));
278
279#ifdef DEBUG
280  fprintf(stderr,"(%.3f , %.3f , %.3f)\n",
281          photon->px[X],photon->px[Y],photon->px[Z]);
282#endif
283
284  if(m!=0){
285    if(param->DOE_efficiency<EsafRandom::Get()->Rndm()){
286      return (RT_ABSORBED);
287    }
288  }
289
290  return (0);
291}
292
293/*
294 * ========================================================================
295 * functions for diffractive optics END
296 * ========================================================================
297 */
298
299/*
300 * ------------------------------------------------------------------------
301 * calculate reflection probability assuming no polariation
302 *   n1,n2: refractive index of the media of this side and the other side
303 *   ci,ct: cosine of the incident light and reflected light to the surface
304 *          normal vector
305 * ------------------------------------------------------------------------
306 */
307static double ref(double n1, double n2, double ci, double ct)
308{
309  double rp,rs;
310
311  rp=(n2*ci-n1*ct)/(n2*ci+n1*ct);
312  rs=(n1*ci-n2*ct)/(n1*ci+n2*ct);
313
314  return((rp*rp+rs*rs)/2.0);
315}
316
317/*
318 * ------------------------------------------------------------------------
319 *  read lens surface shapes, and optical indices
320 * ------------------------------------------------------------------------
321 */
322int NTraceLens::ReadLensData(tel_param *param, FILE *fplog)
323{
324  FILE *fp;
325  char filename[256];
326  int i;
327  double tr,tsz,mat_wl,mat_n,mat_k,*tptr;
328
329  for(i=0;i<NSURFACES;i++)
330    {
331      param->surface[i].ndata=0;
332      param->surface[i].size=0;
333      if(!(i==1 || i==2 || i==4 || i==5 || i==6 || i==7)) continue;
334
335#ifndef QUIET
336  if(fplog)
337    fprintf(fplog,"Reading Surface(%d) Data (%s) ....  ",
338           i,param->surface[i].filename);
339#endif
340
341      sprintf(filename,"%.128s/%.128s",
342              param->lens_dir,param->surface[i].filename);
343      if((fp=fopen(filename,"r"))==NULL)
344        {
345          if(fplog)
346            fprintf(fplog,"%s not found.\n",filename);
347        return -1;
348        }
349
350      param->surface[i].size=NSURFPTS;
351      param->surface[i].r=(double *)malloc(sizeof(double)*NSURFPTS);
352      param->surface[i].z=(double *)malloc(sizeof(double)*NSURFPTS);
353      if(param->surface[i].r==NULL || param->surface[i].z==NULL)
354        {
355          if(fplog)
356            fprintf(fplog,"cannot allocate memory (S%d)\n",i);
357          fclose(fp);
358        }
359
360      while(fscanf(fp,"%lf %lf",&tr, &tsz)!=EOF)
361        {
362          if(param->surface[i].ndata >= param->surface[i].size)
363            {
364              param->surface[i].size+=NSURFPTS;
365              tptr=(double *)realloc(param->surface[i].r,
366                                     sizeof(double)*param->surface[i].size);
367              if(tptr==NULL)
368                {
369                  free(param->surface[i].r);
370                  free(param->surface[i].z);
371                  param->surface[i].r=NULL;
372                  param->surface[i].z=NULL;
373                  param->surface[i].size=0;
374                  param->surface[i].ndata=0;
375                  return (-1);
376                }
377              else
378                {
379                  param->surface[i].r=tptr;
380                  tptr=(double *)realloc(param->surface[i].z,
381                                         sizeof(double)*param->surface[i].size);
382                  if(tptr==NULL)
383                    {
384                      if(fplog)
385                        fprintf(fplog,"cannot allocate memory (S%d,z)\n",i);
386                      free(param->surface[i].r);
387                      free(param->surface[i].z);
388                      param->surface[i].r=NULL;
389                      param->surface[i].z=NULL;
390                      param->surface[i].ndata=0;
391                      fclose(fp);
392                      return -1;
393                    }
394                  else
395                    {
396                      param->surface[i].z=tptr;
397                    }
398                }
399            }
400          param->surface[i].r[param->surface[i].ndata]=tr;
401          param->surface[i].z[param->surface[i].ndata]=tsz;
402
403          param->surface[i].z[param->surface[i].ndata]=param->surface[i].Sz-param->surface[i].z[param->surface[i].ndata];
404          param->surface[i].ndata++;
405        }
406      fclose(fp);
407
408      if(param->surface[i].r[param->surface[i].ndata-1] < param->surface[i].r_lens)
409    param->surface[i].r_lens=param->surface[i].r[param->surface[i].ndata-1];
410
411#ifndef QUIET
412  if(fplog)
413    fprintf(fplog,"Done.\n");
414#endif
415    }
416
417  /* Read Material data */
418  for(i=2;i<param->nmat;i++)
419    {
420#ifndef QUIET
421  if(fplog)
422    fprintf(fplog,"Reading Optical Data, %s (%s) ....   ",
423            param->material[i].name, param->material[i].filename);
424#endif
425
426      sprintf(filename,"%.128s/%.128s",param->lens_dir,param->material[i].filename);
427      if((fp=fopen(filename,"r"))==NULL)
428        {
429          if(fplog)
430            fprintf(fplog,"Cannot find %s\n",filename);
431          return -10;
432        }
433
434      param->material[i].ndata=0;
435      param->material[i].size=NMATPTS;
436      param->material[i].wl=(double *)malloc(sizeof(double)*NMATPTS);
437      param->material[i].n=(double *)malloc(sizeof(double)*NMATPTS);
438      param->material[i].k=(double *)malloc(sizeof(double)*NMATPTS);
439      if(param->material[i].wl==NULL || param->material[i].n==NULL ||
440         param->material[i].k==NULL)
441        {
442          if(fplog)
443            fprintf(fplog,"cannot allocate memory (%s)\n",param->material[i].name);
444          fclose(fp);
445          return (-1);
446        }
447
448      while(fscanf(fp,"%lf %lf %lf",&mat_wl,&mat_n,&mat_k)!=EOF)
449        {
450          if(param->material[i].ndata >= param->material[i].size)
451            {
452              param->material[i].size+=NMATPTS;
453              tptr=(double *)realloc(param->material[i].wl,
454                                     sizeof(double)*param->material[i].size);
455              if(tptr==NULL)
456                {
457                  if(fplog)
458                    fprintf(fplog,"cannot allocate memory (%s,wl)\n",param->material[i].name);
459                  fclose(fp);
460                  free(param->material[i].wl);
461                  free(param->material[i].n);
462                  free(param->material[i].k);
463                  param->material[i].wl=NULL;
464                  param->material[i].n=NULL;
465                  param->material[i].k=NULL;
466                  param->material[i].size=0;
467                  param->material[i].ndata=0;
468                  return (-1);
469                }
470              else
471                {
472                  param->material[i].wl=tptr;
473                  tptr=(double *)realloc(param->material[i].n,
474                                         sizeof(double)*param->material[i].size);
475                  if(tptr==NULL)
476                    {
477                      if(fplog)
478                        fprintf(fplog,"cannot allocate memory (%s,n)\n",param->material[i].name);
479                      fclose(fp);
480                      free(param->material[i].wl);
481                      free(param->material[i].n);
482                      free(param->material[i].k);
483                      param->material[i].wl=NULL;
484                      param->material[i].n=NULL;
485                      param->material[i].k=NULL;
486                      param->material[i].size=0;
487                      param->material[i].ndata=0;
488                      return (-1);
489                    }
490                  else
491                    {
492                      param->material[i].n=tptr;
493
494                      tptr=(double *)realloc(param->material[i].k,
495                                             sizeof(double)*param->material[i].size);
496                      if(tptr==NULL)
497                        {
498                          fclose(fp);
499                          free(param->material[i].wl);
500                          free(param->material[i].n);
501                          free(param->material[i].k);
502                          param->material[i].wl=NULL;
503                          param->material[i].n=NULL;
504                          param->material[i].k=NULL;
505                          param->material[i].size=0;
506                          param->material[i].ndata=0;
507                          return (-1);
508                        }
509                      else
510                        param->material[i].k=tptr;
511                    }
512                }
513            }
514          param->material[i].wl[param->material[i].ndata]=mat_wl;
515          param->material[i].n[param->material[i].ndata]=mat_n;
516          param->material[i].k[param->material[i].ndata]=mat_k;
517          param->material[i].ndata++;
518        }
519      fclose(fp);
520
521      param->material[i].zn=(double *)malloc(sizeof(double)*param->material[i].ndata);
522      param->material[i].zk=(double *)malloc(sizeof(double)*param->material[i].ndata);
523      if(param->material[i].zn==NULL || param->material[i].zk==NULL)
524        {
525          if(fplog)
526            fprintf(fplog,"Canoot allocate memory for optical data (%s)\n",
527                    param->material[i].name);
528          return (-10);
529        }
530      csp_maketable(param->material[i].wl, param->material[i].n,
531                    param->material[i].zn, param->material[i].ndata);
532      csp_maketable(param->material[i].wl, param->material[i].k,
533                    param->material[i].zk, param->material[i].ndata);
534#ifndef QUIET
535  if(fplog)
536    fprintf(fplog,"Done.\n");
537#endif
538    }
539
540  return 0;
541}
542
543/*
544 * ------------------------------------------------------------------------
545 *  get refractive index of material at wavelength lambda
546 *    using cubic spline interpolation
547 * ------------------------------------------------------------------------
548 */
549double NTraceLens::Get_n(int matid, tel_param *param, double lambda)
550{
551  if(matid==0 || matid==1)
552    return (*(param->material[matid].n));
553  else if(matid<param->nmat)
554    return cspline(lambda, param->material[matid].wl,
555                   param->material[matid].n, param->material[matid].zn,
556                   param->material[matid].ndata);
557  else
558    return 0.;
559}
560
561/*
562 * ------------------------------------------------------------------------
563 *  get extinction coefficient of material at wavelength lambda
564 *    using cubic spline interpolation
565 * ------------------------------------------------------------------------
566 */
567double NTraceLens::Get_k(int matid, tel_param *param, double lambda)
568{
569  if(matid==0 || matid==1)
570    return (*(param->material[0].k));
571  else if(matid<param->nmat)
572    return cspline(lambda, param->material[matid].wl,
573                   param->material[matid].k, param->material[matid].zk,
574                   param->material[matid].ndata);
575  else
576    return 0.;
577}
578/*
579 * ------------------------------------------------------------------------
580 *  get the refracted ray
581 *    if reflected it is thrown away.
582 *
583 *  *photon: info on photon (only incident photon direction is used)
584 *   nnn[3]  : surface normal vector (direct to this side)
585 *   n1,n2   : refractive indices of the media of this side and the other side
586 * ------------------------------------------------------------------------
587 */
588int NTraceLens::refract(NPhoton *photon, double nnn[3],double n1,double n2)
589{
590  int i, TotalReflection=0;
591  double CosIN, CosOUT=1.0, ppp[3], norm, out_phdir[3];
592
593  CosIN=vdot(photon->px,nnn);
594  for(i=0;i<3;i++)
595    ppp[i]=(photon->px[i]-CosIN*nnn[i])*n1/n2;
596
597  norm=vdot(ppp,ppp);
598#ifdef NO_REFLECTION
599  if(norm>1.0)
600    return(RT_NO_REFRACTED_LIGHT);  /* total reflection */
601#else
602  if(norm>1.0)
603    TotalReflection=1;
604  else
605    {
606#endif
607
608      norm=sqrt(1.0-norm);
609      for(i=0;i<3;i++)
610        out_phdir[i]=ppp[i]+norm*nnn[i]*(CosIN>0.0?(1.0):(-1.0));
611
612#ifndef NO_REFLECTION
613      /* check reflection */
614      CosOUT=vdot(out_phdir,nnn);
615    }
616
617  if(TotalReflection==1 || fGen==1?EsafRandom::Get()->Rndm():RANDOM() < ref(n1,n2,CosIN,CosOUT))
618    {
619#ifdef DEBUG
620      fprintf(stderr,"Reflection: (%.3f,%.3f,%.3f)=>",
621             photon->px[X],photon->px[Y],photon->px[Z]);
622#endif
623      for(i=0;i<3;i++)
624        photon->px[i]-=2.0*CosIN*nnn[i];
625#ifdef DEBUG
626      fprintf(stderr,"(%.3f,%.3f,%.3f)\n",
627             photon->px[X],photon->px[Y],photon->px[Z]);
628#endif
629      return(RT_REFLECTED);
630    }
631#endif
632
633  /* photon->px[]: direction of incident ray -> refracted ray */
634#ifdef DEBUG
635      fprintf(stderr,"Refraction: (%.3f,%.3f,%.3f)=>",
636             photon->px[X],photon->px[Y],photon->px[Z]);
637#endif
638
639  for(i=0;i<3;i++)
640    photon->px[i]=out_phdir[i];
641
642#ifdef DEBUG
643  fprintf(stderr,"(%.3f,%.3f,%.3f)\n",
644          photon->px[X],photon->px[Y],photon->px[Z]);
645  fprintf(stderr,"normal vect=(%.3f,%.3f,%.3f)\n",nnn[0],nnn[1],nnn[2]);
646#endif
647
648  return(RT_REFRACTED);
649}
650
651/*
652 * ------------------------------------------------------------------------
653 *  calculate the differential coefficient at lens_r[i]
654 *   Input:
655 *      i                : the radius index of lens_r[i]
656 *      lens_r[],lens_z[]: lens surface data sets of the radius and the height
657 *      n_lens_elem      : number of data elements of lens_r[]
658 *   Output:
659 *     (return value)    : differential coefficient at lens_r[i]
660 * ------------------------------------------------------------------------
661 */
662double NTraceLens::differential_r(int i, double *lens_r,double *lens_z, int n_lens_elem)
663{
664  double dr1,dr2;
665  double dz1,dz2;
666  int j1,j2;
667
668  dz1=dz2=0.;
669  j1=i+1;
670  j2=i-1;
671
672  if(i>=n_lens_elem-1)
673    return 0.0;
674  else if(i==n_lens_elem-2)
675    {
676      if(lens_r[i+1]-lens_r[i]<0.001)
677        return 0.0;
678      else
679        return (lens_z[i+1]-lens_z[i])/(lens_r[i+1]-lens_r[i]);
680    }
681
682  dr1=lens_r[j1]-lens_r[i];
683  if(dr1<0.001)
684    {
685      dz1+=lens_z[j1]-lens_z[i];
686      j1++;
687    }
688
689  if(j2<0)
690    {
691      j2=-j2;
692      dr2=lens_r[i]+lens_r[j2];
693    }
694  else
695    {
696      dr2=lens_r[i]-lens_r[j2];
697    }
698  if(dr2<0.001)
699    {
700      dz2+=lens_z[j2]-lens_z[i];
701      j2--;
702    }
703
704  return (lens_z[j1]-dz1-lens_z[j2]+dz2)/(dr1+dr2);
705}
706
707/*
708 *
709 * check if there are any intersections with a cylinder
710 *
711 */
712static int check_intersect_cylinder(double r1, double z1, double r2, double z2,
713                                    NPhoton *photon, double *dis)
714{
715  double aa[3], xx[2], tmp, z;
716  int    nroot, i;
717
718  aa[0] = (photon->px[X])*(photon->px[X]) + (photon->px[Y])*(photon->px[Y]);
719  aa[1] = 2.0*((photon->x[X])*(photon->px[X]) + (photon->x[Y])*(photon->px[Y]));
720  aa[2] = ((photon->x[X])*(photon->x[X]) + (photon->x[Y])*(photon->x[Y])) - r1*r1;
721
722  nroot=SolveQuadratic(aa,xx);
723  if(nroot<2)
724    return(RT_NO_INTERSECTION1);
725
726  /* solution is the smaller positive value */
727  /* swap */
728  if(xx[0]>xx[1])
729    {
730      tmp=xx[0];
731      xx[0]=xx[1];
732      xx[1]=tmp;
733    }
734  if(xx[1]<0)
735    return(RT_NO_INTERSECTION1);
736
737  for(i=0;i<2;i++)
738    {
739      if(xx[i]>=0.0)
740        {
741          z=photon->x[Z]+photon->px[Z]*xx[i];
742          /* check if the solution is in the specified range */
743          if((z-z1)*(z-z2)<0)
744            {
745              *dis=xx[i];
746              return 0;
747            }
748        }
749    }
750
751  *dis=(xx[0]>0 ? xx[0]:xx[1]);
752  return(RT_NO_INTERSECTION2);
753}
754
755/*
756 *
757 * check if there are any intersections with a cone
758 *
759 */
760static int check_intersect_cone(double r1, double z1, double r2, double z2,
761                                NPhoton *photon, double *dis)
762{
763  double a00, a0, aa[3], xx[2], z0, dbuf1, tmp, x, y, r, zz[2];
764  int    nroot, i;
765
766  a00 = (z2-z1)/(r2-r1);
767  z0  = z1-a00*r1;
768  a0  = a00*a00;
769
770  dbuf1 = (photon->x[Z])-z0;
771  aa[0] = a0*(photon->px[X])*(photon->px[X]) + a0*(photon->px[Y])*(photon->px[Y] )- (photon->px[Z])*(photon->px[Z]);
772  aa[1] = 2.0*(a0*(photon->px[X])*(photon->x[X]) + a0*(photon->px[Y])*(photon->x[Y]) - (photon->px[Z])*dbuf1);
773  aa[2] = a0*(photon->x[X])*(photon->x[X]) + a0*(photon->x[Y])*(photon->x[Y]) - dbuf1*dbuf1;
774
775  nroot=SolveQuadratic(aa,xx);
776  if(nroot<2)
777    return(RT_NO_INTERSECTION1);
778
779  /* solution is the smaller positive value */
780  /* swap */
781  if(xx[0]>xx[1])
782    {
783      tmp=xx[0];
784      xx[0]=xx[1];
785      xx[1]=tmp;
786    }
787  if(xx[1]<0.0)
788    return(RT_NO_INTERSECTION1);
789
790  /* We usually have a solution on a virtual cone with -a00 */
791  /* check this condition */
792  for(i=0;i<2;i++)
793    if(xx[i]>=0.0)
794      {
795        x=photon->px[X]*xx[i]+photon->x[X];
796        y=photon->px[Y]*xx[i]+photon->x[Y];
797        zz[i]=photon->px[Z]*xx[i]+photon->x[Z];
798        r=sqrt(x*x+y*y);
799        if((r-r1)*(r-r2)<0.0)
800          {
801            *dis=xx[i];
802            return 0;
803          }
804      }
805
806  if(xx[0]<0.0)
807    *dis=xx[1];
808  else
809    {
810      /* the nearer one to the surface will be what we want... */
811      if(fabs(zz[0]-z1)<fabs(zz[1]-z1))
812        *dis=xx[0];
813      else
814        *dis=xx[1];
815    }
816
817  return(RT_NO_INTERSECTION2);
818}
819
820/*
821 * ------------------------------------------------------------------------
822 *
823 *  ID:            Lens Surface ID
824 *  phot[8]: input and output photon data
825 *    phot[0,1,2]: x,y,z position in mm
826 *    phot[3,4,5]: l,m,n direction
827 *    phot[6]    : wavelength in nm
828 *    phot[7]    : time in ns
829 *  offset_z     : the offset z-position of sphere center approximating
830 *                 lens surface
831 *  radius       : the radius of sphere approximated lens surface
832 *  lens_r[],lens_z[]: lens surface data sets of the radius and the height
833 *  n_lens_elem  : number of data elements of lens_r[]
834 *  medium1_n, medium2_n: refractive indices of the media in this side
835 *                         and the other side
836 * ------------------------------------------------------------------------
837 */
838int lens_interaction(int ID, NPhoton *photon, tel_param *param,
839                     double medium1_n, double medium2_n, double *dis)
840{
841  int r_idx,iter=0,flag,Found,nroot,dir,status;
842  double dbuf1, aa[3], xx[2], x3,y3,z3,r3, r1,r2, z1,z2;
843  double nnn[3];
844  double a00;
845  double tmp, v[3];
846  double Zc, Rc, lens_r_max,  *lens_r, *lens_z;
847  int n_lens_elem;
848#ifdef DEBUG
849  double approx_x[3];
850#endif
851
852  Zc=param->surface[ID].Zc;
853  Rc=param->surface[ID].Rc;
854  lens_r_max=param->surface[ID].r_lens;
855  lens_r=param->surface[ID].r;
856  lens_z=param->surface[ID].z;
857  n_lens_elem=param->surface[ID].ndata;
858
859  if(photon->px[Z]==0)
860    return(RT_HIT_SIDE_WALL);
861
862  if(fabs(Rc)>LARGECURV)
863    {
864      dbuf1 = (param->surface[ID].Sz-photon->x[Z])/photon->px[Z];
865      x3=photon->px[X]*dbuf1+photon->x[X];
866      y3=photon->px[Y]*dbuf1+photon->x[Y];
867      z3=photon->px[Z]*dbuf1+photon->x[Z];
868      r3=sqrt(x3*x3+y3*y3);
869      nnn[0]=0.0;
870      nnn[1]=0.0;
871      nnn[2]=1.0;
872      if(r3>lens_r_max)
873        return(RT_HIT_SIDE_WALL);
874      *dis=dbuf1;
875    }
876  else
877    {
878  /* as an initial guess, intersection with an approximated sphere
879     is calculated */
880
881      dbuf1 = photon->x[Z]-Zc;
882      aa[0] = photon->px[X]*photon->px[X] + photon->px[Y]*photon->px[Y] + photon->px[Z]*photon->px[Z];
883      aa[1] = 2.0*(photon->px[X]*photon->x[X] + photon->px[Y]*photon->x[Y] + photon->px[Z]*dbuf1);
884      aa[2] = photon->x[X]*photon->x[X] + photon->x[Y]*photon->x[Y] + dbuf1*dbuf1-Rc*Rc;
885
886      nroot=SolveQuadratic(aa,xx);
887      if(nroot<2)
888        return(RT_NO_INTERSECTION1);
889
890      /* swap */
891      if(xx[0]>xx[1])
892        {
893          tmp=xx[0];
894          xx[0]=xx[1];
895          xx[1]=tmp;
896        }
897
898      /* We need a positive, smaller solution. */
899      if(xx[0]>0.0)
900        dbuf1=xx[0];
901      else
902        {
903          if(xx[1]>0.0)
904            dbuf1=xx[1];
905          else
906            return(RT_NO_INTERSECTION1);
907        }
908
909      x3=photon->px[X]*dbuf1+photon->x[X];
910      y3=photon->px[Y]*dbuf1+photon->x[Y];
911      z3=photon->px[Z]*dbuf1+photon->x[Z];
912      /* distance to the optical axis at the initially guessed intersection */
913      r3=sqrt(x3*x3+y3*y3);
914
915      if(r3>lens_r_max && xx[0]>-100.0) /* -100 is not optimized */
916        {
917          /* The other solution may be acceptable just for an initial guess. */
918          dbuf1=xx[0];
919          x3=photon->px[X]*dbuf1+photon->x[X];
920          y3=photon->px[Y]*dbuf1+photon->x[Y];
921          z3=photon->px[Z]*dbuf1+photon->x[Z];
922
923          r3=sqrt(x3*x3+y3*y3);
924          if(r3>lens_r_max)
925            return(RT_HIT_SIDE_WALL);
926        }
927      *dis=dbuf1;
928
929#ifdef DEBUG
930      fprintf(stderr,"1st est.: %.3f ",r3);
931      approx_x[0]=x3;approx_x[1]=y3;approx_x[2]=z3;
932#endif
933    }
934  {
935      /* search taking account of Fresnel structure */
936      /* search the index near r3 */
937      r_idx=rt_bsearch(lens_r,r3,0,n_lens_elem-1);
938
939#ifdef DEBUG
940      fprintf(stderr,"%.3f\n",r3);
941#endif
942      if(r_idx>n_lens_elem-2)
943        r_idx=n_lens_elem-2;
944
945      /* search intersection point with real fresnel lens */
946      iter=0; /* number of iterations */
947      do{
948        Found=0;
949
950        r1=lens_r[r_idx];
951        z1=lens_z[r_idx];
952        r2=lens_r[r_idx+1];
953        z2=lens_z[r_idx+1];
954
955        if(r1==r2)
956          {
957            /* intersection with a part of the cylinder */
958            status=check_intersect_cylinder(r1, z1, r2, z2, photon, dis);
959
960            x3=photon->px[X]*(*dis)+photon->x[X];
961            y3=photon->px[Y]*(*dis)+photon->x[Y];
962            z3=photon->px[Z]*(*dis)+photon->x[Z];
963
964            if(!status)     /* intersection found */
965              {
966                Found=1;
967#ifdef DEBUG
968                fprintf(stderr,"Hit Backcut\n");
969                /* return (RT_ABSORBED);*/
970#endif
971                /* normal vector */
972                if(z2>z1)
973                  dir=1;
974                else
975                  dir=-1;
976                nnn[X]=dir*x3/r1;
977                nnn[Y]=dir*y3/r1;
978                nnn[Z]=0.0;
979              }
980            else /* intersection with cylinder not found */
981              {
982                if(status==RT_NO_INTERSECTION2 && fabs(z3-z1)<100.0)
983                  {
984                    v[X]=x3;
985                    v[Y]=y3;
986                    v[Z]=0.0;
987
988                    if((z3-z1)*vdot(v,&photon->px[X])*photon->px[Z]>0)
989                      r_idx--;
990                    else
991                      r_idx++;
992                  }
993                else
994                  {
995                    *dis=(z1-photon->x[Z])/photon->px[Z];
996                    z3=z1;
997                    x3=photon->px[X]*(*dis)+photon->x[X];
998                    y3=photon->px[Y]*(*dis)+photon->x[Y];
999                    r3=sqrt(x3*x3+y3*y3);
1000                    r_idx=rt_bsearch(lens_r,r3,0,n_lens_elem-1);
1001                  }
1002              }
1003#ifdef DEBUG
1004            fprintf(stderr,"ridx0(%d): %d (%.1f: %.3f %.3f %.3f (%.1f %.1f %.1f)(%.3f %.3f %.3f))\n",
1005                    status,r_idx,z1,r1,r2,r3,photon->x[X],photon->x[Y],photon->x[Z],photon->px[X],photon->px[Y],photon->px[Z]);
1006#endif
1007          }
1008        else  /* r1 != r2 */
1009          {
1010            /* intersection with a part of the cone */
1011            if(z1==z2)
1012              {
1013                *dis = (z1-photon->x[Z])/photon->px[Z];
1014                status=0;
1015              }
1016            else
1017              status=check_intersect_cone(r1, z1, r2, z2, photon, dis);
1018
1019            x3=photon->px[X]*(*dis)+photon->x[X];
1020            y3=photon->px[Y]*(*dis)+photon->x[Y];
1021            z3=photon->px[Z]*(*dis)+photon->x[Z];
1022
1023            r3=sqrt(x3*x3+y3*y3);
1024
1025            if(!status)   /* intersection found */
1026              {
1027                Found=1;
1028                /* normal vector */
1029                if(r3<kTolerance){
1030                  nnn[0]=0.;
1031                  nnn[1]=0.;
1032                  nnn[2]=1.;
1033                }else{
1034                  if(r_idx>n_lens_elem-2)
1035                    r_idx=n_lens_elem-2;
1036                  a00=(differential_r(r_idx,lens_r,lens_z,n_lens_elem)
1037                       *(r3-lens_r[r_idx])-
1038                       differential_r(r_idx+1,lens_r,lens_z,n_lens_elem)
1039                       *(r3-lens_r[r_idx+1]))/(lens_r[r_idx+1]-lens_r[r_idx]);
1040                  nnn[X]=a00*x3;
1041                  nnn[Y]=a00*y3;
1042                  nnn[Z]=-r3;
1043                  dbuf1=sqrt(vdot(nnn,nnn));
1044                  nnn[X]/=dbuf1;
1045                  nnn[Y]/=dbuf1;
1046                  nnn[Z]/=dbuf1;
1047                }
1048              }
1049            else /* not found */
1050              {
1051                if(status==RT_NO_INTERSECTION2 &&
1052                   fabs(z3-z1)<10. && fabs(r3-r1)<5.)
1053                  {
1054                    if(r3<r1)
1055                      r_idx--;
1056                    if(r3>r2)
1057                      r_idx++;
1058                  }
1059                else
1060                  {
1061                    *dis=(z1-photon->x[Z])/photon->px[Z];
1062                    z3=z1;
1063                    x3=photon->px[X]*(*dis)+photon->x[X];
1064                    y3=photon->px[Y]*(*dis)+photon->x[Y];
1065                    r3=sqrt(x3*x3+y3*y3);
1066                    r_idx=rt_bsearch(lens_r,r3,0,n_lens_elem-1);
1067                  }
1068              }
1069#ifdef DEBUG
1070            fprintf(stderr,"ridx1(%d): %d (%.1f: %.3f %.3f %.3f (%.1f %.1f %.1f)(%.3f %.3f %.3f))\n",
1071                    status,r_idx,z1,r1,r2,r3,photon->x[X],photon->x[Y],photon->x[Z],photon->px[X],photon->px[Y],photon->px[Z]);
1072#endif
1073          }
1074
1075        if(r_idx>=n_lens_elem)
1076          return(RT_OUT_OF_SURFACE2);
1077        if(r_idx<0 && r3>lens_r_max && fabs(z3-z1)<10.0)
1078          return(RT_OUT_OF_SURFACE2);
1079        if(r_idx<0)
1080          {
1081            r_idx=1;
1082            /* return(RT_OUT_OF_SURFACE1); */
1083          }
1084
1085        /* iteration doesn't converge */
1086        if(iter>RT_SEARCH_ITER_MAX)
1087          {
1088#ifdef DEBUG
1089            fprintf(stderr,"app.(%.1f,%.1f,%.1f)=>end(%.1f,%.1f,%.1f)\n",
1090                    approx_x[0],approx_x[1],approx_x[2],x3,y3,z3);
1091#endif
1092            return(RT_EXCEED_ITER_MAX);
1093          }
1094
1095        iter++;
1096
1097      } while(!Found);
1098    }
1099
1100  /* refraction */
1101  if(photon->px[Z]<0.0)
1102    {
1103      tmp=medium2_n;
1104      medium2_n=medium1_n;
1105      medium1_n=tmp;
1106    }
1107#ifdef DEBUG
1108  fprintf(stderr,"ridx=%d n_med1=%.5f n_med2=%.5f\n",r_idx,medium1_n,medium2_n);
1109#endif
1110  flag=refract(photon, nnn, medium1_n, medium2_n);
1111#if 0
1112  if(flag==RT_REFLECTED)
1113    return(RT_REFLECTION_DISCARDED); /* reflected photons are thrown away */
1114#endif
1115
1116  /* distance between injection point and the Surface */
1117  photon->t+=(*dis)* medium1_n/EConst::Clight();
1118
1119  photon->x[X]=x3;
1120  photon->x[Y]=y3;
1121  photon->x[Z]=z3;
1122
1123  return flag;
1124}
1125
1126int NTraceLens::CheckRcut(tel_param *param, NPhoton *photon, int flag)
1127{
1128  int surfid;
1129  surfid=(-flag)/1000;
1130  if(surfid>8)
1131    surfid=8;
1132
1133  if(fabs(photon->x[Y])>param->surface[surfid].r_cut)
1134    return RT_OUT_OF_SURFACE1;
1135  else if(fabs(photon->x[X])>param->surface[surfid].r_cutx)
1136    return RT_OUT_OF_SURFACE1;
1137  else if(photon->x[X]*photon->x[X]+photon->x[Y]*photon->x[Y]
1138          >param->surface[surfid].r_lens*param->surface[surfid].r_lens) /* r_lens should be r_wall ...? */
1139    return RT_OUT_OF_SURFACE1;
1140  else
1141    return 0;
1142}
1143
1144static int CheckIris(tel_param *param, NPhoton *photon, int flag)
1145{
1146  double x3, y3, z3=0.0, r3=0.0, dist, radius;
1147  int surfid;
1148
1149  surfid=-flag/1000;
1150
1151  z3=param->surface[surfid].Sz;
1152  r3=param->surface[surfid].r_lens;
1153
1154  dist=(z3 - photon->x[Z])/photon->px[Z];
1155  if(dist>=0.0)
1156    {
1157      x3=photon->px[X]*dist+photon->x[X];
1158      y3=photon->px[Y]*dist+photon->x[Y];
1159
1160      photon->x[X]=x3;
1161      photon->x[Y]=y3;
1162      photon->x[Z]=z3;
1163      photon->t+=dist*param->material[param->surface[3].matid].n[0]/EConst::Clight();
1164
1165      radius=sqrt(x3*x3+y3*y3);
1166      if(radius > param->r_wall)
1167        return RT_HIT_SIDE_WALL;
1168      else if(radius > r3)
1169        {
1170          return RT_OUT_OF_IRIS;
1171        }
1172    }
1173
1174    return 0;
1175}
1176
1177/* -------------------------------------------------------------------------
1178 * int trace_lens(tel_param *param, NPhoton *photon, FILE *fp)
1179 *
1180 *   --- trace photons in the EUSO lens system (the first lens to the last one)
1181 *
1182 * Input:
1183 *        param        tel_param* to store the optics parameters
1184 *                       (see also in tracemain_optF1v2.h)
1185 *        photon       info on photon (position[mm], direction, time[ns],wavelength[m]
1186 *        fp           FILE* for (error) message logging
1187 *
1188 * Output:
1189 *  return value
1190 *        0             normal end
1191 *       <0             error
1192 *                      (see trace_optF1v2.h for more detail
1193 *                        from RT_NO_SPINTERSECTION to RT_ESCAPE_FROM_ENTRANCE)
1194 * -------------------------------------------------------------------------
1195 */
1196
1197
1198
1199//#define DEBUG
1200int NTraceLens::trace_lens(tel_param *param, NPhoton *photon, FILE *fplog)
1201{
1202  int    status=0, n_int, WithinOptics, i, matid;
1203  int    PrevPoint, CurPoint, NextPoint;
1204  double nn[8],kk[8];
1205  double AbsorbCoeff, TransProb;
1206  double dis;
1207
1208  /* Opt. index  for this lambda */
1209  for(i=0;i<param->nmat;i++)
1210    {
1211      nn[i]=Get_n(i,param,photon->lambda);
1212      kk[i]=Get_k(i,param,photon->lambda);
1213    }
1214
1215  if(photon->x[Z]>=param->surface[7].Sz && photon->px[Z]<0.0)
1216    {
1217      PrevPoint=SURFACE_FS;
1218      CurPoint=SURFACE_FS;
1219      NextPoint=SURFACE7;
1220    }
1221  else
1222    {
1223      PrevPoint=SURFACE0;
1224      CurPoint=SURFACE0;
1225      NextPoint=SURFACE1;
1226    }
1227  n_int=0;
1228  WithinOptics=1;
1229
1230#ifdef DEBUG
1231  fprintf(stderr,"0 %.1f %.1f %.1f  %.3f %.3f %.3f\n",
1232          photon->x[X],photon->x[Y],photon->x[Z],
1233          photon->px[X],photon->px[Y],photon->px[Z]);
1234#else
1235
1236
1237  if(param->flag_printray)
1238    printf("0 %.1f %.1f %.1f\n",photon->x[X],photon->x[Y],photon->x[Z]);
1239#endif
1240
1241  while(n_int++ < RT_INTERACTION_MAX && WithinOptics)
1242    {
1243      PrevPoint=CurPoint;
1244      CurPoint=NextPoint;
1245
1246      switch(CurPoint)
1247        {
1248        case SURFACE0:
1249        case SURFACE_FS:
1250          WithinOptics=0;
1251          break;
1252
1253        case SURFACE1:
1254          status=lens_interaction(1, photon, param,
1255                                  nn[param->surface[1].matid],
1256                                  nn[param->surface[2].matid],
1257                                  &dis);
1258#ifdef DEBUG
1259      printf("1 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1260
1261#endif
1262          if(status==RT_REFLECTED)
1263            {
1264              NextPoint=PrevPoint;
1265              status=0;
1266            }
1267          else if(PrevPoint==SURFACE0)
1268            NextPoint=SURFACE2;
1269          else
1270            NextPoint=SURFACE0;
1271          break;
1272
1273        case SURFACE2:
1274          status=lens_interaction(2, photon, param,
1275                                  nn[param->surface[2].matid],
1276                                  nn[param->surface[3].matid],
1277                                  &dis);
1278#ifdef DEBUG
1279      printf("2 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1280
1281#endif
1282          if(status==RT_REFLECTED)
1283            {
1284              NextPoint=PrevPoint;
1285              status=0;
1286            }
1287          else if(PrevPoint==SURFACE1)
1288            NextPoint=SURFACE3;
1289          else
1290            NextPoint=SURFACE1;
1291          break;
1292
1293        case SURFACE3:
1294          if((status=CheckIris(param,photon,SURFACE3))<0){
1295            return (status+SURFACE3);}
1296          if(photon->px[Z]>0.0)
1297            NextPoint=SURFACE4;
1298          else
1299            NextPoint=SURFACE2;
1300          break;
1301
1302        case SURFACE4:
1303          status=lens_interaction(4, photon, param,
1304                                  nn[param->surface[4].matid],
1305                                  nn[param->surface[5].matid],
1306                                  &dis);
1307#ifdef DEBUG
1308      printf("4 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1309
1310#endif
1311          if(status==RT_REFLECTED)
1312            {
1313              NextPoint=PrevPoint;
1314              status=0;
1315            }
1316          else if(photon->px[Z]>0.0)
1317            NextPoint=SURFACE5;
1318          else
1319            NextPoint=SURFACE3;
1320          break;
1321
1322        case SURFACE5:
1323#ifndef NO_DIFFRACTIVE
1324          if(photon->px[Z]<0.0)
1325            {
1326              status=diffract(5, param, photon,
1327                              param->surface[5].Zc, 1.0003, 1.0003,
1328                              photon->lambda*1e-6, param->lambda0*1e-6, param->DOE_order);
1329              if(status!=0)
1330                break;
1331            }
1332#endif
1333          status=lens_interaction(5, photon, param,
1334                                  nn[param->surface[5].matid],
1335                                  nn[param->surface[6].matid],
1336                                  &dis);
1337#ifdef DEBUG
1338      printf("5 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1339
1340#endif
1341#ifndef NO_DIFFRACTIVE
1342          if(photon->px[Z]>0.0)
1343            {
1344              status=diffract(5, param, photon,
1345                              param->surface[5].Zc, 1.0003, 1.0003,
1346                              photon->lambda*1e-6, param->lambda0*1e-6, param->DOE_order);
1347              if(status!=0)
1348                break;
1349            }
1350#endif
1351
1352          if(status==RT_REFLECTED)
1353            {
1354              NextPoint=PrevPoint;
1355              status=0;
1356            }
1357          else if(photon->px[Z]>0.0)
1358            NextPoint=SURFACE6;
1359          else
1360            NextPoint=SURFACE4;
1361          break;
1362
1363        case SURFACE6:
1364          status=lens_interaction(6, photon, param,
1365                                  nn[param->surface[6].matid],
1366                                  nn[param->surface[7].matid],
1367                                  &dis);
1368#ifdef DEBUG
1369      printf("6 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1370
1371#endif
1372          if(status==RT_REFLECTED)
1373            {
1374              NextPoint=PrevPoint;
1375              status=0;
1376            }
1377          else if(PrevPoint==SURFACE5)
1378            NextPoint=SURFACE7;
1379          else
1380            NextPoint=SURFACE5;
1381          break;
1382
1383        case SURFACE7:
1384          status=lens_interaction(7, photon, param,
1385                                  nn[param->surface[7].matid],
1386                                  nn[param->surface[8].matid],
1387                                  &dis);
1388#ifdef DEBUG
1389      printf("7 %f %f %f stat=%i\n",photon->x[X],photon->x[Y],photon->x[Z], status);
1390
1391#endif
1392          if(status==RT_REFLECTED)
1393            {
1394              NextPoint=PrevPoint;
1395              status=0;
1396            }
1397          else if(PrevPoint==SURFACE6)
1398            NextPoint=SURFACE_FS;
1399          else
1400            NextPoint=SURFACE6;
1401          break;
1402        }
1403//      if(WithinOptics)
1404//      StorePhotonHistory(photon);
1405#ifndef DEBUG
1406      if(param->flag_printray && CurPoint>SURFACE_FS)
1407        printf("%d %.1f %.1f %.1f\n",
1408               -CurPoint/1000,photon->x[X],photon->x[Y],photon->x[Z]);
1409#else
1410      fprintf(stderr,"%d %.1f %.1f %.1f  %.3f %.3f %.3f\n",
1411              -CurPoint/1000,
1412              photon->x[X],photon->x[Y],photon->x[Z],
1413              photon->px[X],photon->px[Y],photon->px[Z]);
1414#endif
1415      if(status<0){
1416        return (status+CurPoint);}
1417      status=CheckRcut(param,photon,CurPoint);
1418      if(status<0){
1419        return (status+CurPoint);}
1420
1421#ifndef NO_ABSORPTION
1422      /* absorption between S1 and S2, S4 and S5, or between S6 and S7 */
1423
1424      if(PrevPoint>CurPoint)
1425        matid=param->surface[-CurPoint/1000].matid;
1426      else
1427        matid=param->surface[-CurPoint/1000+1].matid;
1428
1429      if(matid>1)
1430        {
1431          AbsorbCoeff=kk[matid]*M_PI*4.0/(photon->lambda*1.e-6);
1432          TransProb=exp(-AbsorbCoeff*dis);
1433
1434          if(fGen==1?EsafRandom::Get()->Rndm():RANDOM()>TransProb){
1435          return(RT_ABSORBED+CurPoint);}
1436        }
1437#endif
1438    }
1439
1440  if(n_int>=RT_INTERACTION_MAX){
1441      return (RT_STRAY_LIGHT);}
1442
1443
1444  return 0;
1445}
1446
1447
1448
Note: See TracBrowser for help on using the repository browser.