source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/optics/src/Mtracemain_optF1v4.cc @ 117

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

ESAF version compilable on mac OS
File size: 10.1 KB
Line 
1/* -------------------------------------------------------------------------
2 *   tracemain_optF1v6.c
3 *
4 *   ---  main program for raytracing
5 *
6 * Copyright (c) 2000-2007 N.Sakaki, Y.Takizawa, Y.Kawasaki
7 * All rights reserved.
8 * $Id$
9 * -------------------------------------------------------------------------
10 */
11#include <stdio.h>
12#include <stdlib.h>
13//#include <time.h>
14#include <math.h>
15#include <string.h>
16#include "EConst.hh"
17#include "EsafRandom.hh"
18#include "Mtrace_optF1v4.hh"
19#include "Mtracemain_optF1v4.hh"
20#define X      0
21#define Y      1
22#define Z      2
23#define PX     3
24#define PY     4
25#define PZ     5
26#define L      6
27#define T      7
28
29/* -------------------------------------------------------------------------
30 * int trace_main(tel_param *param, double ph_in[8], double ph_out[8],
31 *            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 * ph_in[]: info on injected photon
39 * ph_out[]: info on output photon
40 * ph_xx[0,1,2] : x,y,z position in mm
41 * ph_xx[3,4,5] : x,y,z direction (unit vector)
42 * ph_xx[6] : wavelength in nm
43 * ph_xx[7] : time in ns
44 * fp       : FILE* for (error) message logging
45 *
46 * [return value]
47 * 0: normal end
48 * <0: error (see trace_optF1v2.h for more detail)
49 * -------------------------------------------------------------------------
50 */
51int Mtrace_main(Mtel_param *param,
52           double ph_in[8], double ph_out[8], FILE *fplog)
53{
54  static double *z0;
55  static int N_FOCAL_SURFACE;
56  static int fs_init=0;
57  int j;
58  int flg, fsloop;
59  double r,r2,dist,xx,yy,zz;
60
61  if(!fs_init)
62    {
63      fs_init=1;
64      N_FOCAL_SURFACE=(int)(param->Mr_fs/MFS_DEF_STEP);
65      z0=(double *)malloc(sizeof(double)*N_FOCAL_SURFACE);
66      if(!z0)
67    {
68      fprintf(fplog,"Cannot allocate memory for array for focalsurface\n");
69      exit(-10);
70    }
71      for(j=0;j<N_FOCAL_SURFACE;j++)
72    {
73      r=(double)j*MFS_DEF_STEP;
74      z0[j]=MFocalSurface(r,param);
75    }
76    }
77
78  /* the contents of ph_in[] will be overwritten in trace_lens() */
79  flg=trace_lens(param,ph_in,ph_out,fplog);
80  if(flg!=0)
81    {
82#ifdef DEBUG
83      if(fplog)
84    fprintf(fplog,"flg=%d in trace_main\n",flg);
85#endif
86      return flg;
87    }
88  else if(ph_out[PZ]<=0.0)
89    {
90      return (RT_ESCAPE_FROM_ENTRANCE+SURFACE1);
91    }
92  else /* ph_out[PZ]>0.0 */
93    {
94      j=0;
95      zz=z0[j];
96      r=1.0e30;
97      fsloop=0;
98      do {
99    fsloop++;
100    r2=r;
101
102    dist=(zz-ph_out[Z])/ph_out[PZ];
103    xx=ph_out[PX]*dist+ph_out[X];
104    yy=ph_out[PY]*dist+ph_out[Y];
105    r=sqrt(xx*xx+yy*yy);         /* distance to optical axis */
106
107    j=int(r/MFS_DEF_STEP);
108
109    if(j>=N_FOCAL_SURFACE || j < 0)/*changed by francesco 15-01-2010*/
110      j=N_FOCAL_SURFACE-1; /* should be error? */
111
112    if(fabs(z0[j]-zz)<MFS_DEF_STEP)
113      {
114        if(r > param->Mr_fs)
115          flg=1;       /* out of focal surface */
116        else/*changed by francesco 15-01-2010*/
117        zz=MFocalSurface(r,param);
118      }
119    else
120      zz=z0[j];
121      }while(fabs(r-r2)>MFS_TOLERANCE && !flg && fsloop<MFSLOOP_MAX);
122
123      if(flg)
124    {
125#ifdef DEBUG
126      fprintf(stderr,"curpos(%.1f,%.1f,%.1f) vec(%.2f,%.2f,%.2f) "
127          "est pos on FS(%.1f,%.1f,%.1f) r=%.1f\n",
128          ph_out[X],ph_out[Y],ph_out[Z],
129          ph_out[PX],ph_out[PY],ph_out[PZ],
130          xx,yy,zz,r
131          );
132#endif
133      return (SURFACE_FS+RT_OUT_OF_SURFACE1);
134    }
135      else if(fsloop>=MFSLOOP_MAX)
136    {
137      return (SURFACE_FS+RT_STRAY_LIGHT);
138    }
139      else
140    {
141      ph_out[X]=xx;
142      ph_out[Y]=yy;
143      ph_out[Z]=zz;
144      ph_out[T]+=dist/EConst::Clight();
145          if(CheckRcut(param,ph_out,SURFACE_FS))
146            return (SURFACE_FS+RT_OUT_OF_SURFACE1);
147    }
148    }
149
150  return 0;
151}
152
153/* -------------------------------------------------------------------------
154 * int read_tel_param(char *filename, tel_param *param, FILE *fplog)
155 *
156 *  --- read optics parameters from a file
157 *
158 *  [input value]
159 *  filename: the name of the file for optics parameters
160 *  param   : pointer to tel_param structure (see also tracemain_optF1v2.h)
161 *  fplog   : FILE* for (error) message logging
162 *
163 *  [return value]
164 *  -1 : "(char *)filename" not found
165 *   0 : normal end
166 * -------------------------------------------------------------------------
167 */
168int Mread_tel_param(const char *filename, Mtel_param *param, FILE *fplog)
169{
170  FILE *fr;
171  char buf[256],item[256];
172  double val;
173
174  if(( fr=fopen(filename,"r") )==NULL)
175    {
176      if(fplog)
177    fprintf(fplog,"Telescope parameter file %s cannot be opened.\n",filename);
178      return -1;
179    }
180  while(fgets(buf,255,fr))
181    {
182      if(buf[0]!=0 || buf[0]!='#')
183    {
184      sscanf(buf,"%255s",item);
185      switch(item[0])
186        {
187        case 'R':
188          sscanf(buf,"%*s %lf",&val);
189          switch(item[1])
190        {
191        case '_':
192          switch(item[2])
193            {
194            case 'L':
195              param->Mr_lens=val;
196              break;
197            case 'W':
198              param->Mr_wall=val;
199              break;
200            case 'F':
201              param->Mr_fs=val;
202              break;
203                    case 'c':
204                      param->Mr_cut[item[5]-'0']=val;
205                      break;
206            }
207          break;
208        case 'c':
209          switch(item[2])
210            {
211            case '1':
212              param->MRc1=val;
213              break;
214            case '2':
215              param->MRc2=val;
216              break;
217            case '3':
218              param->MRc3=val;
219              break;
220            case '6':
221              param->MRc6=val;
222              break;
223            case '7':
224              param->MRc7=val;
225              break;
226            case 'D':
227              param->MRcD=val;
228              break;
229            }
230          break;
231        }
232          break;
233
234        case 'S':
235          sscanf(buf,"%*s %lf",&val);
236          if(item[1]=='z')
237        {
238          switch(item[2])
239            {
240            case '0':
241              param->MSz0=val;
242              break;
243            case '1':
244              param->MSz1=val;
245              break;
246            case '2':
247              param->MSz2=val;
248              break;
249            case '3':
250              param->MZc3=val;
251              break;
252            case '4':
253              param->MSz4=val;
254              break;
255            case '6':
256              param->MSz6=val;
257              break;
258            case '7':
259              param->MSz7=val;
260              break;
261            case 'D':
262              param->MZcD=val;
263              break;
264            }
265        }
266          break;
267
268        case 'F':
269          sscanf(buf,"%*s %lf",&val);
270          switch(item[2])
271        {
272        case '_':
273          if(item[3]=='C')
274            param->MFS_C=val;
275          else if(item[3]=='K')
276            param->MFS_K=val;
277          else if(item[3]=='O')
278            param->MFS_OFFSET=val;
279          break;
280
281        case 'A':
282          param->MFSA=val;
283          break;
284        case 'B':
285          param->MFSB=val;
286          break;
287        case 'C':
288          param->MFSC=val;
289          break;
290        case 'D':
291          param->MFSD=val;
292          break;
293        }
294          break;
295
296        case 'l':
297          switch(item[1])
298        {
299        case 'e':
300          sscanf(buf,"%*s %255s",param->Mlens_dir);
301          break;
302        case 'a':
303          sscanf(buf,"%*s %lf",&param->Mlambda0);
304          break;
305        }
306          break;
307        }
308    }
309    }
310  fclose(fr);
311
312  param->MSz1 += param->MSz0;
313  param->MSz2 += param->MSz0;
314  param->MZc3 += param->MSz0;
315  param->MSz4 += param->MSz0;
316  param->MSz6 += param->MSz0;
317  param->MSz7 += param->MSz0;
318  param->MZcD += param->MSz0;
319
320  param->MZc1 =  param->MRc1+param->MSz1;
321  param->MZc2 =  param->MRc2+param->MSz2;
322  param->MZc4 =             param->MSz4;
323  param->MZc6 = -param->MRc6+param->MSz6;
324  param->MZc7 = -param->MRc7+param->MSz7;
325
326  param->MFS_OFFSET += param->MSz0;
327
328  param->Mr_cut[0]=3000.0;
329  param->Mr_cut[5]=param->Mr_cut[4];
330  param->Mr_cut[2]=param->Mr_cut[1];
331  param->Mr_cut[7]=param->Mr_cut[6];
332
333  return 0;
334}
335
336/* -------------------------------------------------------------------------
337 * int print_param(tel_param *p)
338 *
339 *  --- print optics parameters for debbuging purpose
340 *
341 *  [input value]
342 *  p: pointer to tel_param structure
343 *
344 *  [return value]
345 *  always 0
346 *
347 * -------------------------------------------------------------------------
348 */
349int Mprint_param(Mtel_param *p)
350{
351  int i;
352
353  printf("Mr_lens=%.3f ",p->Mr_lens);
354  printf("Mr_wall=%.3f ",p->Mr_wall);
355  printf("Mr_fs=%.3f\n",p->Mr_fs);
356  printf("MRc1=%.5f ",p->MRc1);
357  printf("MSz1=%.5f ",p->MSz1);
358  printf("MZc1=%.5f\n",p->MZc1);
359
360  printf("MRc2=%.5f ",p->MRc2);
361  printf("MSz2=%.5f ",p->MSz2);
362  printf("MZc2=%.5f\n",p->MZc2);
363
364  printf("MSz4=%.5f ",p->MSz4);
365
366  printf("MRc6=%.5f ",p->MRc6);
367  printf("MSz6=%.5f ",p->MSz6);
368  printf("MZc6=%.5f\n",p->MZc6);
369
370  printf("MRc7=%.5f ",p->MRc7);
371  printf("MSz7=%.5f ",p->MSz7);
372  printf("MZc7=%.5f\n",p->MZc7);
373
374  printf("Rc3=%.5f ",p->MRc3);
375  printf("MSz1=%.5f ",p->MSz1);
376  printf("MZc3=%.5f\n",p->MZc3);
377
378  printf("dir=%s\n",p->Mlens_dir);
379
380  printf("FS_C=%g ",p->MFS_C);
381  printf("FS_K=%.2f\n",p->MFS_K);
382  printf("FSA=%g ",p->MFSA);
383  printf("FSB=%g ",p->MFSB);
384  printf("FSC=%g ",p->MFSC);
385  printf("FSD=%g ",p->MFSD);
386  printf("FS_offset=%.2f\n",p->MFS_OFFSET);
387  for(i=1;i<9;i++)
388    printf("Rcut%d=%.2f\n",i,p->Mr_cut[i]);
389  return 0;
390}
391
392/* -------------------------------------------------------------------------
393 * double FocalSurface(double r, tel_param *p)
394 *
395 *    --- calculate the Focal surface height as a function of the radius
396 *        from the center
397 *        As photons move toward z>0, FS_C should always be negative.
398 *
399 *  [input value]
400 *  r: radius in mm
401 *  p: pointer to tel_param structure (optics parameter)
402 *
403 *  [return value]
404 *  relative height of FS to that at r=0 in mm
405 *
406 * -------------------------------------------------------------------------
407 */
408/* photons move to the direction z>0 */
409double MFocalSurface(double r /* mm */, Mtel_param *p)  /* mm */
410{
411  double r2;
412
413  r2=r*r;
414  return r2*(p->MFS_C /(1.0+sqrt(1.0-(1.0 + p->MFS_K)* p->MFS_C * p->MFS_C *r2))
415         +r2*(p->MFSA +r2*(p->MFSB+r2*(p->MFSC + r2* p->MFSD))))+ p->MFS_OFFSET;
416}
Note: See TracBrowser for help on using the repository browser.