/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2009, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Lens
*
* %Identification
*
* Written by: C. Monzat/E. Farhi/S. Desert/G. Euzen
* Date: 2009
* Origin: ILL/LLB
*
* Refractive lens with absorption, incoherent scattering and surface imperfection.
*
* %Description
* Refractive Lens with absorption, incoherent scattering and surface imperfection.
* Geometry may be:
*   spherical (use r1 and r2 to specify radius of curvature),
*   planar    (use phiy1 and phiy2 to specify rotation angle of plane w.r.t beam)
*   parabolic (use focus1 and focus2 as optical focal length).
*
* Optionally, you can specify the 'geometry' parameter as a OFF/PLY file name.
*   The complex geometry option handles any closed non-convex polyhedra.
*   It computes the intersection points of the neutron ray with the object
*   transparently, so that it can be used like a regular sample object.
*   It supports the PLY, OFF and NOFF file format but not COFF (colored faces).
*   Such files may be generated from XYZ data using:
*     qhull < coordinates.xyz Qx Qv Tv o > geomview.off
*   or
*     powercrust coordinates.xyz
*   and viewed with geomview or java -jar jroff.jar (see below).
*
* The lens cross-section is seen as a 2*radius disk from the beam Z axis, except
*   when a 'geometry' file is given.
* Usually, you should stack more than one of these to get a significant effect
* on the neutron beam, so-called 'compound refractive lens'.
*
* The focal length for N lenses with focal 'f' is f/N, where f=R/(1-n)
*   and R = r/2   for a spherical lens with curvature radius 'r'
*       R = focus for a parabolic lens with focal of the parabola
*   and n = sqrt(1-(lambda*lambda*rho*bc/PI)) with
*     bc  = sqrt(fabs(sigma_coh)*100/4/PI)*1e-5
*     rho = density*6.02214179*1e23*1e-24/weight
*
* Common materials: Should have high coherent, and low incoherent and absorption cross sections
*   Be:            density=1.85,  weight=9.0121, sigma_coh=7.63,  sigma_inc=0.0018,sigma_abs=0.0076
*   Pb:            density=11.115,weight=207.2,  sigma_coh=11.115,sigma_inc=0.003, sigma_abs=0.171
*     Pb206:                                     sigma_coh=10.68, sigma_inc=0    , sigma_abs=0.03
*     Pb208:                                     sigma_coh=11.34, sigma_inc=0    , sigma_abs=0.00048
*   Zr:            density=6.52,  weight=91.224, sigma_coh=6.44,  sigma_inc=0.02,  sigma_abs=0.185
*     Zr90:                                      sigma_coh=5.1,   sigma_inc=0    , sigma_abs=0.011
*     Zr94:                                      sigma_coh=8.4,   sigma_inc=0    , sigma_abs=0.05
*   Bi:            density=9.78,  weight=208.98, sigma_coh=9.148, sigma_inc=0.0084,sigma_abs=0.0338
*   Mg:            density=1.738, weight=24.3,   sigma_coh=3.631, sigma_inc=0.08,  sigma_abs=0.063
*   MgF2:          density=3.148, weight=62.3018,sigma_coh=11.74, sigma_inc=0.0816,sigma_abs=0.0822
*   diamond:       density=3.52,  weight=12.01,  sigma_coh=5.551, sigma_inc=0.001, sigma_abs=0.0035
*   Quartz/silica: density=2.53,  weight=60.08,  sigma_coh=10.625,sigma_inc=0.0056,sigma_abs=0.1714
*   Si:            density=2.329, weight=28.0855,sigma_coh=2.1633,sigma_inc=0.004, sigma_abs=0.171
*   Al:            density=2.7,   weight=26.98,  sigma_coh=1.495, sigma_inc=0.0082,sigma_abs=0.231
*   perfluoropolymer(PTFE/Teflon/CF2):
*                  density=2.2,   weight=50.007, sigma_coh=13.584,sigma_inc=0.0026,sigma_abs=0.0227
*   Organic molecules with C,O,H,F
*
* Example: Lens(r1=0.025,r2=0.025, thickness=0.001,radius=0.0150)
*
* %BUGS
* parabolic shape is not fully validated yet, but should do still...
*
* %Parameters
* INPUT PARAMETERS:
* r1: [m]            radius of the first circle describing the lens. r>0 means concave face, r<0 means convex, r=0 means plane
* r2: [m]            radius of the second circle describing the lens r>0 means concave face, r<0 means convex, r=0 means plane
* focus1: [m]        focal of the first parabola describing the lens
* focus2: [m]        focal of the second parabola describing the lens
* phiy1: [deg]       angle of plane1 (r1=0) around y vertical axis
* phiy2: [deg]       angle of plane2 (r2=0) around y vertical axis
* thickness: [m]     thickness of the lens between its two  surfaces
* radius: [m]        radius of the lens section, e.g. beam size.
* density:  [g/cm3] density of the material for the lens
* weight: [g/mol]    molar mass of the material
* focus_aw: [deg]    vertical angle to forward focus after scattering
* focus_ah: [deg]    horizontal angle to forward focus after scattering
* sigma_coh: [barn]  coherent cross section
* sigma_inc: [barn]  incoherent cross section
* sigma_abs: [barn]  thermal absorption cross section
* p_interact: [1]    MC Probability for scattering the ray; otherwise transmit
* RMS: [Angs]        root mean square roughness of the surface
* geometry: [str]    Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust
*
* %L
* M. L. Goldberger et al, Phys. Rev. 71, 294 - 310 (1947)
* %L
* Sears V.F. Neutron optics. An introduction to the theory of neutron optical phenomena and their applications. Oxford University Press, 1989.
* %L
* H. Park et al. Measured operationnal neutron energies of compound refractive lenses. Nuclear Instruments and Methods B, 251:507-511, 2006.
* %L
* J. Feinstein and R. H. Pantell. Characteristics of the thick, compound refractive lens. Applied Optics, 42 No. 4:719-723, 2001.
* %L
* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a>
* %L
* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a>
* %L
* <a href="http://qhull.org">qhull</a>
* %L
* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a>
* %E
*******************************************************************************/

DEFINE COMPONENT Lens

SETTING PARAMETERS (r1=0,r2=0,focus1=0,focus2=0,phiy1=0,phiy2=0,
thickness=0.001,radius=0.015,
sigma_coh=11.74,sigma_inc=0.0816,sigma_abs=0.0822,
density=3.148,weight=62.3018,
p_interact=0.1,focus_aw=10,focus_ah=10,RMS=0, string geometry=0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
/* support for OFF/PLY geometry */
%include "read_table-lib"
%include "interoff-lib"

/* Lens_roughness: function to rotate normal vector around axis for roughness
*                 with specified tilt angle (deg)
   * RETURNS: rotated nx,ny,nz coordinates
   */
#pragma acc routine seq
  void Lens_roughness(double *nx, double *ny, double *nz, double tilt
#ifdef OPENACC
    , _class_particle* _particle
#endif
    )
  {
    double nt_x, nt_y, nt_z;  /* transverse vector */
    double n1_x, n1_y, n1_z;  /* normal vector (tmp) */

    /* normal vector n_z = [ 0,1,0], n_t = n x n_z; */
    vec_prod(nt_x,nt_y,nt_z, *nx,*ny,*nz, 0,1,0);

    /* rotate n with angle wav_z around n_t -> n1 */
    tilt  *= DEG2RAD/(sqrt(8*log(2)))*randnorm();
    rotate(n1_x,n1_y,n1_z, *nx,*ny,*nz, tilt, nt_x,nt_y,nt_z);

    /* rotate n1 with angle phi around n -> nt */
    rotate(nt_x,nt_y,nt_z, n1_x,n1_y,n1_z, 2*PI*rand01(), *nx,*ny,*nz);

    *nx=nt_x;
    *ny=nt_y;
    *nz=nt_z;
  }

/* parabola_intersect: Calculate intersection between a line and a parabola with
*  axis along z, focal length f, centered at (0,0,0)
 * RETURNS 0 when no intersection is found
 *      or 1 in case of intersection with resulting times t0 and t1 */
#pragma acc routine seq
  int parabola_intersect(double *t0, double *t1, double x, double y, double z,
                 double vx, double vy, double vz, double f)
  {
     /* equation of line: (x,y,z) = (vx,vy,vz)*t+(x,y,z)_0 */
     /* equation of parabola:  z  = (x^2+y^2)/4/f
        that is: 4*f*(vz*t+z0) = (vx*t+x0)^2 + (vy*t+y0)^2 */

     double A = vx*vx+vy*vy;
     double B = 2*vx*x+2*vy*y-4*vz*f;
     double C = x*x+y*y-4*f*z;

     return(solve_2nd_order(t0, t1, A,B,C));
  }

  /* Lens_intersect: function to compute intersection with lens surface
   * of section 2*radius (beam size)
   * which can be
   *   type="spherical" (arg=curvature radius),
   *   type="parabolic" (arg=focus) or
   *   type="planar"    (arg=phi around y).
   * The Surface intersection along the Z axis is 'shift' (e.g. +/-thickness/2)
   * RETURNS: intersection time to surface, or negative for no intersection
   *          neutron must then be propagated on the surface with PROP_DT,
   *          and checked for actual intersection.
   */
#pragma acc routine seq
  double Lens_intersect(double x,double y,double z,
    double vx,double vy,double vz,
    double shift, double radius,
    int type, double arg)
  {
    double t0, t1;

    if (!vz || !type || !radius) return -1;
    if (type==1 && arg) {   /* spherical */
      if (sphere_intersect(&t0,&t1, x,y,z+shift+arg, vx,vy,vz, fabs(arg))) {
        return ( arg > 0 ? t1 : t0 ); /* concave : convex surface */
      } else
        return(-1);
    } else if (type==2 && arg) { /* parabolic */
      if (parabola_intersect(&t0, &t1, x,y,z+shift, vx,vy,vz, arg)) {
        if (t0 < 0 && 0 < t1) return(t1);
        if (t1 < 0 && 0 < t0) return(t0);
        return(t1 < t0 ? t1 : t0);
      } else
        return(-1);
      /* end parabola */
    } else if (type==3) {    /* planar */
      if (plane_intersect(&t0, x,y,z+shift, vx,vy,vz, -sin(arg),0,-cos(arg), 0,0,0)>0)
        return(t0);
      else
        return(-1);
    } else  /* unsupported geometry */
      return(-1);
  }
%}

DECLARE
%{
  double rho;
  double bc;
  double my_s;
  double my_a_v;
  double mean_n;
  double events;
  off_struct offdata;
%}

INITIALIZE
%{
  /* density: Number of atoms by AA^3, pow(10,-24) stands for conversion from cm^3 to AA^3 */
  if (density<=0 || weight<=0)
    exit(printf("Lens: %s: FATAL: invalid material density or molar weight: density=%g weight=%g\n",
      NAME_CURRENT_COMP, density, weight));
  rho= density*6.02214179*1e23*1e-24/weight;
  if (sigma_coh==0)
    exit(printf("Lens: %s: FATAL: invalid material coherent cross section: sigma_coh=%g\n",
      NAME_CURRENT_COMP, sigma_coh));
  bc=sqrt(fabs(sigma_coh)*100/4/PI)*1e-5;      /* bound coherent cross section */
  if (sigma_coh<0) bc *= -1;
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    #ifndef USE_OFF
    fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n");
    exit(-1);
    #else
    /* init the OFF/PLY without centering and scaling */
	  if (!off_init(geometry, 0, 0, 0, 1, &offdata))
      exit(printf("Lens: %s: FATAL: could not initialize geometry file %s [OFF/PLY]\n",
      NAME_CURRENT_COMP, geometry));
    #endif
  }

  phiy1 *= DEG2RAD; /* planes orientation */
  phiy2 *= DEG2RAD;

  focus_aw *= DEG2RAD;
  focus_ah *= DEG2RAD;
  my_s  = rho * 100 *(sigma_inc+sigma_coh);
  my_a_v= rho * 100 * sigma_abs;
%}

TRACE
%{

double n; /* refractive index */
double v, lambda; /* neutron speed and wavelength */
double dt;        /* time to intersection */
double theta_RMS;
double nx=0, ny=0, nz=0; /* normal vector to surface */
double ax, ay, az; /* temporary vector for surface refraction */
double alpha1, beta1, theta1, theta2; /* angles for refraction computation */

#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
#endif
#else
#define thread_offdata offdata
#endif

/* iterate two faces intersection until no more is possible */
do {
  /* ======================== First face of the lens ========================== */
  /* determine intersection time */
  nx=ny=nz=0;
  #ifdef USE_OFF
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    double t1=0,t0=0;
    Coords n, n0, n1;
    int    intersect=off_intersect  (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
    if (!intersect) dt=-1;
    else if (t0 < 0 && 0 < t1)      { dt = t1; n=n1; }
    else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
    else {
      if (t1 < t0) { dt=t1; n=n1; }
      else  { dt=t0; n=n0; }
    }
    coords_get(n, &nx, &ny, &nz);
  } else
  #endif
  if (r1 && !focus1) {
    dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 1, r1);
  } else if (!r1 && focus1){
      if (focus1>0){
          dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 2, -focus1);
      }else{
          dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2 + radius*radius/(4*fabs(focus1)) , radius, 2, -focus1);
      }
  } else {
    dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 3, phiy1);
  }
  if (dt <= 0) break; /* no intersection: exit from main loop */

  /* propagate to surface (1) */
  PROP_DT(dt);

  /* check for lens cross section (radius) */
  if (radius && x*x+y*y > radius*radius) break; /* outside from cross section: exit from main loop */
  SCATTER;

  v=sqrt(vx*vx+vy*vy+vz*vz);

  if (v) lambda=3956.0032/v; else ABSORB;

  /* compute refractive index */
  /* without inc nor abs: n  = sqrt(1-(lambda*lambda*rho*bc/PI));   */
  /* M. L. Goldberger et al, Phys. Rev. 71, 294 - 310 (1947) */
  /* alpha1 = (sigma_inc+sigma_abs*2200/v)/2/lambda;
  n      = 1-lambda*lambda*rho/2/PI*sqrt(bc*bc-alpha1*alpha1); */
  n   = sqrt(1-(lambda*lambda*rho*bc/PI));

#pragma acc atomic
  mean_n += n*p;
#pragma acc atomic
  events += p;

  theta_RMS=atan(2*RMS/lambda);           /* cone angle from RMS roughness */

  /* compute normal vector (1) when not given by OFF/PLY */
  if (!nx && !ny && !nz) {
    dt = 0;
    if (r1 && !focus1) {
      double sign=r1/fabs(r1);
      nx=-sign*x;
      ny=-sign*y;
      nz=sign*(-z-thickness/2-r1);
    } else if (!r1 && focus1) {
      nx=-x/2/focus1;
      ny=-y/2/focus1;
      nz=-1;
    } else {
      nx=-sin(phiy1);
      ny=0;
      nz=-cos(phiy1);
    }
  }
  /* tilt normal vector for roughness, in cone theta_RMS */
  if (RMS>0) Lens_roughness(&nx, &ny, &nz, theta_RMS
#ifdef OPENACC
			    , _particle
#endif
			    );

  /* compute incoming angles w.r.t surface */
  NORM(nx,ny,nz);
  vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz);
  /* if n and v are parallel - no refraction.*/
  if (ax!=0 || ay!=0 || az!=0){
      theta1 = atan2(sqrt(ax*ax+ay*ay+az*az),scalar_prod(nx,ny,nz,-vx,-vy,-vz));

      /* Fresnel formula for refraction */
      theta2 = asin(sin(theta1)/n);

      /* compute new velocity after refraction */
      alpha1 = (sin(theta2))/(sin(theta1));
      beta1  = (alpha1*v*cos(theta1))-(v*cos(theta2));
      vx = beta1*nx + alpha1*vx;
      vy = beta1*ny + alpha1*vy;
      vz = beta1*nz + alpha1*vz;
  }
  /* ======================== Second face of the lens ========================= */
  /* determine intersection time to go to second surface */
  nx=ny=nz=0;
  #ifdef USE_OFF
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    double t1=0,t0=0;
    Coords n, n0, n1;
    int    intersect=off_intersect  (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
    if (!intersect) dt=-1;
    else if (t0 < 0 && 0 < t1)      { dt = t1; n=n1; }
    else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
    else {
      if (t1 < t0) { dt=t1; n=n1; }
      else  { dt=t0; n=n0; }
    }
    coords_get(n, &nx, &ny, &nz);
  } else
  #endif
  if (r2 && !focus2)
    dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 1, -r2);
  else if (!r2 && focus2){
      if (focus2>0){
          dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 2, focus2);
      }else{
          dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2 - radius*radius/(4*fabs(focus2)) , radius, 2, focus2);
      }
  }else{
    dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 3, phiy2);
  }
  if (dt <= 0) break; /* no intersection: exit from main loop */

  /* =============== Absorption and scattering between the two faces =========== */
  double my_a,my_t,p_trans,p_scatt,mc_trans,mc_scatt,ws,d_path;
  double p_mult=1;
  int    flag=0;
  double l_i=0;
  double solid_angle=0;

  my_a = my_a_v*(2200/v);
  my_t = my_a + my_s;
  ws = my_s/my_t;  /* (inc+coh)/(inc+coh+abs) */

  d_path = v * dt;

  /* Proba of transmission along length d_path */
  p_trans = exp(-my_t*d_path);
  p_scatt = 1 - p_trans;

  flag = 0; /* flag used for propagation to exit point before ending */
  /* are we next to the exit ? probably no scattering (avoid rounding errors) */
  if (my_s*d_path <= 4e-7) {
    flag = 1;           /* No interaction before the exit */
  }

  /* force a given fraction of the beam to scatter */
  if (p_interact>0 && p_interact<=1) {
    /* we force a portion of the beam to interact */
    /* This is used to improve statistics on single scattering (and multiple) */
    mc_trans = 1-p_interact;
  } else {
    mc_trans = p_trans; /* 1 - p_scatt */
  }

  mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */
  if (mc_scatt <= 0 || mc_scatt>1) flag=1;

  /* MC choice: Interaction or transmission ? */
  if (!flag && mc_scatt > 0 && (mc_scatt >= 1 || (rand01()) < mc_scatt)) { /* Interaction neutron/sample */
    p_mult *= ws; /* Update weight ; account for absorption and retain scattered fraction */
    /* we have chosen portion mc_scatt of beam instead of p_scatt, so we compensate */
    p_mult *= fabs(p_scatt/mc_scatt); /* lower than 1 */
  } else {
    flag = 1; /* Transmission : no interaction neutron/sample */
    p_mult *= fabs(p_trans/mc_trans);  /* attenuate beam by portion which is scattered (and left along) */
  }

  if (flag) { /* propagate directly to secound surface */
    if (!isnan(p_mult)) {
	p *= p_mult; /* apply absorption correction */
    } else {
        ABSORB;
    }
  }
  else
  {
    double a;
    if (my_t*d_path < 1e-6){
      /* For very weak scattering, use simple uniform sampling of scattering
      point to avoid rounding errors. */
      dt = rand0max(d_path); /* length */
    } else {
      a = rand0max((1 - exp(-my_t*d_path)));
      dt = -log(1 - a) / my_t; /* length */
    }

    l_i = dt;/* Penetration in sample: scattering+abs */
    dt /= v; /* Time from present position to scattering point */
    PROP_DT(dt); /* Point of scattering */

    randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,0,0,1, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
    NORM(vx, vy, vz);

    vx*=v;
    vy*=v;
    vz*=v;

    p_mult *= solid_angle/4/PI;
    if (!isnan(p_mult)) {
	p *= p_mult;
    } else {
        ABSORB;
    }
    SCATTER;

    /* recompute new intersection time to go to second surface (velocity changed during scattering event) */
    nx=ny=nz=0;
    #ifdef USE_OFF
    if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
      double t1=0,t0=0;
      Coords n, n0, n1;
      int    intersect=off_intersect  (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
      if (!intersect) dt=-1;
      else if (t0 < 0 && 0 < t1)      { dt = t1; n=n1; }
      else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
      else {
        if (t1 < t0) { dt=t1; n=n1; }
        else  { dt=t0; n=n0; }
      }
      coords_get(n, &nx, &ny, &nz);
    } else
    #endif
    if (r2 && !focus2)
      dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 1, -r2);
    else if (!r2 && focus2)
      if (focus2>0){
          dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 2, focus2);
      }else{
          dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2 - radius*radius/(4*fabs(focus2)) , radius, 2, focus2);
      }
    else
      dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 3, phiy2);

    if (dt <= 0) break; /* no intersection: exit from main loop */

  } /* end scattering handling in material*/

  /* propagate to surface */
  PROP_DT(dt);

  /* check for lens cross section (radius) */
  if (radius && x*x+y*y > radius*radius) break; /* outside from cross section: exit from main loop */
  SCATTER;

  v=sqrt(vx*vx+vy*vy+vz*vz);

  /* compute normal vector (2) when not given by OFF/PLY */
  if (!nx && !ny && !nz) {
    dt = 0;
    if (r2 && !focus2) {
      double sign=r2/fabs(r2);
      nx=sign*x;
      ny=sign*y;
      nz=sign*(z-thickness/2-r2);
    } else if (!r2 && focus2){
      nx=x/2/focus2;
      ny=y/2/focus2;
      nz=-1;
    } else {
      nx=-sin(phiy2);
      ny=0;
      nz=-cos(phiy2);
    }
  }
  /* tilt normal vector for roughness, in cone theta_RMS */
  if (RMS>0) Lens_roughness(&nx, &ny, &nz, theta_RMS
#ifdef OPENACC
			    , _particle
#endif
			    );

  /* compute incoming angles w.r.t surface */
  NORM(nx,ny,nz);
  vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz);
  /* if n and v are parallel - no refraction.*/
  if (ax!=0 || ay!=0 || az!=0){
      theta1 = atan2(sqrt(ax*ax+ay*ay+az*az),scalar_prod(nx,ny,nz,-vx,-vy,-vz));

      /* Fresnel formula for refraction */
      theta2=asin(sin(theta1)*n);

      /* compute new velocity after refraction */
      alpha1 = (sin(theta2))/(sin(theta1));
      beta1 = (alpha1*v*cos(theta1))-(v*cos(theta2));
      vx = beta1*nx + alpha1*vx;
      vy = beta1*ny + alpha1*vy;
      vz = beta1*nz + alpha1*vz;
  }
} while (dt>0); /* loop until no more intersection */

%}

FINALLY %{
  /* compute focal length */
  double f1=0, f2=0, f=0;
  double n = mean_n / events;

  if (n != 1) {
    if (r1 && !focus1)      f1 = r1/2/(1-n);
    else if (!r1 && focus1) f1 = focus1/(1-n);

    if (r2 && !focus2)      f2 = r2/2/(1-n);
    else if (!r2 && focus2) f2 = focus2/(1-n);
  }

  if (f1) f+= 1/f1;
  if (f2) f+= 1/f2;

  if (f)  {
    f = 1/f;
    printf("Lens: %s: focal length f=%g [m]. Focal length for N lenses is f/N.\n", NAME_CURRENT_COMP, f);
  }

%}

MCDISPLAY
%{
  magnify("xy");
  double theta1,theta00=0,theta01=0,eps,eps2,theta_line,distance,height,height2,dist_parab1,dist_parab2;
  height=radius/2;
  height2=radius/2;
  dist_parab1=0.0;
  dist_parab2=0.0;

  /* draw circles to describe the 2 surfaces */
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {	/* OFF file */
    off_display(offdata);
  }
  else {
    if (r1 && !focus1) { /* sphere1 */
      theta00=asin(fabs((radius/2)/r1));
      if (r1<0 && r2<=0 && (-r1+r1*cos(theta00)>thickness/2) ){
        theta00=acos((fabs(r1)-thickness/2)/fabs(r1));
        height=fabs(r1)*sin(theta00);
      }

      theta1=-theta00;
      eps=theta00/4;
      eps2=2*PI/4;
      while (theta1<0) {
        circle("xy",0,0,-thickness/2-r1+r1*cos(theta1),fabs(r1)*sin(theta1));
        theta1+=eps;
        theta_line=0;
        while (theta_line<2*PI){
          line(fabs(r1)*sin(theta1-eps)*cos(theta_line),fabs(r1)*sin(theta1-eps)*sin(theta_line),-thickness/2-r1+r1*cos(theta1-eps),fabs(r1)*sin(theta1)*cos(theta_line),fabs(r1)*sin(theta1)*sin(theta_line),-thickness/2-r1+r1*cos(theta1));
          theta_line+=eps2;
        }
      }

    } else if (!r1 && focus1) { /* parabola1 */
      if (focus1>0){
        dist_parab1=-radius*radius/(4*focus1);
      } else {
        dist_parab1=0.0;
      }

      distance=-(radius*radius)/(4*focus1);
      eps=-distance/4;
      eps2=2*PI/4;
      while (focus1*distance<0) {
        if (focus1>0){
          circle("xy",0,0,distance-thickness/2,sqrt((4*focus1*fabs(distance))));
          distance+=eps;
          theta_line=0;
          while (theta_line<2*PI){
             line(sqrt(4*fabs(focus1)*(fabs(distance-eps)))*cos(theta_line),sqrt(4*fabs(focus1)*(fabs(distance-eps)))*sin(theta_line),distance-eps-thickness/2,sqrt(4*fabs(focus1)*fabs(distance))*cos(theta_line),sqrt(4*fabs(focus1)*fabs(distance))*sin(theta_line),distance-thickness/2);
          theta_line+=eps2;
          }
        } else {
          theta_line=0;
          double rp,rpe,zz;
          rp=sqrt(4*fabs(focus1)*fabs(distance));
          zz=-thickness/2-radius*radius/(4*fabs(focus1))+distance;
          circle("xy",0,0,zz,rp);
          distance+=eps;
          zz=-thickness/2-radius*radius/(4*fabs(focus1))+distance;
          rpe=sqrt(4*fabs(focus1)*fabs(distance));

          theta_line=0;
          while (theta_line<2*PI){
            line(rpe*cos(theta_line),rpe*sin(theta_line),zz, rp*cos(theta_line), rp*sin(theta_line), zz-eps);
            theta_line+=eps2;
          }
        }
      }
    } else  { /* plane 1 */
      /* phiy1 is rotation angle around 'y' at z=-thickness. width/height is radius. */
      double x1= height2*cos(phiy1);
      double y1= height2;
      double z1=-height2*sin(phiy1)-thickness/2;
      multiline(5,  x1,  y1, z1,
                    x1, -y1, z1,
                   -x1, -y1,-z1,
                   -x1,  y1,-z1,
                    x1,  y1, z1);
    }

    if (r2 && !focus2) {  /* sphere 2 */
      theta01=asin(fabs((radius/2)/r2));
      if (r1<=0&&r2<0&&(r2-r2*cos(theta01)<-thickness/2)){
        theta01=acos((fabs(r2)-thickness/2)/fabs(r2));
        height2=fabs(r2)*sin(theta01);
      }
      theta1=PI-theta01;

      eps=(PI-theta1)/4;
      eps2=2*PI/4;
      while (theta1<PI) {
        circle("xy",0,0,thickness/2+r2+r2*cos(theta1),fabs(r2)*sin(theta1));
        theta1+=eps;
        theta_line=0;
        while (theta_line<2*PI) {
          line(fabs(r2)*sin(theta1-eps)*cos(theta_line),fabs(r2)*sin(theta1-eps)*sin(theta_line),thickness/2+r2+r2*cos(theta1-eps),fabs(r2)*sin(theta1)*cos(theta_line),fabs(r2)*sin(theta1)*sin(theta_line),thickness/2+r2+r2*cos(theta1));
          theta_line+=eps2;
        }
      }
    } else if (!r2 && focus2) { /* parabola 2 */
      if (focus2>0){
        dist_parab2=radius*radius/(4*focus2);
      }else{
        dist_parab2=0.0;
      }

      distance=(radius*radius)/(4*focus2);
      height2=sqrt((4*focus2*fabs(distance)));
      eps=-distance/4;
      eps2=2*PI/4;
      while (focus2*distance>0){
        if (focus2>0) {
          circle("xy",0,0,distance+thickness/2,sqrt((4*fabs(focus2)*fabs(distance))));
          distance+=eps;
          theta_line=0;
          while (theta_line<2*PI) {
            line(sqrt(4*focus2*(fabs(distance-eps)))*cos(theta_line),sqrt(4*focus2*(fabs(distance-eps)))*sin(theta_line),distance-eps+thickness/2,sqrt(4*focus2*fabs(distance))*cos(theta_line),sqrt(4*focus2*fabs(distance))*sin(theta_line),distance+thickness/2);
            theta_line+=eps2;
          }
        } else {
          double rp,rpe,zz;
          rp=sqrt(4*fabs(focus2)*fabs(distance));
          zz=thickness/2+(radius*radius)/(4*fabs(focus2))+distance;
          circle("xy",0,0,zz,rp);
          distance+=eps;
          zz=thickness/2+(radius*radius)/(4*fabs(focus2))+distance;
          rpe=sqrt(4*fabs(focus2)*fabs(distance));
          
          theta_line=0;
          while (theta_line<2*PI){
            line(rpe*cos(theta_line),rpe*sin(theta_line),zz, rp*cos(theta_line), rp*sin(theta_line), zz-eps);
            theta_line+=eps2;
          }
        }
      }
    } else { /* plane 2 */
      /* phiy2 is rotation angle around 'y' at z=+thickness. width/height is radius. */
      double x1= height2*cos(phiy2);
      double y1= height2;
      double z1=-height2*sin(phiy2)+thickness/2;
      multiline(5,  x1,  y1, z1,
                    x1, -y1, z1,
                   -x1, -y1,-z1,
                   -x1,  y1,-z1,
                    x1,  y1, z1);
    }

    /* draw outer containing cylinder */
    if (r1 && r2 && !focus1 && !focus2){
		  line(-fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),-fabs(r2)*sin(theta01),fabs(r2)*sin(theta01),thickness/2+r2-r2*cos(theta01));
		  line(fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
		  line(-fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),-fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
		  line(-fabs(r1)*sin(theta00),fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
	  }
	  else if (r1 && !r2) {
		  if (!focus2){
			  line(-height,0,-thickness/2-r1+r1*cos(theta00),-height,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(height,0,-thickness/2-r1+r1*cos(theta00),height,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
			  line(0,-height,-thickness/2-r1+r1*cos(theta00),0,-height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(0,height,-thickness/2-r1+r1*cos(theta00),0,height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
		  }
		  else
		  {
			  line(-fabs(r1)*sin(theta00),0,-thickness/2-r1+r1*cos(theta00),-radius/2,0,thickness/2+radius*radius/(16*focus2));
			  line(fabs(r1)*sin(theta00),0,-thickness/2-r1+r1*cos(theta00),radius/2,0,thickness/2+radius*radius/(16*focus2));
		  }
	  }
	  else if (r2 && !r1) {
		  if (!focus1){
			  line(-height2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-height2,0,thickness/2+r2-r2*cos(theta01));
			  line(height2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),height2,0,thickness/2+r2-r2*cos(theta01));
			  line(0,height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,height2,thickness/2+r2-r2*cos(theta01));
			  line(0,height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,height2,thickness/2+r2-r2*cos(theta01));
		  }
		  else
		  {
			  line(-height2,0,-thickness/2-radius*radius/(16*focus1),-height2,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(height2,0,-thickness/2-radius*radius/(16*focus1),height2,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
			  line(0,height2,-thickness/2-radius*radius/(16*focus1),0,height2,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(0,height2,-thickness/2-radius*radius/(16*focus1),0,height2,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
		  }
	  }
	  else {
		  if (!focus1 && !focus2){
			  line(-radius/2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-radius/2,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(radius/2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),radius/2,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
			  line(0,-radius/2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,-radius/2,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(0,radius/2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,radius/2,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
		  }
		  else if (!focus1){
			  line(-radius/2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-radius/2,0,thickness/2+dist_parab2);
			  line(radius/2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),radius/2,0,thickness/2+dist_parab2);
			  line(0,-radius/2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,-radius/2,thickness/2+dist_parab2);
			  line(0,radius/2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,radius/2,thickness/2+dist_parab2);
		  }
		  else if (!focus2){
			  line(-height,0,-thickness/2+dist_parab1,-height,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(height,0,-thickness/2+dist_parab1,height,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
			  line(0,-height,-thickness/2+dist_parab1,0,-height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
			  line(0,height,-thickness/2+dist_parab1,0,height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
		  }
		  else {
			  line(-radius,0,-thickness/2+dist_parab1,-radius,0,thickness/2+dist_parab2);
			  line(radius,0,-thickness/2+dist_parab1,radius,0,thickness/2+dist_parab2);
			  line(0,-radius,-thickness/2+dist_parab1,0,-radius,thickness/2+dist_parab2);
			  line(0,radius,-thickness/2+dist_parab1,0,radius,thickness/2+dist_parab2);
		  }
	  }
  } /* end of not a OFF/PLY */
%}

END