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