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

/* -------------------------------------------------------------------------
 *   trace_optF1v1.c
 *
 *   --- main program of raytracing
 *
 * $Id: Ktrace_optF1v1.cc 2924 2011-06-12 20:22:13Z mabl $
 * -------------------------------------------------------------------------
 */
#include <stdio.h>
#include <stdlib.h>
#include "TMath.h"
#include "EsafRandom.hh"
#include "Kspline.hh"
#include "KEUSO_optF1v1.hh"
#include "Kutils.hh"
#include "Ktrace_optF1v1.hh"

#define RT_SEARCH_ITER_MAX        100
#define RT_INTERACTION_MAX         20

#define X      0
#define Y      1
#define Z      2
#define PX     3
#define PY     4
#define PZ     5
#define L      6
#define T      7

/* lens surface data */
static Int_t n_s1=0,n_s2=0,n_s5=0,n_s6=0;
static Double_t sr1[S_DATA],sz1[S_DATA];
static Double_t sr2[S_DATA],sz2[S_DATA];
static Double_t sr5[S_DATA],sz5[S_DATA];
static Double_t sr6[S_DATA],sz6[S_DATA];

/* lens material data */
static Int_t n_acryl=0, n_bg3=0;
static Double_t ac_la[Acryl_N],ac_n[Acryl_N],ac_k[Acryl_N];
static Double_t bg3_la[Bg3_N],bg3_n[Bg3_N],bg3_k[Bg3_N];

/* diffractive structure data */
static Int_t n_sd1=0, n_sd6=0;
#define N_DIFFR_DATA   13000
static Double_t SD1_R[N_DIFFR_DATA];
static Double_t SD1_D[N_DIFFR_DATA];
static Double_t SD6_R[N_DIFFR_DATA];
static Double_t SD6_D[N_DIFFR_DATA];

static Double_t Kref(Double_t n1, Double_t n2, Double_t ci, Double_t ct);
static Int_t Kcheck_intersect_cylinder(Double_t r1, Double_t z1, Double_t r2, Double_t z2,
                    Double_t phot[8], Double_t *dis);
static Int_t Kcheck_intersect_cone    (Double_t r1, Double_t z1, Double_t r2, Double_t z2,
                    Double_t phot[8], Double_t *dis);
static Int_t KCheckIris(Double_t phot[6]);

/*
 * ========================================================================
 * functions for diffractive optics start
 * ========================================================================
 */
/*
 * ------------------------------------------------------------------------
 *  Calculate normal vector to the sphere approximating the lens surface
 *  Input:
 *    aX, aY, aZ: photon position on the lens surface
 *    aCZ       : center of the approximated sphere
 *
 *  Output:
 *    nnnDIFF[3]: normal unit vector
 * ------------------------------------------------------------------------
 */
Int_t KCalSpheVec(Double_t aX, Double_t aY, Double_t aZ, Double_t aCZ, Double_t nnnDIFF[3])
{
        
  Double_t N_;
        
  nnnDIFF[0] = aX;
  nnnDIFF[1] = aY;
  nnnDIFF[2] = aZ - aCZ;
        
  N_ = TMath::Sqrt(Kvdot(nnnDIFF,nnnDIFF));
        
  nnnDIFF[0] = nnnDIFF[0]/N_;
  nnnDIFF[1] = nnnDIFF[1]/N_;
  nnnDIFF[2] = nnnDIFF[2]/N_;

  return 0;
}


/*
 * ------------------------------------------------------------------------
 *  read parameters of diffractive optics
 *    Input:
 *        dirname    the name of the directory name
 *                        where data files are stored
 *    Output:
 *        (return value) status
 * ------------------------------------------------------------------------
 */
Int_t KReadSD(const char *dirname)
{
  FILE *fr;
  char Infname[128];
        
#ifdef DEBUG
  fprintf(stderr,"Reading Diffractive Surface Data....  ");
#endif

  /* read diffractive optics data for S1 */
  sprintf(Infname,"%.100s/SD1_F1v1.dat",dirname);
        
  if(( fr = fopen(Infname,"r" )) == NULL)
    {
      fprintf(stderr,"%s not OPEN\n",Infname );
      return (-1);
    }
  n_sd1=0;
  while(n_sd1<N_DIFFR_DATA 
    && fscanf(fr,"%lf %lf",&SD1_R[n_sd1],&SD1_D[n_sd1])==2)
    n_sd1++;
  fclose(fr);
  if(n_sd1==N_DIFFR_DATA)
    return (-2);

  /* read diffractive optics data for S6 */
  sprintf(Infname,"%.100s/SD6_F1v1.dat",dirname);
        
  if(( fr = fopen(Infname,"r" )) == NULL)
    {
      fprintf(stderr,"%s not OPEN\n",Infname);
      return (-1);
    }

  n_sd6=0;
  while(n_sd6<N_DIFFR_DATA && 
    fscanf(fr,"%lf %lf",&SD6_R[n_sd6],&SD6_D[n_sd6])==2)
    n_sd6++;
  fclose(fr);

  if(n_sd6==N_DIFFR_DATA)
    return (-2);
        
#ifdef DEBUG
  fprintf(stderr,"Done.\n");
#endif
  return (0);
}

/*
 * ------------------------------------------------------------------------
 *  calculate the diffractive structure
 *   Input:
 *      r   the distance in mm to the optical axis at the intersection point
 *      s   surface ID (i.e. S1,S2,...)
 *   Output:
 *     (return value)
 * ------------------------------------------------------------------------
 */
Double_t KCal_D(Double_t r, Int_t s) 
{
  Int_t i;

  if(s==1)
    {
      i=Krt_bsearch(SD1_R,r,0,n_sd1-1);
      return (SD1_D[i]);
    }

  if(s==6)
    {
      i=Krt_bsearch(SD6_R,r,0,n_sd6-1);
      return (SD6_D[i]);
    }

  return 0.0;
}

/*
 * ------------------------------------------------------------------------
 *  calculate the diffracted ray
 *  Input:
 *        in_ph[3]      array x,y,z of photon
 *        nnn[3]        normal vector to the diffractive surface
 *                         at the intersection
 *        n1            refractive index of xxx
 *        n2            refractive index of xxx
 *        d
 *        wl            wavelength in mm
 *        m             index
 *  Output:
 *        in_ph[3]
 * ------------------------------------------------------------------------
 */
Int_t Kdiffract(Int_t ID, Double_t in_phdir[3], 
         Double_t x3, Double_t y3, Double_t z3, Double_t Zc,
         Double_t n1, Double_t n2, Double_t wl, Int_t m)
{
  Int_t    i;
  Double_t CosIN, CosD, dir, nnnR[3], r3, ppp[3], nnn[3], d, norm;

  r3=TMath::Sqrt(x3*x3+y3*y3);
  nnnR[0]= x3/r3;
  nnnR[1]= y3/r3;
  nnnR[2]=0.0;

  KCalSpheVec(x3, y3, z3, Zc, nnn);
  d=KCal_D(r3,ID);

  if(in_phdir[Z]>0.)
    dir=1.0;
  else
    dir=-1.0;

  CosIN=Kvdot(in_phdir,nnn);
  CosD=Kvdot(nnn,nnnR);

  for(i=0;i<3;i++)
    ppp[i]=(n1*(in_phdir[i]-CosIN*nnn[i])
      +(Double_t)m*wl/357.0e-6*d*(nnnR[i]-CosD*nnn[i])*dir)/n2;

  norm=Kvdot(ppp,ppp);
  if(norm>1.0)
    return -1;

  norm=TMath::Sqrt(1.0-norm);

#ifdef DEBUG
  fprintf(stderr,"Diffractive: (%.3f , %.3f , %.3f)=>",
      in_phdir[0],in_phdir[1],in_phdir[2]);
#endif

  for(i=0;i<3;i++)
    in_phdir[i]=ppp[i]+norm*nnn[i]*(CosIN>0.0?(1.0):(-1.0));

#ifdef DEBUG
  fprintf(stderr,"(%.3f , %.3f , %.3f)\n",in_phdir[0],in_phdir[1],in_phdir[2]);
#endif

  return (0);
}
        
/*
 * ========================================================================
 * functions for diffractive optics END
 * ========================================================================
 */

/*
 * ------------------------------------------------------------------------
 * calculate reflection probability assuming no polariation
 *   n1,n2: refractive index of the media of this side and the other side
 *   ci,ct: cosine of the incident light and reflected light to the surface
 *          normal vector
 * ------------------------------------------------------------------------
 */
static Double_t Kref(Double_t n1, Double_t n2, Double_t ci, Double_t ct)
{
  Double_t rp,rs;

  rp=(n2*ci-n1*ct)/(n2*ci+n1*ct);
  rs=(n1*ci-n2*ct)/(n1*ci+n2*ct);

  return((rp*rp+rs*rs)/2.0);
}

/*
 * ------------------------------------------------------------------------
 *  read parameters of lens surface shapes, and optical indices
 * ------------------------------------------------------------------------
 */
Int_t Kread_para(const char *dir)
{
  FILE *fp;
  char filename[128];

#ifdef DEBUG
  fprintf(stderr,"Reading Surface Data....  ");
#endif

  /* Read S1 surface */
  sprintf(filename,"%.100s/S1_F1v1.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -1;
  n_s1=0;
  while(n_s1<S_DATA && 
        fscanf(fp,"%lf %*f %*f %lf",&sr1[n_s1], &sz1[n_s1])!=EOF)
    {
      sz1[n_s1]=Sz1-sz1[n_s1];
      n_s1++;
    }
  fclose(fp);

  /* Read S2 surface */
  sprintf(filename,"%.100s/S2_F1v1.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -2;
  n_s2=0;
  while(n_s2<S_DATA && 
        fscanf(fp,"%lf %*f %*f %lf",&sr2[n_s2], &sz2[n_s2])!=EOF)
    {
      sz2[n_s2]=Sz2-sz2[n_s2];
      n_s2++;
    }
  fclose(fp);

  /* Read S5 surface */
  sprintf(filename,"%.100s/S5_F1v1.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -5;
  n_s5=0;
  while(n_s5<S_DATA && 
        fscanf(fp,"%lf %*f %*f %lf",&sr5[n_s5], &sz5[n_s5])!=EOF)
    {
      sz5[n_s5]=Sz5-sz5[n_s5];
      n_s5++;
    }
  fclose(fp);

  /* Read S6 surface */
  sprintf(filename,"%.100s/S6_F1v1.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -6;
  n_s6=0;
  while(n_s6<S_DATA &&
        fscanf(fp,"%lf %*f %*f %lf",&sr6[n_s6], &sz6[n_s6])!=EOF)
    {
      sz6[n_s6]=Sz6-sz6[n_s6];
      n_s6++;
    }
  fclose(fp);
#ifdef DEBUG
  fprintf(stderr,"Done.\n");

  fprintf(stderr,"Reading Optical Data....   ");
#endif

  /* Read Acryl data */
  sprintf(filename,"%.100s/Mitsubishi_000.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -10;
  n_acryl=0;
  while(n_acryl<Acryl_N &&
    fscanf(fp,"%lf %lf %lf",&ac_la[n_acryl],
           &ac_n[n_acryl], &ac_k[n_acryl])!=EOF)
    {
      n_acryl++;
    }
  fclose(fp);

  /* Read BG3 data */
  sprintf(filename,"%.100s/bg3.dat",dir);
  if((fp=fopen(filename,"r"))==NULL)
    return -11;
  n_bg3=0;
  while(n_bg3<Bg3_N &&
    fscanf(fp,"%lf %lf %lf",&bg3_la[n_bg3],
           &bg3_n[n_bg3], &bg3_k[n_bg3])!=EOF)
    {
      n_bg3++;
    }
  fclose(fp);

#ifdef DEBUG
  fprintf(stderr,"Done.\n");
#endif
  return 0;
}

/*
 * ------------------------------------------------------------------------
 *  get refractive index of the acryl at wavelength lambda
 *    using cubic spline interpolation
 * ------------------------------------------------------------------------
 */
Double_t KGet_n_Acryl(Double_t lambda)
{
  static Int_t init=0;
  static Double_t *z;

  if(!init)
    {
      init=1;
      if((z=(Double_t *)malloc(sizeof(Double_t)*n_acryl))==NULL)
    return -1000.0;
      Kcsp_maketable(ac_la, ac_n, z, n_acryl);
    }

  return Kcspline(lambda, ac_la, ac_n, z, n_acryl);
}

/*
 * ------------------------------------------------------------------------
 *  get extinction coefficient of the acryl at wavelength lambda
 *    using cubic spline interpolation
 * ------------------------------------------------------------------------
 */
Double_t KGet_k_Acryl(Double_t lambda)
{
  static Int_t init=0;
  static Double_t *z;

  if(!init)
    {
      init=1;
      if((z=(Double_t *)malloc(sizeof(Double_t)*n_acryl))==NULL)
    return -1000.0;
      Kcsp_maketable(ac_la, ac_k, z, n_acryl);
    }

  return Kcspline(lambda, ac_la, ac_k, z, n_acryl);
}

/*
 * ------------------------------------------------------------------------
 *  get the refracted ray
 *    if reflected it is thrown away.
 *
 *   in_ph[3]: incident photon vector (normalized)
 *   nnn[3]  : surface normal vector (direct to this side)
 *   n1,n2   : refractive indices of the media of this side and the other side
 * ------------------------------------------------------------------------
 */
Int_t Krefract(Double_t in_phdir[3],Double_t nnn[3],Double_t n1,Double_t n2)
{
  Int_t i, TotalReflection=0;
  Double_t CosIN, CosOUT=1.0, ppp[3], norm, out_phdir[3];

  CosIN=Kvdot(in_phdir,nnn);
  for(i=0;i<3;i++)
    ppp[i]=(in_phdir[i]-CosIN*nnn[i])*n1/n2;

  norm=Kvdot(ppp,ppp);
#ifdef NO_REFLECTION
  if(norm>1.0)
    return(RT_NO_REFRACTED_LIGHT);  /* total reflection */
#else
  if(norm>1.0)
    TotalReflection=1;
  else
    {
#endif

      norm=TMath::Sqrt(1.0-norm);
      for(i=0;i<3;i++)
    out_phdir[i]=ppp[i]+norm*nnn[i]*(CosIN>0.0?(1.0):(-1.0));

#ifndef NO_REFLECTION
      /* check reflection */
      CosOUT=Kvdot(out_phdir,nnn);
    }
  TRandom *rndm = EsafRandom::Get();
  if(TotalReflection==1 || rndm->Rndm() < Kref(n1,n2,CosIN,CosOUT))
    { 
      for(i=0;i<3;i++)
    in_phdir[i]-=2.0*CosIN*nnn[i];
      return(RT_REFLECTED);
    }
#endif

  /* in_phdir[]: direction of incident ray -> refracted ray */
  for(i=0;i<3;i++)
    in_phdir[i]=out_phdir[i];

  return(RT_REFRACTED);
}

/*
 *
 * check if there are any intersections with a cylinder
 *
 */
static Int_t Kcheck_intersect_cylinder(Double_t r1, Double_t z1, Double_t r2, Double_t z2,
                    Double_t phot[8], Double_t *dis)
{
  Double_t aa[3], xx[2], tmp, z;
  Int_t    nroot, i;

  aa[0] = phot[PX]*phot[PX] + phot[PY]*phot[PY];
  aa[1] = 2.0*(phot[X]*phot[PX] + phot[Y]*phot[PY]);
  aa[2] = (phot[X]*phot[X] + phot[Y]*phot[Y]) - r1*r1;

  nroot=KSolveQuadratic(aa,xx);
  if(nroot<2)
    return(RT_NO_INTERSECTION1);

  /* solution is the smaller positive value */
  /* swap */
  if(xx[0]>xx[1])
    {
      tmp=xx[0];
      xx[0]=xx[1];
      xx[1]=tmp;
    }
  if(xx[1]<0)
    return(RT_NO_INTERSECTION1);

  for(i=0;i<2;i++)
    {
      if(xx[i]>0.0)
    {
      z=phot[Z]+phot[PZ]*xx[i];
      /* check if the solution is in the specified range */
      if((z-z1)*(z-z2)<0)
        {
          *dis=xx[i];
          return 0;
        }
    }
    }

  *dis=(xx[0]>0 ? xx[0]:xx[1]);
  return(RT_NO_INTERSECTION2);
}

/*
 *
 * check if there are any intersections with a cone
 *
 */
static Int_t Kcheck_intersect_cone(Double_t r1, Double_t z1, Double_t r2, Double_t z2,
                Double_t phot[8], Double_t *dis)
{
  Double_t a00, a0, aa[3], xx[2], z0, dbuf1, tmp, x, y, r, zz[2];
  Int_t    nroot, i;

  a00 = (z2-z1)/(r2-r1);
  z0  = z1-a00*r1;
  a0  = a00*a00;

  dbuf1 = phot[Z]-z0;
  aa[0] = a0*phot[PX]*phot[PX] + a0*phot[PY]*phot[PY] - phot[PZ]*phot[PZ];
  aa[1] = 2.0*(a0*phot[PX]*phot[X] + a0*phot[PY]*phot[Y] - phot[PZ]*dbuf1);
  aa[2] = a0*phot[X]*phot[X] + a0*phot[Y]*phot[Y] - dbuf1*dbuf1;

  nroot=KSolveQuadratic(aa,xx);
  if(nroot<2)
    return(RT_NO_INTERSECTION1);

  /* solution is the smaller positive value */
  /* swap */
  if(xx[0]>xx[1])
    {
      tmp=xx[0];
      xx[0]=xx[1];
      xx[1]=tmp;
    }
  if(xx[1]<=0.0)
    return(RT_NO_INTERSECTION1);

  /* We usually have a solution on a virtual cone with -a00 */
  /* check this condition */
  for(i=0;i<2;i++)
    if(xx[i]>0.0)
      {
    x=phot[PX]*xx[i]+phot[X];
    y=phot[PY]*xx[i]+phot[Y];
    zz[i]=phot[PZ]*xx[i]+phot[Z];
    r=TMath::Sqrt(x*x+y*y);
    if((r-r1)*(r-r2)<0.0)
      {
        *dis=xx[i];
        return 0;
      }
      }

  if(xx[0]<=0.0)
    *dis=xx[1];
  else
    {
      /* the nearer one to the surface will be what we want... */
      if(fabs(zz[0]-z1)<fabs(zz[1]-z1))
    *dis=xx[0];
      else
    *dis=xx[1];
    }

  return(RT_NO_INTERSECTION2);
}

/*
 * ------------------------------------------------------------------------
 *
 *  ID:            Lens Surface ID
 *  phot[8]: input and output photon data
 *    phot[0,1,2]: x,y,z position in mm
 *    phot[3,4,5]: l,m,n direction
 *    phot[6]    : wavelength in nm
 *    phot[7]    : time in ns
 *  offset_z     : the offset z-position of sphare center approximating
 *                 lens surface
 *  radius       : the radius of sphare approximated lens surface
 *  lens_r[],lens_z[]: lens surface data sets of the radius and the height
 *  n_lens_elem  : number of data elements of lens_r[]
 *  medium1_n, medium2_n: refractive indices of the media in this side
 *                         and the other side
 * ------------------------------------------------------------------------
 */
Int_t Klens_interaction(Int_t ID, Double_t phot[8],
             Double_t Zc, Double_t radius,
             Double_t *lens_r, Double_t *lens_z, Int_t n_lens_elem,
             Double_t medium1_n, Double_t medium2_n, Double_t *dis)
{
  Int_t r_idx,iter,flag,Found,nroot,dir,status,idx_dir;
  Double_t dbuf1, aa[3], xx[2], x3,y3,z3,r3, r1,r2, z1,z2;
  Double_t nnn[3];
  Double_t a00;
  Double_t tmp, tmp_dis, v[3], tmp_x[3];

  /* as an initial guess, intersection with an approximated sphere
     is calculated */

  dbuf1 = phot[Z]-Zc;
  aa[0] = phot[PX]*phot[PX] + phot[PY]*phot[PY] + phot[PZ]*phot[PZ];
  aa[1] = 2.0*(phot[PX]*phot[X] + phot[PY]*phot[Y] + phot[PZ]*dbuf1);
  aa[2] = phot[X]*phot[X] + phot[Y]*phot[Y] + dbuf1*dbuf1-radius*radius;

  nroot=KSolveQuadratic(aa,xx);
  if(nroot<2)
    return(RT_NO_INTERSECTION1);

  /* swap */
  if(xx[0]>xx[1])
    {
      tmp=xx[0];
      xx[0]=xx[1];
      xx[1]=tmp;
    }

  /* We need a positive, smaller solution. */
  if(xx[0]>0.0)
    dbuf1=xx[0];
  else
    {
      if(xx[1]>0.0)
    dbuf1=xx[1];
      else
    return(RT_NO_INTERSECTION1);
    }

  x3=phot[PX]*dbuf1+phot[X];
  y3=phot[PY]*dbuf1+phot[Y];
  z3=phot[PZ]*dbuf1+phot[Z];
   /* distance to the optical axis at the initially guessed intersection */
  r3=TMath::Sqrt(x3*x3+y3*y3);

  if(r3>R_LENS)
    {
      /* The other solution may be acceptable just for an initial guess. */
      dbuf1=xx[0];
      x3=phot[PX]*dbuf1+phot[X];
      y3=phot[PY]*dbuf1+phot[Y];
      z3=phot[PZ]*dbuf1+phot[Z];

      r3=TMath::Sqrt(x3*x3+y3*y3);
      if(r3>R_LENS)
    return(RT_HIT_SIDE_WALL);
    }

  /* search the index near r3 */
  r_idx=Krt_bsearch(lens_r,r3,0,n_lens_elem-1);

  if(r_idx>n_lens_elem-2)
    r_idx=n_lens_elem-2;

  /* search intersection point with real fresnel lens */
  iter=0; /* number of iterations */
  do{
    Found=0;

    r1=lens_r[r_idx];
    z1=lens_z[r_idx];
    r2=lens_r[r_idx+1];
    z2=lens_z[r_idx+1];

    if(r1==r2)
      {
        /* intersection with a part of the cylinder */
    status=Kcheck_intersect_cylinder(r1, z1, r2, z2, phot, dis);

    x3=phot[PX]*(*dis)+phot[X];
    y3=phot[PY]*(*dis)+phot[Y];
    z3=phot[PZ]*(*dis)+phot[Z];

    if(!status)     /* intersection found */
      {
        Found=1;
        /* normal vector */
        if(z2>z1)
          dir=1;
        else
          dir=-1;
        nnn[X]=dir*x3/r1;
        nnn[Y]=dir*y3/r1;
        nnn[Z]=0.0;
      }
    else /* intersection with cylinder not found */
      {
        if(status==RT_NO_INTERSECTION2 && fabs(z3-z1)<100.0)
          {
        v[X]=x3;
        v[Y]=y3;
        v[Z]=0.0;

        if((z3-z1)*Kvdot(v,&phot[PX])*phot[PZ]>0)
          r_idx--;
        else
          r_idx++;
          }
        else
          {
        *dis=(z1-phot[Z])/phot[PZ];
        z3=z1;
        x3=phot[PX]*(*dis)+phot[X];
        y3=phot[PY]*(*dis)+phot[Y];
        r3=TMath::Sqrt(x3*x3+y3*y3);
        r_idx=Krt_bsearch(lens_r,r3,0,n_lens_elem-1);
          }
      }
#ifdef DEBUG
    fprintf(stderr,"ridx0(%d): %d (%.1f: %.3f %.3f %.3f (%.1f %.1f %.1f)(%.3f %.3f %.3f))\n",
    status,r_idx,z1,r1,r2,r3,phot[X],phot[Y],phot[Z],phot[PX],phot[PY],phot[PZ]);
#endif
      }
    else  /* r1 != r2 */
      {
    /* intersection with a part of the cone */
    status=Kcheck_intersect_cone(r1, z1, r2, z2, phot, dis);

    x3=phot[PX]*(*dis)+phot[X];
    y3=phot[PY]*(*dis)+phot[Y];
    z3=phot[PZ]*(*dis)+phot[Z];

    r3=TMath::Sqrt(x3*x3+y3*y3);

    if(!status)   /* intersection found */
      {
        Found=1;
        /* normal vector */
        a00=(z2-z1)/(r2-r1);
        nnn[X]=a00*x3;
        nnn[Y]=a00*y3;
        nnn[Z]=-r3;
        dbuf1=TMath::Sqrt(Kvdot(nnn,nnn));
        nnn[X]/=dbuf1;
        nnn[Y]/=dbuf1;
        nnn[Z]/=dbuf1;
      }
    else /* not found */
      {
        if(status==RT_NO_INTERSECTION2 &&
           fabs(z3-z1)<10. && fabs(r3-r1)<5.)
          {
        if(r3<r1)
          r_idx--;
        if(r3>r2)
          r_idx++;
          }
        else
          {
        *dis=(z1-phot[Z])/phot[PZ];
        z3=z1;
        x3=phot[PX]*(*dis)+phot[X];
        y3=phot[PY]*(*dis)+phot[Y];
        r3=TMath::Sqrt(x3*x3+y3*y3);
        r_idx=Krt_bsearch(lens_r,r3,0,n_lens_elem-1);
          }
      }
#ifdef DEBUG
fprintf(stderr,"ridx1(%d): %d (%.1f: %.3f %.3f %.3f (%.1f %.1f %.1f)(%.3f %.3f %.3f))\n",
    status,r_idx,z1,r1,r2,r3,phot[X],phot[Y],phot[Z],phot[PX],phot[PY],phot[PZ]);
#endif
      }

    if(r_idx>=S_DATA)
      return(RT_OUT_OF_SURFACE2);
    if(r_idx<0 && r3>R_LENS && fabs(z3-z1)<10.0)
      return(RT_OUT_OF_SURFACE2);
    if(r_idx<0)
      {
    r_idx=1;
    /* return(RT_OUT_OF_SURFACE1); */
      }

    /* iteration doesn't converge */
    if(iter>RT_SEARCH_ITER_MAX)
      return(-iter);

    iter++;

  } while(!Found);

  /*
   * check no intersections at a nearer point
   */
  v[X]=x3;
  v[Y]=y3;
  v[Z]=0.0;
  if(Kvdot(v,&phot[PX])<0)
    idx_dir=1;
  else
    idx_dir=-1;

  do {
    r_idx+=idx_dir;

    r1=lens_r[r_idx];
    z1=lens_z[r_idx];
    r2=lens_r[r_idx+1];
    z2=lens_z[r_idx+1];

    if(r1==r2)
      status=Kcheck_intersect_cylinder(r1, z1, r2, z2, phot, &tmp_dis);
    else
      status=Kcheck_intersect_cone(r1, z1, r2, z2, phot, &tmp_dis);

    if(!status)  /* intersection found */
      {
    tmp_x[X]=phot[PX]*tmp_dis+phot[X];
    tmp_x[Y]=phot[PY]*tmp_dis+phot[Y];
    tmp_x[Z]=phot[PZ]*tmp_dis+phot[Z];
    r3=TMath::Sqrt(tmp_x[X]*tmp_x[X]+tmp_x[Y]*tmp_x[Y]);

    if(r1==r2)
      {
        if(z2>z1)
          dir=1;
        else
          dir=-1;
        nnn[X]=dir*tmp_x[X]/r3;
        nnn[Y]=dir*tmp_x[Y]/r3;
        nnn[Z]=0.0;
      }
    else
      {
        a00=(z2-z1)/(r2-r1);
        /* normal vector */
        nnn[X]=a00*tmp_x[X];
        nnn[Y]=a00*tmp_x[Y];
        nnn[Z]=-r3;
        dbuf1=TMath::Sqrt(Kvdot(nnn,nnn));
        nnn[X]/=dbuf1;
        nnn[Y]/=dbuf1;
        nnn[Z]/=dbuf1;
      }

    if(Kvdot(nnn,&phot[PX])>=0) /* wrong intersection; back side */
      status=1;
    else                      /* update intersection point */
      {
        x3=tmp_x[X];
        y3=tmp_x[Y];
        z3=tmp_x[Z];
        *dis=tmp_dis;
      }
      }
  }while(!status);

#ifndef NO_DIFFRACTIVE
  if(ID==1 && phot[PZ]>0.0 && nnn[Z]!=0.0)
    {
      Kdiffract(ID, &phot[PX], x3, y3, z3, Zc, 1.0003, 1.0003, 
           phot[L]*1e-6, 1);
    }
  if(ID==6 && phot[PZ]<0.0 && nnn[Z]!=0.0)
    {
      Kdiffract(ID, &phot[PX], x3, y3, z3, Zc, 1.0003, 1.0003, 
           phot[L]*1e-6, 1);
    }
#endif

  /* refraction */
  if(phot[PZ]<0.0)
    {
      tmp=medium2_n;
      medium2_n=medium1_n;
      medium1_n=tmp;
    }
  flag=Krefract(&phot[PX], nnn, medium1_n, medium2_n);
#if 0
  if(flag==RT_REFLECTED) 
    return(RT_REFLECTION_DISCARDED); /* reflected photons are thrown away */
#endif

#ifndef NO_DIFFRACTIVE
  if(ID==1 && phot[PZ]<0.0 && nnn[Z]!=0.0)
    {
      Kdiffract(ID, &phot[PX], x3, y3, z3, Zc, 1.0003, 1.0003, 
           phot[L]*1e-6, 1);
    }
  if(ID==6 && phot[PZ]>0.0 && nnn[Z]!=0.0)
    {
      Kdiffract(ID, &phot[PX], x3, y3, z3, Zc, 1.0003, 1.0003, 
           phot[L]*1e-6, 1);
    }
#endif

  /* distance between injection point and the Surface */
  phot[T]+=(*dis)* medium1_n/C0;

  phot[X]=x3;
  phot[Y]=y3;
  phot[Z]=z3;

  return 0;
}

static Int_t KCheckIris(Double_t in_ph[6])
{
  Double_t x3, y3, dist, radius;

  dist=(Zc3-in_ph[Z])/in_ph[PZ];
  if(dist>=0.0)
    {
      x3=in_ph[PX]*dist+in_ph[X];
      y3=in_ph[PY]*dist+in_ph[Y];

      radius=TMath::Sqrt(x3*x3+y3*y3);
      if(radius>Rc3)
        return (RT_OUT_OF_IRIS);
    }

    return 0;
}

/*******************************************/
/*         main routine                    */
/*******************************************/

Int_t Ktrace(Double_t ph_in[], Double_t ph_out[])
     /* ph_in[0] : x0 mm */
     /* ph_in[1] : y0 mm */
     /* ph_in[2] : z0 mm */
     /* ph_in[3] : l     */
     /* ph_in[4] : m     */
     /* ph_in[5] : n     */
     /* ph_in[6] : lambda nm */
     /* ph_in[7] : time ns */
{
  Int_t    status=0, n_int, WithinOptics;
  Int_t    PrevPoint, CurPoint, NextPoint;
  Double_t nn,kk;
  Double_t AbsorbCoeff, TransProb;
  Double_t dis;

  ph_out[L]=ph_in[L];


  /* Opt. index  for this lambda */ 
  nn=KGet_n_Acryl(ph_in[L]);
  kk=KGet_k_Acryl(ph_in[L]);

  if(ph_in[Z]>=Sz6 && ph_in[PZ]<0.0)
    {
      PrevPoint=SURFACE7;
      CurPoint=SURFACE7;
      NextPoint=SURFACE6;
    }
  else
    {
      PrevPoint=SURFACE0;
      CurPoint=SURFACE0;
      NextPoint=SURFACE1;
    }
  n_int=0;
  WithinOptics=1;

  while(n_int++ < RT_INTERACTION_MAX && WithinOptics)
    {
      PrevPoint=CurPoint;
      CurPoint=NextPoint;

      switch(NextPoint)
    {
    case SURFACE0:
    case SURFACE7:
      WithinOptics=0;
      break;

    case SURFACE1:
      status=Klens_interaction(1, ph_in, Zc1, Rc1, sr1, sz1, n_s1,
                  REFRACTIVE_IDX_MEDIA,nn,&dis);
#ifdef DEBUG
      printf("1 %f %f %f\n",ph_in[X],ph_in[Y],ph_in[Z]);
#endif
      if(ph_in[PZ]>0.0)
        NextPoint=SURFACE2;
      else
        NextPoint=SURFACE0;
      break;

    case SURFACE2:
      status=Klens_interaction(2, ph_in, Zc2, Rc2, sr2, sz2, n_s2,
                  nn,REFRACTIVE_IDX_MEDIA, &dis);
#ifdef DEBUG
      printf("2 %f %f %f\n",ph_in[X],ph_in[Y],ph_in[Z]);
#endif
      if(ph_in[PZ]>0.0)
        NextPoint=SURFACE5;
      else
        NextPoint=SURFACE1;
      break;

    case SURFACE5:
      status=Klens_interaction(5, ph_in, Zc5, Rc5, sr5, sz5, n_s5,
                  REFRACTIVE_IDX_MEDIA, nn, &dis);
#ifdef DEBUG
      printf("5 %f %f %f\n",ph_in[X],ph_in[Y],ph_in[Z]);
#endif
      if(ph_in[PZ]>0.0)
        NextPoint=SURFACE6;
      else
        NextPoint=SURFACE2;
      break;

    case SURFACE6:
      status=Klens_interaction(6, ph_in, Zc6, Rc6, sr6, sz6, n_s6,
                  nn, REFRACTIVE_IDX_MEDIA, &dis);
#ifdef DEBUG
      printf("6 %f %f %f\n",ph_in[X],ph_in[Y],ph_in[Z]);
#endif
      if(ph_in[PZ]>0.0)
        NextPoint=SURFACE7;
      else
        NextPoint=SURFACE5;
      break;
    }
#ifdef DEBUG
      fprintf(stderr,"%d %.1f %.1f %.1f  %.3f %.3f %.3f\n",
          -CurPoint/1000,
          ph_in[X],ph_in[Y],ph_in[Z],ph_in[PX],ph_in[PY],ph_in[PZ]);
#endif
    // thea debug
    // even if it dies, I want to know where it does
    ph_out[X]  = ph_in[X];
    ph_out[Y]  = ph_in[Y];
    ph_out[Z]  = ph_in[Z];
    ph_out[PX] = ph_in[PX];
    ph_out[PY] = ph_in[PY];
    ph_out[PZ] = ph_in[PZ];
    ph_out[T]  = ph_in[T];

    if (status == RT_NO_INTERSECTION1 || status == RT_HIT_SIDE_WALL) {
        // track the photons up to the borders
        if ( TMath::Sqrt(ph_in[X]*ph_in[X]+ph_in[Y]*ph_in[Y]) > R_WALL){
            //ph is outside, invert it's dir for tracing
            ph_in[PX] = -ph_in[PX];
            ph_in[PY] = -ph_in[PY];
            ph_in[PZ] = -ph_in[PZ];
        }
        Double_t tmp_dis;
        if (!Kcheck_intersect_cylinder(R_WALL, Sz1, R_WALL, Sz6, ph_in, &tmp_dis)){
            // the photon crosses the boundaries of the optics
            ph_out[X]=ph_in[X]+ph_in[PX]*tmp_dis;
            ph_out[Y]=ph_in[Y]+ph_in[PY]*tmp_dis;
            ph_out[Z]=ph_in[Z]+ph_in[PZ]*tmp_dis;
        }
        return 0;
    }

    if(status)
        return (status+CurPoint);

#ifndef NO_ABSORPTION
      /* absorption between S1 and S2 or between S5 and S6 */
      if((PrevPoint==SURFACE1 && CurPoint==SURFACE2) ||
     (PrevPoint==SURFACE2 && CurPoint==SURFACE1) ||
     (PrevPoint==SURFACE5 && CurPoint==SURFACE6) ||
     (PrevPoint==SURFACE6 && CurPoint==SURFACE5))
   {
      AbsorbCoeff=kk*TMath::Pi()*4.0/(ph_in[L]*1.e-6);
      TransProb=exp(-AbsorbCoeff*dis);
          TRandom *rndm = EsafRandom::Get();
      if(rndm->Rndm()>TransProb)
        return(RT_ABSORBED+CurPoint);
    }
#endif
    if((CurPoint==SURFACE2 && NextPoint==SURFACE5) ||
     (CurPoint==SURFACE5 && NextPoint==SURFACE2))
        if(KCheckIris(ph_in)<0){
            Double_t tmp_dis, radius, x3, y3;
            
            ph_out[PX] = ph_in[PX];
            ph_out[PY] = ph_in[PY];
            ph_out[PZ] = ph_in[PZ];
            ph_out[T]  = ph_in[T];

            tmp_dis=(Zc3-ph_in[Z])/ph_in[PZ];
            x3=ph_in[PX]*tmp_dis+ph_in[X];
            y3=ph_in[PY]*tmp_dis+ph_in[Y];
            radius = TMath::Sqrt(x3*x3+y3*y3);
            if ( radius < R_WALL) {
                // intersection with the pupil
                ph_out[X]=x3;
                ph_out[Y]=y3;
                ph_out[Z]=Zc3;
                return RT_OUT_OF_IRIS;
            } else {
                // the photon hit the wall, not the iris
                Kcheck_intersect_cylinder(R_WALL, ph_in[Z], R_WALL, Zc3, ph_in, &tmp_dis);
                ph_out[X]=ph_in[X]+ph_in[PX]*tmp_dis;
                ph_out[Y]=ph_in[Y]+ph_in[PY]*tmp_dis;
                ph_out[Z]=ph_in[Z]+ph_in[PZ]*tmp_dis;
                return RT_HIT_SIDE_WALL;
            }

        }
    }

  if(n_int>=RT_INTERACTION_MAX)
    return (RT_STRAY_LIGHT);

  return 0;
}