/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Incoherent
*
* %I
* Written by: Kim Lefmann and Kristian Nielsen
* Date: 15.4.98
* Origin: Risoe
* Modified by: Aziz Daoud-aladine, ISIS, 2007: Added option to handle a spherical sample shape
* Modified by: Peter Christiansen, Risoe: Added outgoing polarization: P' = 1/3*P-2/3P = -1/3P NB! As above multiple scattering is ignored .
* Modified by: Reynald Arnerin, ILL, 2008: Added option to handle a complex geometry (OFF files)
*
* Incoherent sample (such as Vanadium) sample, with quasielastic component OR or global energy transfer.
*
* %D
* A Double-cylinder shaped incoherent scatterer (like Vanadium)
*   with both elastic and quasielastic (Lorentzian) components.
* No multiple scattering (but approximation available). Absorption included.
* <b>Sample focusing:</b>
*   The area to scatter to is a disk of radius 'focus_r' situated at the target.
*   This target area may also be rectangular if specified focus_xw and focus_yh
*   or focus_aw and focus_ah, respectively in meters and degrees.
*   The target itself is either situated according to given coordinates (x,y,z),
*   or defined with the relative target_index of the component to focus
*   to (next is +1).
*   This target position will be set to its AT position. When targeting to
*   centered components, such as spheres or cylinders, define an Arm component
*   where to focus to.
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or any other shape
*   box/plate:       xwidth x yheight x zdepth (thickness=0)
*   hollow box/plate:xwidth x yheight x zdepth and thickness>0
*   cylinder:        radius x yheight (thickness=0)
*   hollow cylinder: radius x yheight and thickness>0
*   sphere:          radius (yheight=0 thickness=0)
*   hollow sphere:   radius and thickness>0 (yheight=0)
*   any shape:       geometry=OFF file
*
*   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 default size of the object depends of the OFF file data, but its
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* Example: Incoherent(radius=0.05,focus_r=0.035, pack=1, target_index=1)
*          Incoherent(geometry="socket.off",focus_r=0.035, pack=1, target_index=1)
*
* %P
* INPUT PARAMETERS:
* radius: [m]         Outer radius of sample in (x,z) plane 
* xwidth: [m]         Horiz. dimension of sample (bounding box if off file), as a width 
* yheight: [m]        Vert.  dimension of sample (bounding box if off file), as a height. A sphere shape is used when 0 and radius is set 
* zdepth: [m]         Depth of sample (bounding box if off file) 
* thickness: [m]      Thickness of hollow sample 
* focus_r: [m]        Radius of disk containing target. Use 0 for full space 
* target_index: [1]   Relative index of component to focus at, e.g. next is +1 
*
* Optional parameters:
* pack: [1]           Packing factor 
* p_interact: [1]     MC Probability for scattering the ray; otherwise transmit 
* focus_xw: [m]       horiz. dimension of a rectangular area 
* focus_yh: [m]       vert.  dimension of a rectangular area 
* focus_aw: [deg]     horiz. angular dimension of a rectangular area 
* focus_ah: [deg]     vert.  angular dimension of a rectangular area 
* sigma_abs: [barns]  Absorption cross section pr. unit cell at 2200 m/s 
* sigma_inc: [barns]  Incoherent scattering cross section pr. unit cell 
* Vc: [AA^3]          Unit cell volume
* f_QE: [1]           Fraction of quasielastic scattering (rest is elastic) 
* gamma: [1]          Lorentzian width of quasielastic broadening (HWHM) 
* Etrans: [meV]       Global energy-transfer, for use in inelastic settings
* deltaE: [meV]       Width in energy around Etrans, for use in inelastic settings
* 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 
* concentric: [1]     Indicate that this component has a hollow geometry and may contain other components. It should then be duplicated after the inside part (only for box, cylinder, sphere) 
*
* order: [1]          Limit multiple scattering up to given order 0 means all (default), 1 means single, 2 means double, ... 
* target_x: []        
* target_y: [m]       position of target to focus at 
* target_z: []        
*
* %L
* <a href="http://www.ncnr.nist.gov/resources/n-lengths/">Cross sections for single elements</a>
* %L
* <a href="http://www.ncnr.nist.gov/resources/sldcalc.html>Cross sections for compounds</a>
* %L
* <a href="http://www.webelements.com/">Web Elements</a>
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/tests/Incoherent/">Test
* results</A> (not up-to-date).
* %L
* The test/example instrument <a href="../examples/vanadium_example.instr">vanadium_example.instr</a>.
* %L
* The test/example instrument <a href="../examples/QENS_test.instr">QENS_test.instr</a>.
* %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 Incoherent

SETTING PARAMETERS (string geometry=0, radius=0, xwidth=0, yheight=0, zdepth=0,
thickness=0, nx=0, ny=1, nz=0,
target_x = 0, target_y = 0, target_z = 0, focus_r = 0,
focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, int target_index=0,
pack = 1, p_interact=1, f_QE=0, gamma=0, Etrans=0,deltaE=0,
sigma_abs=5.08, sigma_inc=5.08, Vc=13.827, concentric=0, order=0)

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

SHARE
%{
%include "read_table-lib"
%include "interoff-lib"
struct StructVarsInc
{
double  sigma_a; /* Absorption cross section per atom (barns) */
    double  sigma_i; /* Incoherent scattering cross section per atom (barns) */
    double  rho;     /* Density of atoms (AA-3) */
    double  my_s;
    double  my_a_v;
    int     shape;       /* 0 cylinder, 1 box, 2 sphere, 3 OFF file */
    double  aw,ah;       /* rectangular angular dimensions */
    double  xw,yh;       /* rectangular metrical dimensions */
    double  tx,ty,tz;    /* target coords */
    };
%}

DECLARE
%{
  struct StructVarsInc VarsInc;
  off_struct offdata;
%}

INITIALIZE
%{
  VarsInc.shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape  */
  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
    if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
      VarsInc.shape=3; thickness=0; concentric=0;
    }
    #endif
  }
  else if (xwidth && yheight && zdepth)  VarsInc.shape=1; /* box */
  else if (radius > 0 &&  yheight)       VarsInc.shape=0; /* cylinder */
  else if (radius > 0 && !yheight)       VarsInc.shape=2; /* sphere */

  if (VarsInc.shape < 0)
    exit(fprintf(stderr,"Incoherent: %s: sample has invalid dimensions.\n"
                        "ERROR       Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
  if (thickness) {
    if (radius && (radius < thickness || ( yheight && (yheight < 2*thickness)))) {
      fprintf(stderr,"Incoherent: %s: hollow sample thickness is larger than its volume (sphere/cylinder).\n"
                     "WARNING     Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP);
      thickness=0;
    }
    else if (!radius && (xwidth < 2*thickness || yheight < 2*thickness || zdepth < 2*thickness)) {
      fprintf(stderr,"Incoherent: %s: hollow sample thickness is larger than its volume (box).\n"
                     "WARNING     Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP);
      thickness=0;
    }
  }

  if (concentric && thickness<=0) {
    printf("Incoherent: %s:Can not use concentric mode\n"
           "WARNING     on non hollow shape. Ignoring.\n",
           NAME_CURRENT_COMP);
    concentric=0;
  }

  VarsInc.sigma_a= sigma_abs;
  VarsInc.sigma_i= sigma_inc;
  VarsInc.rho    = (pack/Vc);
  VarsInc.my_s   = (VarsInc.rho * 100 * VarsInc.sigma_i);
  VarsInc.my_a_v = (VarsInc.rho * 100 * VarsInc.sigma_a);

  /* now compute target coords if a component index is supplied */
  VarsInc.tx= VarsInc.ty=VarsInc.tz=0;
  if (!target_index && !target_x && !target_y && !target_z) target_index=1;
  if (target_index)
  {
    Coords ToTarget;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
    coords_get(ToTarget, &VarsInc.tx, &VarsInc.ty, &VarsInc.tz);
  }
  else
  { VarsInc.tx = target_x; VarsInc.ty = target_y; VarsInc.tz = target_z; }

  if (!(VarsInc.tx || VarsInc.ty || VarsInc.tz)) {
    MPI_MASTER(
    printf("Incoherent: %s: The target is not defined. Using direct beam (Z-axis).\n",
      NAME_CURRENT_COMP);
    );
    VarsInc.tz=1;
  }

  /* different ways of setting rectangular area */
  VarsInc.aw  = VarsInc.ah = 0;
  if (focus_xw) { VarsInc.xw = focus_xw; }
  if (focus_yh) { VarsInc.yh = focus_yh; }
  if (focus_aw) { VarsInc.aw = DEG2RAD*focus_aw; }
  if (focus_ah) { VarsInc.ah = DEG2RAD*focus_ah; }

  MPI_MASTER(
  printf("Incoherent: %s: Vc=%g [Angs] sigma_abs=%g [barn] sigma_inc=%g [barn]\n",
      NAME_CURRENT_COMP, Vc, VarsInc.sigma_a, VarsInc.sigma_i);
  );

%}

TRACE
%{
  double t0, t3;                /* Entry/exit time for outer surface */
  double t1, t2;                /* Entry/exit time for inner surface */
  double dt0, dt1, dt2, dt;     /* Flight times through sample */
  double v=0;                   /* Neutron velocity */
  double d_path;                /* Flight path length for non-scattered neutron */
  double l_i, l_o=0;            /* Flight path lenght in/out for scattered neutron */
  double my_a=0,my_t=0;         /* Velocity-dependent attenuation factor and total Xsec */
  double solid_angle=0;         /* Solid angle of target as seen from scattering point */
  double aim_x=0, aim_y=0, aim_z=1;   /* Position of target relative to scattering point */
  double v_i, v_f, E_i, E_f; /* initial and final energies and velocities */
  double dE;                 /* Energy transfer */
  int    intersect=0;
  int    flag_concentric=0;
  int    flag=0;
  double mc_trans, p_trans, mc_scatt, p_scatt, ws;
  double p_mult=1;

  #ifdef OPENACC
  #ifdef USE_OFF
  off_struct thread_offdata = offdata;
  #endif
  #else
  #define thread_offdata offdata
  #endif
  
  do { /* Main interaction loop. Ends with intersect=0 */

    /* Intersection neutron trajectory / sample (sample surface) */
    if (VarsInc.shape == 0)
      intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
    else if (VarsInc.shape == 1)
      intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
    else if (VarsInc.shape == 2)
      intersect = sphere_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius);
    #ifdef USE_OFF
    else if (VarsInc.shape == 3)
      intersect = off_intersect(&t0, &t3, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, thread_offdata );
    #endif
    
    if (intersect) {
      int flag_ishollow = 0;
      if (thickness>0) {
        if (VarsInc.shape==0 && cylinder_intersect(&t1,&t2, x,y,z,vx,vy,vz, radius-thickness,yheight-2*thickness))
          flag_ishollow=1;
        else if (VarsInc.shape==2 && sphere_intersect   (&t1,&t2, x,y,z,vx,vy,vz, radius-thickness))
          flag_ishollow=1;
        else if (VarsInc.shape==1 && box_intersect(&t1,&t2, x,y,z,vx,vy,vz, xwidth-2*thickness, yheight-2*thickness, zdepth-2*thickness))
          flag_ishollow = 1;
      }
      if (!flag_ishollow) t1 = t2 = t3; /* no empty space inside */

      dt0 = t1-t0;                /* Time in sample, ingoing */
      dt1 = t2-t1;                /* Time in hole */
      dt2 = t3-t2;                /* Time in sample, outgoing */

      if (t0 > 0) {  /* we are before the sample */
        PROP_DT(t0); /* propagates neutron to the entry of the sample */
      } else if (t1 > 0 && t1 > t0) { /* we are inside first part of the sample */
        /* no propagation, stay inside */
      } else if (t2 > 0 && t2 > t1) { /* we are in the hole */
        PROP_DT(t2); /* propagate to inner surface of 2nd part of sample */
      } else if (t3 > 0 && t3 > t2) { /* we are in the 2nd part of sample */
        /* no propagation, stay inside */
      }

      dt0=t1-(t0 > 0 ? t0 : 0); /* Time in first part of hollow/cylinder/box */
      dt1=t2-(t1 > 0 ? t1 : 0); /* Time in hole */
      dt2=t3-(t2 > 0 ? t2 : 0); /* Time in 2nd part of hollow cylinder */

      if (dt0 < 0) dt0 = 0;
      if (dt1 < 0) dt1 = 0;
      if (dt2 < 0) dt2 = 0;

      /* initialize concentric mode */
      if (concentric && !flag_concentric && t0 >= 0
       && VarsInc.shape==0 && thickness>0) {
        flag_concentric=1;
      }

      if (flag_concentric == 1) {
        dt1=dt2=0; /* force exit when reaching hole/2nd part */
      }

      if (!dt0 && !dt2) {
        intersect = 0; /* the sample was passed entirely */
        break;
      }

      p_mult = 1;
      if (!v) v = sqrt(vx*vx + vy*vy + vz*vz);
      if (v) my_a = VarsInc.my_a_v*(2200/v);
      else {
        printf("Incoherent: %s: ERROR: Null velocity\n",NAME_CURRENT_COMP);
        ABSORB; /* should never occur */
      }

      my_t = my_a + VarsInc.my_s;  /* total scattering Xsect (tmp var) */
      if (my_t <= 0) {
        printf("Incoherent: %s: ERROR: Null total cross section %g. Removing event.\n",
          NAME_CURRENT_COMP, my_t);
        ABSORB; /* should never occur */
      }
      d_path = v * (dt0 + dt2);   /* Length of full path through sample */
      /* Proba of scattering vs absorption (integrating along the whole trajectory) */
      ws = VarsInc.my_s/my_t;  /* (inc+coh)/(inc+coh+abs) */
      /* Proba of transmission along length d_path */
      p_trans = exp(-my_t*d_path);
      p_scatt = 1 - p_trans; /* portion of beam which scatters */
      flag = 0; /* flag used for propagation to exit point before ending */
      /* are we next to the exit ? probably no scattering (avoid rounding errors) */
      if (VarsInc.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) */
        if (!SCATTERED) mc_trans = 1-p_interact;
        else            mc_trans = 1-p_interact/(4*SCATTERED+1); /* reduce effect on multi scatt */
      } 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 */
        if (!mc_scatt) ABSORB;
        /* 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 */
        if (!mc_trans) ABSORB;
        p_mult *= fabs(p_trans/mc_trans);  /* attenuate beam by portion which is scattered (and left along) */
      }

      if (flag) { /* propagate to exit of sample and finish */
        intersect = 0;
        p *= p_mult; /* apply absorption correction */
        PROP_DT(dt0+dt2);
        break; /* exit main multi scatt while loop */
      }
      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
        dt = -log(1 - rand0max((1 - exp(-my_t*d_path)))) / my_t; /* length */
      l_i = dt;/* Penetration in sample: scattering+abs */
      dt /= v; /* Time from present position to scattering point */

      /* If t0 is in hole, propagate to next part of the hollow cylinder */
      if (dt1 > 0 && dt0 > 0 && dt > dt0) dt += dt1;

      PROP_DT(dt);                /* Point of scattering */

      if ((VarsInc.tx || VarsInc.ty || VarsInc.tz)) {
        aim_x = VarsInc.tx-x;       /* Vector pointing at target (anal./det.) */
        aim_y = VarsInc.ty-y;
        aim_z = VarsInc.tz-z;
      }
      if(VarsInc.aw && VarsInc.ah) {
        randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
          aim_x, aim_y, aim_z, VarsInc.aw, VarsInc.ah, ROT_A_CURRENT_COMP);
      } else if(VarsInc.xw && VarsInc.yh) {
        randvec_target_rect(&vx, &vy, &vz, &solid_angle,
          aim_x, aim_y, aim_z, VarsInc.xw, VarsInc.yh, ROT_A_CURRENT_COMP);
      } else {
        randvec_target_circle(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
      }
      NORM(vx, vy, vz);

      v_i = v;          /* Store initial velocity in case of quasielastic */
      E_i = VS2E*v_i*v_i;
      if (deltaE==0) {
	if (rand01()<f_QE)	/* Quasielastic contribution */
	  {
	    dE = gamma*tan(PI/2*randpm1());
	    E_f = E_i + dE;
	    if (E_f <= 0)
	      ABSORB;
	    v_f = SE2V*sqrt(E_f);
	    v = v_f;
	    /*          printf("vi: %g Ei: %g dE: %g Ef %g vf: %g v: %g \n",
			v_i,E_i,dE,E_f,v_f,v); */
	  }
      } else {
	E_f = E_i - Etrans + deltaE*randpm1(); // E_f=E0;
	v_f = SE2V*sqrt(E_f);
	v = v_f;
      }

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

      /* We do not consider scattering from 2nd part (outgoing) */
      p_mult *= solid_angle/4/PI;
      p *= p_mult;

      /* Polarisation part (1/3 NSF, 2/3 SF) */
      sx *= -1.0/3.0;
      sy *= -1.0/3.0;
      sz *= -1.0/3.0;

      SCATTER;

      /* test for a given multiple order */
      if (order && SCATTERED >= order) {
        intersect=0; /* reached required number of SCATTERing */
        break;       /* finish multiple scattering loop */
      }
    } /* end if intersect */
  } while (intersect); /* end do (intersect) (multiple scattering loop) */
%}

MCDISPLAY
%{
      if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {	/* OFF file */
            off_display(offdata);
      }
      else
        if (radius > 0 &&  yheight) {	/* cylinder */
          cylinder(0,0,0,radius,yheight,thickness, nx, ny, nz);
      }
      else if (xwidth && yheight) { 	/* box/rectangle */
          box(0,0,0,xwidth,yheight,zdepth,thickness, nx, ny, nz);
      }
      else if (radius > 0 && !yheight) {	/* sphere */
          sphere(0,0,0,radius);
      }
%}


END