source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/Ntrace_optics.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: 12.2 KB
Line 
1/* -------------------------------------------------------------------------
2 *   trace_optics.c
3 *
4 *   ---  main program for raytracing
5 *
6 * Copyright (c) 2000-2010 N.Sakaki, Y.Takizawa, Y.Kawasaki
7 * All rights reserved.
8 * $Id$
9 * -------------------------------------------------------------------------
10 */
11#define NTRACE
12
13#include <stdio.h>
14#include <stdlib.h>
15#include <ctype.h>
16#include <time.h>
17#include <math.h>
18#include <string.h>
19#include "EConst.hh"
20#include "Nrand_gen.hh"
21#include "Ntrace_lenses.hh"
22#include "Ntrace_optics.hh"
23#include "Nphoton.hh"
24
25//#define DEBUG
26
27using namespace NTraceLens;
28
29
30/* -------------------------------------------------------------------------
31 * int trace_main(tel_param *param, NPhoton *photon, FILE *fp)
32 *
33 *   --- trace photons in the EUSO optics (Entrance to FS)
34 *
35 * [input value]
36 * param  : tel_param* to store the optics parameters
37 *            (see also in tracemain_optF1v2.h)
38 * photon : info of photon (position, vector, time, wavelength,...)
39 * fp       : FILE* for (error) message logging
40 *
41 * [return value]
42 * 0: normal end
43 * <0: error (see trace_optF1v2.h for more detail)
44 * -------------------------------------------------------------------------
45 */
46int NTraceLens::trace_main(tel_param *param, NPhoton *photon, FILE *fplog)
47{
48  static double *z0;
49  static int N_FOCAL_SURFACE;
50  static int fs_init=0;
51  int j;
52  int flg, fsloop;
53  double r,r2,dist,xx,yy,zz;
54
55  if(!fs_init)
56    {
57      fs_init=1;
58      N_FOCAL_SURFACE=(int)(param->r_fs/FS_DEF_STEP);
59      z0=(double *)malloc(sizeof(double)*N_FOCAL_SURFACE);
60      if(!z0)
61        {
62          fprintf(fplog,"Cannot allocate memory for array for focalsurface\n");
63          exit(-10);
64        }
65      for(j=0;j<N_FOCAL_SURFACE;j++)
66        {
67          r=(double)j*FS_DEF_STEP;
68          z0[j]=FocalSurface(r,param);
69        }
70    }
71
72  /* the contents of photon will be overwritten in trace_lens() */
73  flg=trace_lens(param,photon,fplog);
74  if(flg!=0)
75    {
76#ifdef DEBUG
77      if(fplog)
78        fprintf(fplog,"flg=%d in trace_main\n",flg);
79#endif
80      return flg;
81    }
82  else if(photon->px[Z]<=0.0)
83    {
84      return (RT_ESCAPE_FROM_ENTRANCE+SURFACE1);
85    }
86  else /* photon->px[Z]>0.0 */
87    {
88      j=0;
89      zz=z0[j];
90      r=1.0e30;
91      fsloop=0;
92      do {
93        fsloop++;
94        r2=r;
95
96        dist=(zz-photon->x[Z])/photon->px[Z];
97        xx=photon->px[X]*dist+photon->x[X];
98        yy=photon->px[Y]*dist+photon->x[Y];
99        r=sqrt(xx*xx+yy*yy);         /* distance to optical axis */
100
101        j=int(r/FS_DEF_STEP);
102        if(j>=N_FOCAL_SURFACE || j < 0) /*changed by francesco 15-01-2010*/
103          j=N_FOCAL_SURFACE-1; /* should be error? */
104
105        if(fabs(z0[j]-zz)<FS_DEF_STEP)
106          {
107            if(r > param->r_fs)
108              flg=1;       /* out of focal surface */
109            else/*changed by francesco 15-01-2010*/
110              zz=FocalSurface(r,param);
111          }
112        else
113          zz=z0[j];
114      }while(fabs(r-r2)>FS_TOLERANCE && !flg && fsloop<FSLOOP_MAX);
115
116      if(flg)
117        {
118#ifdef DEBUG
119          fprintf(stderr,"curpos(%.1f,%.1f,%.1f) vec(%.2f,%.2f,%.2f) "
120                  "est pos on FS(%.1f,%.1f,%.1f) r=%.1f\n",
121                  photon->x[X],photon->x[Y],photon->x[Z],
122                  photon->px[X],photon->px[Y],photon->px[Z],
123                  xx,yy,zz,r
124                  );
125
126#endif
127          return (SURFACE_FS+RT_OUT_OF_SURFACE1);
128        }
129      else if(fsloop>=FSLOOP_MAX)
130        {
131          return (SURFACE_FS+RT_STRAY_LIGHT);
132        }
133      else
134        {
135          photon->x[X]=xx;
136          photon->x[Y]=yy;
137          photon->x[Z]=zz;
138          photon->t+=dist*param->material[param->surface[3].matid].n[0]/EConst::Clight();
139
140//DN      StorePhotonHistory(photon);
141          if(CheckRcut(param,photon,SURFACE_FS))
142            return (SURFACE_FS+RT_OUT_OF_SURFACE1);
143        }
144    }
145
146
147  return 0;
148}
149
150/* -------------------------------------------------------------------------
151 * int read_tel_param(char *filename, tel_param *param, FILE *fplog)
152 *
153 *  --- read optics parameters from a file
154 *
155 *  [input value]
156 *  filename: the name of the file for optics parameters
157 *  param   : pointer to tel_param structure (see also tracemain_optF1v2.h)
158 *  fplog   : FILE* for (error) message logging
159 *
160 *  [return value]
161 *  -1 : "(char *)filename" not found
162 *   0 : normal end
163 * -------------------------------------------------------------------------
164 */
165int NTraceLens::read_tel_param(char *filename, tel_param *param, FILE *fplog)
166{
167  FILE *fr;
168  char buf[256],item[256],tname[256],tfile[256];
169  double val;
170  int i, ival;
171
172  /* initialize tel_param structure */
173  param->flag_printray=0;
174
175  param->nmat=2;
176  param->material[0].name=strdup("VACUUM");
177  param->material[0].ndata=1;
178  param->material[0].n=(double *)malloc(sizeof(double));
179  param->material[0].k=(double *)malloc(sizeof(double));
180  param->material[0].n[0]=1.0;
181  param->material[0].k[0]=0.0;
182
183  param->material[1].name=strdup("AIR");
184  param->material[1].ndata=1;
185  param->material[1].n=(double *)malloc(sizeof(double));
186  param->material[1].k=(double *)malloc(sizeof(double));
187  param->material[1].n[0]=1.0003;
188  param->material[1].k[0]=0.0;
189
190  param->DOE_order=1;
191  param->DOE_efficiency=0.95;
192  for(i=0;i<NSURFACES;i++){
193    param->surface[i].r_cutx=param->surface[i].r_cut=1.e10;
194  }
195  param->r_cut_fs=param->r_cutx_fs=1.e10;
196
197
198  if(( fr=fopen(filename,"r") )==NULL)
199    {
200      if(fplog)
201        fprintf(fplog,"Telescope parameter file %s cannot be opened.\n",filename);
202      return -1;
203    }
204  while(fgets(buf,255,fr))
205    {
206      if(buf[0]!=0 || buf[0]!='#')
207        {
208          sscanf(buf,"%255s",item);
209          switch(item[0])
210            {
211            case 'D':
212              sscanf(buf,"%s",tname);
213              if(strcmp(tname,"DOEord")==0){
214                sscanf(buf,"%*s %d",&ival);
215                param->DOE_order=ival;
216              }else if(strcmp(tname,"DOEeff")==0){
217                sscanf(buf,"%*s %lf",&val);
218                if(val>=0. && val<=1.0)
219                  param->DOE_efficiency=val;
220              }
221              break;
222            case 'F':
223              sscanf(buf,"%*s %lf",&val);
224              switch(item[2])
225                {
226                case '_':
227                  if(item[3]=='C')
228                    param->FS_C=val;
229                  else if(item[3]=='K')
230                    param->FS_K=val;
231                  else if(item[3]=='O')
232                    param->FS_OFFSET=val;
233                  break;
234
235                case 'A':
236                case 'B':
237                case 'C':
238                case 'D':
239                  param->FS[item[2]-'A']=val;
240                  break;
241                }
242              break;
243
244            case 'l':
245              switch(item[1])
246                {
247                case 'e':
248                  sscanf(buf,"%*s %255s",param->lens_dir);
249                  break;
250                case 'a':
251                  sscanf(buf,"%*s %lf",&param->lambda0);
252                  break;
253                }
254              break;
255
256            case 'M':
257              if(isdigit(item[1]) && item[1]>='0' && item[1]<'0'+NSURFACES)
258                {
259                  sscanf(buf,"%*s %256s",tname);
260                  for(i=0; i<param->nmat; i++)
261                    if(!strcmp(tname,param->material[i].name))
262                      {
263                        param->surface[item[1]-'0'].matid=i;
264                        break;
265                      }
266                  if(i==param->nmat){
267                    fprintf(stderr,"%s: Material \"%s\" is NOT registered!\n",
268                            item,tname);
269                    exit(-20);
270                  }
271                }
272              else if(item[1]=='a')
273                {
274                  sscanf(buf,"%*s %256s %256s",tname,tfile);
275                  for(i=0;i<param->nmat;i++)
276                    if(!strcmp(tname,param->material[i].name))
277                      {
278                        fprintf(stderr,"%s is already registered!\n",tname);
279                        break;
280                      }
281                  if(i==param->nmat)
282                    {
283                      param->material[param->nmat].name=strdup(tname);
284                      param->material[param->nmat].filename=strdup(tfile);
285                      param->material[param->nmat].zn=NULL;
286                      param->material[param->nmat].zk=NULL;
287                      param->nmat++;
288                    }
289                }
290              break;
291
292            case 'N':
293              sscanf(buf,"%*s %255s",tname);
294              param->name=strdup(tname);
295              break;
296
297            case 'P':
298              if(!strcmp(item,"PrintRay"))
299                param->flag_printray=1;
300              break;
301
302            case 'R':
303              sscanf(buf,"%*s %lf",&val);
304              switch(item[1])
305                {
306                case '_':
307                  switch(item[2])
308                    {
309                    case 'L':
310              if(strncmp(item,"R_LENS",6)==0)
311            {
312              if(item[6]&& item[6]>='1' && item[6]<='7')
313                param->surface[item[6]-'0'].r_lens=val;
314              for(i=0;i<NSURFACES+1;i++)
315                param->surface[i].r_lens=val;
316            }
317                      break;
318                    case 'W':
319                      param->r_wall=val;
320                      break;
321                    case 'F':
322                      param->r_fs=val;
323                      break;
324                    case 'c':
325                    case 'y':
326                      if(item[5]-'0'<NSURFACES)
327                        param->surface[item[5]-'0'].r_cut=val;
328                      else if(item[5]=='8') /* FS */
329                        param->r_cut_fs=val;
330                      break;
331                    case 'x':
332                      if(item[5]-'0'<NSURFACES)
333                        param->surface[item[5]-'0'].r_cutx=val;
334                      else if(item[5]=='8') /* FS */
335                        param->r_cutx_fs=val;
336                      break;
337                    }
338                  break;
339
340                case 'c':
341                  if(isdigit(item[2]) && item[2]>'0' && item[2] < '0'+NSURFACES)
342                    param->surface[item[2]-'0'].Rc=val;
343                  else if(item[2]=='D') /* for DOE */
344                    param->surface[SURFDOE].Rc=val;
345                  break;
346
347                case 's':
348                  if(isdigit(item[2]) && item[2]>'0' && item[2] < '0'+NSURFACES)
349                    param->surface[item[2]-'0'].r_lens=val;
350                  else if(item[2]=='D')
351                    {
352                      param->surface[SURFDOE].r_lens=val;
353                      sscanf(buf,"%*s %*f %256s",tfile);
354                      param->surface[SURFDOE].filename=strdup(tfile);
355                    }
356                  break;
357                }
358
359            case 'S':
360              sscanf(buf,"%*s %lf",&val);
361              if(item[1]=='z')
362                {
363                  if(isdigit(item[2]) && item[2]<'0'+NSURFACES)
364                    param->surface[item[2]-'0'].Sz=val;
365                  else if(item[2]=='D')
366                    param->surface[SURFDOE].Sz=val;
367                  if(item[2]>'0' && item[2]<'0'+NSURFACES)
368                    {
369                      sscanf(buf,"%*s %*f %s",tfile);
370                      param->surface[item[2]-'0'].filename=strdup(tfile);
371                    }
372                }
373              break;
374
375            default:
376              ; /* skip */
377            }
378        }
379  }
380  fclose(fr);
381
382  for(i=1;i<=7;i++)
383    param->surface[i].Sz += param->surface[0].Sz;
384
385  for(i=1;i<8;i++)
386    {
387      param->surface[i].Zc = param->surface[i].Rc+param->surface[i].Sz;
388    }
389
390  param->FS_OFFSET += param->surface[0].Sz;
391
392  param->surface[0].r_cut=3000.0;
393  param->surface[0].r_cutx=3000.0;
394  param->surface[5].r_cut=param->surface[4].r_cut;
395  param->surface[2].r_cut=param->surface[1].r_cut;
396  param->surface[7].r_cut=param->surface[6].r_cut;
397  param->surface[5].r_cutx=param->surface[4].r_cutx;
398  param->surface[2].r_cutx=param->surface[1].r_cutx;
399  param->surface[7].r_cutx=param->surface[6].r_cutx;
400
401  param->surface[2].r_lens=param->surface[1].r_lens;
402  param->surface[5].r_lens=param->surface[4].r_lens;
403  param->surface[7].r_lens=param->surface[6].r_lens;
404
405  return 0;
406}
407
408/* -------------------------------------------------------------------------
409 * int print_param(tel_param *p)
410 *
411 *  --- print optics parameters for debbuging purpose
412 *
413 *  [input value]
414 *  p: pointer to tel_param structure
415 *
416 *  [return value]
417 *  always 0
418 *
419 * -------------------------------------------------------------------------
420 */
421int NTraceLens::print_param(tel_param *p)
422{
423  int i;
424
425  printf("*** JEM-EUSO Optics <%s> ***\n",p->name);
426  for(i=1;i<=7;i++)
427    printf("r_lens(%d)=%.3f\n",i,p->surface[i].r_lens);
428  printf("r_lens(DOE)=%.3f\n\n",p->surface[SURFDOE].r_lens);
429
430  printf("r_wall=%.3f ",p->r_wall);
431  printf("r_fs=%.3f\n",p->r_fs);
432  for(i=1;i<=7;i++)
433    {
434      printf("Rc%d=%.5f \t",i,p->surface[i].Rc);
435      printf("Sz%d=%.5f \t",i,p->surface[i].Sz);
436      printf("Zc%d=%.5f\n",i,p->surface[i].Zc);
437    }
438  printf("dir=%s\n",p->lens_dir);
439
440  printf("FS_C=%g ",p->FS_C);
441  printf("FS_K=%.2f\n",p->FS_K);
442  for(i=0;i<4;i++)
443    printf("FS%c=%g ",i+'A',p->FS[i]);
444
445  printf("FS_offset=%.2f\n",p->FS_OFFSET);
446  for(i=1;i<NSURFACES;i++)
447    printf("Rcut%d(x)=%.2e   Rcut%d(y)=%.2e    Rlens=%.2f\n",
448           i,p->surface[i].r_cutx,
449           i,p->surface[i].r_cut,
450           p->surface[i].r_lens);
451
452  printf("\nDOE order=%d\n",p->DOE_order);
453  printf("DOE efficiency=%.3f\n",p->DOE_efficiency);
454
455  printf("\n\n");
456  for(i=1;i<7;i++)
457    printf("material(S%d):%s\n",i,p->material[p->surface[i].matid].name);
458
459  return 0;
460}
461
462/* -------------------------------------------------------------------------
463 * double FocalSurface(double r, tel_param *p)
464 *
465 *    --- calculate the Focal surface height as a function of the radius
466 *        from the center
467 *        As photons move toward z>0, FS_C should always be negative.
468 *
469 *  [input value]
470 *  r: radius in mm
471 *  p: pointer to tel_param structure (optics parameter)
472 *
473 *  [return value]
474 *  relative height of FS to that at r=0 in mm
475 *
476 * -------------------------------------------------------------------------
477 */
478/* photons move to the direction z>0 */
479double NTraceLens::FocalSurface(double r /* mm */, tel_param *p)  /* mm */
480{
481  double r2;
482
483  r2=r*r;
484  return r2*(p->FS_C /(1.0+sqrt(1.0-(1.0 + p->FS_K)* p->FS_C * p->FS_C *r2))
485             +r2*(p->FS[0] +r2*(p->FS[1]+r2*(p->FS[2] + r2* p->FS[3]))))+ p->FS_OFFSET;
486}
Note: See TracBrowser for help on using the repository browser.