/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Tunneling_sample
*
* %I
* Written by: Kim Lefmann
* Date: 10.05.07
* Origin: Risoe
*
* A Double-cylinder shaped all-incoherent scatterer
* with elastic, quasielastic (Lorentzian), and tunneling (sharp)
* components.
*
* %D
* A Double-cylinder shaped all-incoherent scatterer
* with both elastic, quasielastic (Lorentzian), and tunneling (sharp)
* components. No multiple scattering. Absorbtion included.
* The shape of the sample may be a box with dimensions xwidth, yheight, zdepth.
* 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.
*
* The outgoing polarization is calculated as for nuclear spin incoherence:
*   P' = 1/3*P-2/3P = -1/3P
* As above multiple scattering is ignored .
*
* Example: Tunneling_sample(thickness=0.001,radius=0.01,yheight=0.02,focus_r=0.035,
*           target_index=1)
*
* %P
* INPUT PARAMETERS:
* radius:        [m] Outer radius of sample in (x,z) plane
* yheight:       [m] vert.  dimension of sample, as a height
* thickness:     [m] Thickness of cylindrical sample in (x,z) plane
* 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
* xwidth:        [m] horiz. dimension of sample, as a width
* zdepth:        [m] depth of sample
* 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] Absorbtion cross section pr. unit cell
* sigma_inc: [barns] Total incoherent scattering cross section pr. unit cell
* Vc:         [AA^3] Unit cell volume
* p_interact:    [1] MC Probability for scattering the ray; otherwise transmit
* f_QE:          [1] Fraction of quasielastic scattering
* f_tun:         [1] Fraction of tunneling scattering (f_QE+f_tun < 1)
* gamma:       [meV] Lorentzian width of quasielastic broadening (HWHM)
* E_tun:       [meV] Tunneling energy
* target_x:      [m] X-position of target to focus at
* target_y:      [m] Y-position of target to focus at
* target_z:      [m] Z-position of target to focus at
*
* Variables calculated in the component
*
* V_my_s: [m^-1] Attenuation factor due to scattering
* V_my_a: [m^-1] Attenuation factor due to absorbtion
*
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/tests/v_sample/">Test
* results</A> (not up-to-date).
*
* %E
*******************************************************************************/

DEFINE COMPONENT Tunneling_sample

SETTING PARAMETERS (thickness=0, radius=0.01, focus_r = 0,
p_interact = 1, f_QE=0, f_tun=0, gamma=0, E_tun=0,
target_x = 0, target_y = 0, target_z = 0.235, focus_xw=0, focus_yh=0,
focus_aw=0, focus_ah=0, xwidth=0, yheight=0.05, zdepth=0, sigma_abs=5.08, sigma_inc=4.935, Vc=13.827, int target_index=0)

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

SHARE
%{
  struct StructVarsV {
    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;
    char isrect;       /* true when sample is a box */
    double distance;   /* when non zero, gives rect target distance */
    double aw, ah;     /* rectangular angular dimensions */
    double xw, yh;     /* rectangular metrical dimensions */
    double tx, ty, tz; /* target coords */
  };
%}

DECLARE
%{
  struct StructVarsV VarsV;
  double ftun;
  double fQE;
%}

INITIALIZE
%{
  if (!xwidth || !yheight || !zdepth) /* Cannot define a rectangle */
    if (!radius || !yheight)          /* Cannot define a cylinder either */
      exit (fprintf (stderr, "V_sample: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
    else /* It is a cylinder */
      VarsV.isrect = 0;
  else /* It is a rectangle */
    VarsV.isrect = 1;

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

  /* now compute target coords if a component index is supplied */
  VarsV.tx = VarsV.ty = VarsV.tz = 0;
  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, &VarsV.tx, &VarsV.ty, &VarsV.tz);
  } else {
    VarsV.tx = target_x;
    VarsV.ty = target_y;
    VarsV.tz = target_z;
  }

  if (!(VarsV.tx || VarsV.ty || VarsV.tz))
    printf ("Tunneling_sample: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP);

  VarsV.distance = sqrt (VarsV.tx * VarsV.tx + VarsV.ty * VarsV.ty + VarsV.tz * VarsV.tz);

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

  /* Check that probabilities are positive and do not exceed unity */
  if (f_tun < 0)
    ftun = 0;
  else
    ftun = f_tun;
  if (f_QE < 0)
    fQE = 0;
  else
    fQE = f_QE;
  if ((ftun + fQE) > 1) {
    ftun = 0;
    printf ("Tunneling_sample: Sum of inelastic probabilities > 1. Setting f_tun=0");
    if (fQE > 1) {
      fQE = 0;
      printf ("Tunneling_sample: Probability fQE > 1. Setting fQE=0.");
    }
  }
%}

TRACE
%{
  double t0, t3;                          /* Entry/exit time for outer cylinder */
  double t1, t2;                          /* Entry/exit time for inner cylinder */
  double v;                               /* Neutron velocity */
  double dt0, dt1, dt2, dt;               /* Flight times through sample */
  double l_full;                          /* 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;                        /* Velocity-dependent attenuation factor */
  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 */
  double scatt_choice;                    /* Representing random choice of scattering type */
  int intersect = 0;

  if (VarsV.isrect)
    intersect = box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
  else
    intersect = cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
  if (intersect) {
    if (t0 < 0)
      ABSORB; /* we already passed the sample; this is illegal */
    /* Neutron enters at t=t0. */
    if (VarsV.isrect)
      t1 = t2 = t3;
    else if (!thickness || !cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius - thickness, yheight))
      t1 = t2 = t3;

    dt0 = t1 - t0; /* Time in sample, ingoing */
    dt1 = t2 - t1; /* Time in hole */
    dt2 = t3 - t2; /* Time in sample, outgoing */
    v = sqrt (vx * vx + vy * vy + vz * vz);
    l_full = v * (dt0 + dt2); /* Length of full path through sample */
    if (v)
      my_a = VarsV.my_a_v * (2200 / v);

    if (p_interact >= 1 || rand01 () < p_interact) /* Scattering */
    {
      dt = rand01 () * (dt0 + dt2); /* Time of scattering (relative to t0) */
      l_i = v * dt;                 /* Penetration in sample: scattering+abs */
      if (dt > dt0)
        dt += dt1; /* jump to 2nd side of cylinder */

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

      if ((VarsV.tx || VarsV.ty || VarsV.tz)) {
        aim_x = VarsV.tx - x; /* Vector pointing at target (anal./det.) */
        aim_y = VarsV.ty - y;
        aim_z = VarsV.tz - z;
      }
      if (VarsV.aw && VarsV.ah) {
        randvec_target_rect_angular (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, VarsV.aw, VarsV.ah, ROT_A_CURRENT_COMP);
      } else if (VarsV.xw && VarsV.yh) {
        randvec_target_rect (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, VarsV.xw, VarsV.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);

      scatt_choice = rand01 (); /* chooses type of scattering */
      v_i = v;                  /* Store initial velocity in case of inel. */
      E_i = VS2E * v_i * v_i;
      if (scatt_choice < (fQE + ftun)) /* Inelastic choices */
      {
        if (scatt_choice < fQE) /* Quasielastic */
        {
          dE = gamma * tan (PI / 2 * randpm1 ());
        } else {
          if (randpm1 () > 0)
            dE = E_tun;
          else
            dE = -E_tun;
        }
        E_f = E_i + dE;
        if (E_f <= 0)
          ABSORB;
        v_f = SE2V * sqrt (E_f);
        v = v_f;
      }

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

      if (!VarsV.isrect) {
        if (!cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight)) {
          /* ??? did not hit cylinder */
          printf ("FATAL ERROR: Did not hit cylinder from inside.\n");
          exit (1);
        }
        dt = t3; /* outgoing point */
        if (thickness && cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius - thickness, yheight) && t2 > 0)
          dt -= (t2 - t1); /* Subtract hollow part */
      } else {
        if (!box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth)) {
          /* ??? did not hit box */
          printf ("FATAL ERROR: Did not hit box from inside.\n");
          exit (1);
        }
        dt = t3;
      }
      l_o = v * dt; /* trajectory after scattering point: absorption only */

      p *= v / v_i * l_full * VarsV.my_s * exp (-my_a * (l_i + v_i / v * l_o) - VarsV.my_s * l_i);
      /* We do not consider scattering from 2nd part (outgoing) */
      p /= 4 * PI / solid_angle;
      p /= p_interact;

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

      SCATTER;
    } else /* Transmitting; always elastic */
    {
      p *= exp (-(my_a + VarsV.my_s) * l_full);
      p /= (1 - p_interact);
    }
  }
%}

MCDISPLAY
%{

  if (!VarsV.isrect) {
    circle ("xz", 0, yheight / 2.0, 0, radius);
    circle ("xz", 0, -yheight / 2.0, 0, radius);
    line (-radius, -yheight / 2.0, 0, -radius, +yheight / 2.0, 0);
    line (+radius, -yheight / 2.0, 0, +radius, +yheight / 2.0, 0);
    line (0, -yheight / 2.0, -radius, 0, +yheight / 2.0, -radius);
    line (0, -yheight / 2.0, +radius, 0, +yheight / 2.0, +radius);
    if (thickness) {
      double radius_i = radius - thickness;
      circle ("xz", 0, yheight / 2.0, 0, radius_i);
      circle ("xz", 0, -yheight / 2.0, 0, radius_i);
      line (-radius_i, -yheight / 2.0, 0, -radius_i, +yheight / 2.0, 0);
      line (+radius_i, -yheight / 2.0, 0, +radius_i, +yheight / 2.0, 0);
      line (0, -yheight / 2.0, -radius_i, 0, +yheight / 2.0, -radius_i);
      line (0, -yheight / 2.0, +radius_i, 0, +yheight / 2.0, +radius_i);
    }
  } else {
    double xmin = -0.5 * xwidth;
    double xmax = 0.5 * xwidth;
    double ymin = -0.5 * yheight;
    double ymax = 0.5 * yheight;
    double zmin = -0.5 * zdepth;
    double zmax = 0.5 * zdepth;
    multiline (5, xmin, ymin, zmin, xmax, ymin, zmin, xmax, ymax, zmin, xmin, ymax, zmin, xmin, ymin, zmin);
    multiline (5, xmin, ymin, zmax, xmax, ymin, zmax, xmax, ymax, zmax, xmin, ymax, zmax, xmin, ymin, zmax);
    line (xmin, ymin, zmin, xmin, ymin, zmax);
    line (xmax, ymin, zmin, xmax, ymin, zmax);
    line (xmin, ymax, zmin, xmin, ymax, zmax);
    line (xmax, ymax, zmin, xmax, ymax, zmax);
  }
%}

END