/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Single_crystal
*
* %I
* Written by: Kristian Nielsen
* Date: December 1999
* Version: $Revision: 1.36 $
* Origin: Risoe
* Release: McStas 1.6
* Modified by: EF, 22nd Apr 2003 : now uses Read_Table library
*
* Mosaic single crystal with multiple scattering vectors.
*
* %D
* Single crystal with mosaic. Delta-D/D option for finite-size effects.
* Rectangular geometry. Multiple scattering and secondary extinction included.
* The mosaic may EITHER be specified isotropic by setting the mosaic input
* parameter, OR anisotropic by setting the mosaic_h, mosaic_v, and mosaic_n
* parameters. Geometry is a box.
* Crystal structure is specified with an ascii data file. Each line contains
* 4 or more numbers, separated by white spaces:
*
*       h k l ... F2
*
* The first three numbers are the (h,k,l) indices of the reciprocal lattice
* point, and the last number is the value of the structure factor |F|**2, in
* barns. The rest of the numbers are not used; the file is in the format
* output by the Crystallographica program.
* Lines begining by '#' are read as comments (ignored).
*
* Example: Single_crystal(xwidth=0.01, yheight=0.01, zthick=0.01,
*           delta_d_d=1e-4, mosaic = 5,
*           ax = 3.8186, ay = 0, az = 0,
*           bx = 0, by = 3.8843, bz = 0,
*           cx = 0, cy = 0, cz = 11.6777,
*           reflections="YBaCuO.lau")
*
* Vanadium incoherent elastic scattering with multiple scattering
* Single_crystal(xwidth=0.01, yheight=0.01, zthick=0.01,
*           reflections="", absorbtion=5.08, incoherent=4.935,
*           ax=3.0282, by=3.0282, cz=3.0282/2)
* %P
* INPUT PARAMETERS:
*
* xwidth:    Width of crystal [m]
* yheight:   Height of crystal [m]
* zthick:    Thichness of crystal (no extinction simulated) [m]
* delta_d_d: Lattice spacing variance, gaussian RMS [1]
* mosaic:    Crystal mosaic (isotropic), gaussian RMS [arc minutes]
* mosaic_h:  Horizontal (rotation around Y) mosaic (anisotropic),
*               gaussian RMS [arc minutes]
* mosaic_v:  Vertical (rotation around Z) mosaic (anisotropic),
*               gaussian RMS [arc minutes]
* mosaic_n:  Out-of-plane (Rotation around X) mosaic (anisotropic),
*               gaussian RMS [arc minutes]
* ax:        [AA]
* ay:        Coordinates of first unit cell vector [AA]
* az:        [AA]
* bx:        [AA]
* by:        Coordinates of second unit cell vector [AA]
* bz:        [AA]
* cx:        [AA]
* cy:        Coordinates of third unit cell vector [AA]
* cz:        [AA]
* reflections: File name containing structure factors of reflections. Use
*              empty ("") or NULL for incoherent scattering only [string]
* order:      limit multiple scattering up to given order
*                 (0: all, 1: first, 2: second, ...) [1]
*
* Optional input parameters:
*
* p_transmit: Monte Carlo probability for neutrons to be transmitted
*               without any scattering. Used to improve statistics from
*               weak reflections [1]
* absorbtion: Absorbtion cross-section per unit cell at 2200 m/s [barns]
* incoherent: Incoherent scattering cross-section per unit cell [barns]
* aa:         [deg].
* bb:         unit cell angles alpha, beta and gamma. Then uses norms of
*               vectors a,b and c as lattice parameters [deg]
* cc:         [deg].
*
* OUTPUT PARAMETERS:
*
* hkl_info: Internal [structure]
*
* %L
* See <a href="icsd.ill.fr">ICSD</a> powder diffraction data base
* %L
* Cross sections for single elements: http://www.ncnr.nist.gov/resources/n-lengths/
* %L
* Cross sections for compounds:       http://www.ncnr.nist.gov/resources/sldcalc.html
*
* %E
****************************************************************************/

/*
  %D
  Overview of algorithm:

  (1). The neutron intersects the crystal at (x,y,z) with given
       incoming wavevector ki=(kix,kiy,kiz).

  (2). Every reciprocal lattice point tau of magnitude less than 2*ki
       is considered for scattering. The scattering probability is the
       area of the intersection of the Ewald sphere (approximated by
       the tangential plane) with the 3-D Gaussian mosaic of the point
       tau.

  (3). The total coherent scattering cross section is computed as the
       sum over all tau. Together with the absorption and incoherent
       scattering cross section and known potential flight-length
       l_full through the sample, we can compute the probability of
       the four events absorption, coherent scattering, incoherent
       scattering, and transmission.

  (4). Absorption is never simulated explicitly, just incorporated in
       the neutron weight.

  (5). Transmission in the first event is selected with the Monte
       Carlo probability p_transmit, which defaults to the actual
       transmission probability. After the first event, transmission
       is selected with the correct Monte Carlo probability.

  (6). Incoherent scattering is done simply by selecting a random
       direction for the outgoing wave vector kf.

  (7). For coherent scattering, a reciprocal lattice point is selected
       using the relative probabilities computed in (2), and the
       weight is adjusted with the contribution from the structure
       factors (this way all reflections will get equally good
       statistics in the detector).

  (8). The outgoing wave vector direction is picked at random using
       the intersecting 2-D Gauss computed in (2). The vector is
       normalized to the length of ki (elastic scattering) to account
       for the error caused by the planar approximation of the Ewald
       sphere.

  (9). The process is repeated from (2) with kf as new initial wave
       vector ki.

  %E
 */

DEFINE COMPONENT Single_crystal
DEFINITION PARAMETERS(reflections)
SETTING PARAMETERS(xwidth, yheight, zthick, delta_d_d=1e-4,
            mosaic = -1, mosaic_h = -1, mosaic_v = -1, mosaic_n = -1,
            ax = 3.8186, ay = 0, az = 0,
            bx = 0, by = 3.8843, bz = 0,
            cx = 0, cy = 0, cz = 11.6777,
            p_transmit = -1, absorbtion = 0, incoherent = 0,
            aa=0, bb=0, cc=0, order=0)
OUTPUT PARAMETERS(hkl_info)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

SHARE
%{
  /* used for reading data table from file */
  %include "read_table-lib"
/* Declare structures and functions only once in each instrument. */
#ifndef SINGLE_CRYSTAL_DECL
#define SINGLE_CRYSTAL_DECL

    struct hkl_data
    {
      int h,k,l;                  /* Indices for this reflection */
      double F2;                  /* Value of structure factor */
      double tau_x, tau_y, tau_z; /* Coordinates in reciprocal space */
      double tau;                 /* Length of (tau_x, tau_y, tau_z) */
      double u1x, u1y, u1z;       /* First axis of local coordinate system */
      double u2x, u2y, u2z;       /* Second axis of local coordinate system */
      double u3x, u3y, u3z;       /* Third axis of local coordinate system */
      double sig1, sig2, sig3;    /* RMSs of Gauss axis */
      double sig123;              /* The product sig1*sig2*sig3 */
      double m1, m2, m3;          /* Diagonal matrix representation of Gauss */
      double cutoff;              /* Cutoff value for Gaussian tails */
    };

  struct hkl_info_struct
    {
      struct hkl_data *list;      /* Reflection array */
      int count;                  /* Number of reflections */
      struct tau_data *tau_list;  /* Reflections close to Ewald Sphere */
      double m_delta_d_d;         /* Delta-d/d FWHM */
      double m_ax,m_ay,m_az;      /* First unit cell axis (direct space, AA) */
      double m_bx,m_by,m_bz;      /* Second unit cell axis */
      double m_cx,m_cy,m_cz;      /* Third unit cell axis */
      double asx,asy,asz;         /* First reciprocal lattice axis (1/AA) */
      double bsx,bsy,bsz;         /* Second reciprocal lattice axis */
      double csx,csy,csz;         /* Third reciprocal lattice axis */
      double V0;                  /* Unit cell volume (AA**3) */
    };

  struct tau_data
    {
      int index;                  /* Index into reflection table */
      double refl;
      double xsect;
      double sigma_1, sigma_2;
      /* The following vectors are in local koordinates. */
      double kix, kiy, kiz;       /* Initial wave vector */
      double rho_x, rho_y, rho_z; /* The vector ki - tau */
      double rho;                 /* Length of rho vector */
      double ox, oy, oz;          /* Origin of Ewald sphere tangent plane */
      double nx, ny, nz;          /* Normal vector of Ewald sphere tangent */
      double b1x, b1y, b1z;       /* Spanning vectors of Ewald sphere tangent */
      double b2x, b2y, b2z;
      double l11, l12, l22;       /* Cholesky decomposition L of 2D Gauss */
      double det_L;               /* Determinant of L */
      double y0x, y0y;            /* 2D Gauss center in tangent plane */
    };

  int
  read_hkl_data(char *SC_file, struct hkl_info_struct *info,
    double SC_mosaic, double SC_mosaic_h, double SC_mosaic_v, double SC_mosaic_n)
  {
    struct hkl_data *list = NULL;
    int size = 0;
    t_Table sTable; /* sample data table structure from SC_file */
    int i=0;

    if (!SC_file || !strlen(SC_file)) {
      info->count = 0;
      return(-1);
    }

    Table_Read(&sTable, SC_file, 1); /* read 1st block data from SC_file into sTable*/
    if (sTable.columns < 4) {
      fprintf(stderr, "Single_crystal: Error: The number of columns in %s should be at least %d for [h,k,l,F2]\n", SC_file, 4);
      return(0);
    }
    if (!sTable.rows) {
      fprintf(stderr, "Single_crystal: Error: The number of rows in %s should be at least %d\n", SC_file, 1);
      return(0);
    } else size = sTable.rows;

    /* allocate hkl_data array */
    list = (struct hkl_data*)malloc(size*sizeof(struct hkl_data));

    for (i=0; i<size; i++)
    {
      double h, k, l, F2;
      double b1[3], b2[3];

      /* get data from table */
      h = Table_Index(sTable, i, 0);
      k = Table_Index(sTable, i, 1);
      l = Table_Index(sTable, i, 2);
      F2= Table_Index(sTable, i, sTable.columns-1); /* last column */

      list[i].h = h;
      list[i].k = k;
      list[i].l = l;
      list[i].F2 = F2;
      /* Precompute some values */
      list[i].tau_x = h*info->asx + k*info->bsx + l*info->csx;
      list[i].tau_y = h*info->asy + k*info->bsy + l*info->csy;
      list[i].tau_z = h*info->asz + k*info->bsz + l*info->csz;
      list[i].tau = sqrt(list[i].tau_x*list[i].tau_x +
                         list[i].tau_y*list[i].tau_y +
                         list[i].tau_z*list[i].tau_z);
      list[i].u1x = list[i].tau_x/list[i].tau;
      list[i].u1y = list[i].tau_y/list[i].tau;
      list[i].u1z = list[i].tau_z/list[i].tau;
      list[i].sig1 = FWHM2RMS*info->m_delta_d_d*list[i].tau;
      /* Find two arbitrary axes perpendicular to tau and each other. */
      normal_vec(&b1[0], &b1[1], &b1[2],
                 list[i].u1x, list[i].u1y, list[i].u1z);
      vec_prod(b2[0], b2[1], b2[2],
               list[i].u1x, list[i].u1y, list[i].u1z,
               b1[0], b1[1], b1[2]);
      /* Find the two mosaic axes perpendicular to tau. */
      if(SC_mosaic > 0) {
        /* Use isotropic mosaic. */
        list[i].u2x = b1[0];
        list[i].u2y = b1[1];
        list[i].u2z = b1[2];
        list[i].sig2 = FWHM2RMS*list[i].tau*MIN2RAD*SC_mosaic;
        list[i].u3x = b2[0];
        list[i].u3y = b2[1];
        list[i].u3z = b2[2];
        list[i].sig3 = FWHM2RMS*list[i].tau*MIN2RAD*SC_mosaic;
      } else if(SC_mosaic_h > 0 && SC_mosaic_v > 0 && SC_mosaic_n > 0) {
        /* Use anisotropic mosaic. */
        double mos_D[3][3];
        double mos_a, mos_b, mos_c, mos_det, mos_l1, mos_l2;
        double mos_e1_x, mos_e1_y, mos_e1_len, mos_e2_x, mos_e2_y, mos_e2_len;
        int j1, j2;
        /* Define the 3D Gaussian for the mosaic in Euler coordinates. */
        mos_D[0][0] = 1/((FWHM2RMS*MIN2RAD*SC_mosaic_n)*(FWHM2RMS*MIN2RAD*SC_mosaic_n));
        mos_D[0][1] = 0;
        mos_D[0][2] = 0;
        mos_D[1][0] = 0;
        mos_D[1][1] = 1/((FWHM2RMS*MIN2RAD*SC_mosaic_h)*(FWHM2RMS*MIN2RAD*SC_mosaic_h));
        mos_D[1][2] = 0;
        mos_D[2][0] = 0;
        mos_D[2][1] = 0;
        mos_D[2][2] = 1/((FWHM2RMS*MIN2RAD*SC_mosaic_v)*(FWHM2RMS*MIN2RAD*SC_mosaic_v));
        /* ToDo: This is the place to apply a rotation U to the mosaic Gaussian,
           by setting D := U^T D U. */
        /* Now compute the 2D intersection Gauss (B^T D B) between the 3D SC_mosaic
           Gauss D and the plane spanned by b1 and b2 (matrix B). The matrix
           (B^T D B) is 2x2 symmetric and positive definite.*/
        mos_a = mos_b = mos_c = 0;
        for(j1 = 0; j1 < 3; j1++) {
          for(j2 = 0; j2 < 3; j2++) {
            mos_a += b1[j1]*b1[j2]*mos_D[j1][j2];
            mos_b += b1[j1]*b2[j2]*mos_D[j1][j2];
            mos_c += b2[j1]*b2[j2]*mos_D[j1][j2];
          }
        }
        /* Compute the eigenvalues of the 2D Gaussian. */
        if(mos_b*mos_b < 1e-20*(mos_a*mos_a + mos_c*mos_c)) {
          /* Matrix is diagonal. */
          mos_l1 = mos_a;
          mos_e1_x = 1;
          mos_e1_y = 0;
          mos_l2 = mos_c;
          mos_e2_x = 0;
          mos_e2_y = 1;
        } else {                  /* Not diagonal */
          mos_det = (mos_a - mos_c)*(mos_a - mos_c) + 4*mos_b*mos_b;
          mos_l1 = 0.5*(mos_a + mos_c - sqrt(mos_det));
          mos_l2 = 0.5*(mos_a + mos_c + sqrt(mos_det));
          /* Compute the eigenvectors. */
          mos_e1_x = -mos_b;
          mos_e1_y = mos_a - mos_l1;
          mos_e1_len = sqrt(mos_e1_x*mos_e1_x + mos_e1_y*mos_e1_y);
          mos_e1_x /= mos_e1_len;
          mos_e1_y /= mos_e1_len;
          mos_e2_x = mos_c - mos_l2;
          mos_e2_y = -mos_b;
          mos_e2_len = sqrt(mos_e2_x*mos_e2_x + mos_e2_y*mos_e2_y);
          mos_e2_x /= mos_e2_len;
          mos_e2_y /= mos_e2_len;
        }
        /* Now finally compute the axes (in reciprocal space) and RMSs of the
           reciprocal lattice point Gaussian that are perpendicular to tau. */
        list[i].u2x = b1[0]*mos_e2_x + b2[0]*mos_e2_y;
        list[i].u2y = b1[1]*mos_e2_x + b2[1]*mos_e2_y;
        list[i].u2z = b1[2]*mos_e2_x + b2[2]*mos_e2_y;
        list[i].sig2 = list[i].tau/sqrt(mos_l1);
        list[i].u3x = b1[0]*mos_e1_x + b2[0]*mos_e1_y;
        list[i].u3y = b1[1]*mos_e1_x + b2[1]*mos_e1_y;
        list[i].u3z = b1[2]*mos_e1_x + b2[2]*mos_e1_y;
        list[i].sig3 = list[i].tau/sqrt(mos_l2);
      } else {
        fprintf(stderr,
                "Single_crystal: Error: EITHER mosaic OR (mosaic_h, mosaic_v, mosaic_n)\n"
                "must be given and be >0.\n");
        return(0);
      }
      list[i].sig123 = list[i].sig1*list[i].sig2*list[i].sig3;
      list[i].m1 = 1/(2*list[i].sig1*list[i].sig1);
      list[i].m2 = 1/(2*list[i].sig2*list[i].sig2);
      list[i].m3 = 1/(2*list[i].sig3*list[i].sig3);
      /* Set Gauss cutoff to 5 times the maximal sigma. */
      if(list[i].sig1 > list[i].sig2)
        if(list[i].sig1 > list[i].sig3)
          list[i].cutoff = 5*list[i].sig1;
        else
          list[i].cutoff = 5*list[i].sig3;
      else
        if(list[i].sig2 > list[i].sig3)
          list[i].cutoff = 5*list[i].sig2;
        else
          list[i].cutoff = 5*list[i].sig3;
    }
    Table_Free(&sTable);
    info->list = list;
    info->count = i;
    info->tau_list = malloc(i*sizeof(*info->tau_list));
    if(!info->tau_list)
    {
      fprintf(stderr, "Single_crystal: Error: Out of memory!\n");
      return(0);
    }
    return(info->count = i);
  }
#endif /* !SINGLE_CRYSTAL_DECL */

%}

DECLARE
%{
  struct hkl_info_struct hkl_info;
%}

INITIALIZE
%{
  double tmp_x, tmp_y, tmp_z;
  double as, bs, cs;

  if (aa*bb*cc != 0)
  {
    as = sqrt(ax*ax+ay*ay+az*az);
    bs = sqrt(bx*bx+by*by+bz*bz);
    cs = sqrt(cx*cx+cy*cy+cz*cz);
    ax = as; ay = 0; az = 0;
    bx = bs*cos(cc*DEG2RAD);
    by = bs*sin(cc*DEG2RAD);
    bz = 0;
    cx = cs*cos(bb*DEG2RAD);
    cy = cs*(cos(aa*DEG2RAD)-cos(cc*DEG2RAD)*cos(bb*DEG2RAD))/sin(cc*DEG2RAD);
    cz = sqrt(cs*cs - cx*cx - cy*cy);
  }

  hkl_info.m_delta_d_d = delta_d_d;
  hkl_info.m_ax = ax;
  hkl_info.m_ay = ay;
  hkl_info.m_az = az;
  hkl_info.m_bx = bx;
  hkl_info.m_by = by;
  hkl_info.m_bz = bz;
  hkl_info.m_cx = cx;
  hkl_info.m_cy = cy;
  hkl_info.m_cz = cz;
  /* Compute reciprocal lattice vectors. */
  vec_prod(tmp_x, tmp_y, tmp_z, bx, by, bz, cx, cy, cz);
  hkl_info.V0 = fabs(scalar_prod(ax, ay, az, tmp_x, tmp_y, tmp_z));
  hkl_info.asx = 2*PI/hkl_info.V0*tmp_x;
  hkl_info.asy = 2*PI/hkl_info.V0*tmp_y;
  hkl_info.asz = 2*PI/hkl_info.V0*tmp_z;
  vec_prod(tmp_x, tmp_y, tmp_z, cx, cy, cz, ax, ay, az);
  hkl_info.bsx = 2*PI/hkl_info.V0*tmp_x;
  hkl_info.bsy = 2*PI/hkl_info.V0*tmp_y;
  hkl_info.bsz = 2*PI/hkl_info.V0*tmp_z;
  vec_prod(tmp_x, tmp_y, tmp_z, ax, ay, az, bx, by, bz);
  hkl_info.csx = 2*PI/hkl_info.V0*tmp_x;
  hkl_info.csy = 2*PI/hkl_info.V0*tmp_y;
  hkl_info.csz = 2*PI/hkl_info.V0*tmp_z;
  /* Read in structure factors, and do some pre-calculations. */
  if (!read_hkl_data(reflections, &hkl_info, mosaic, mosaic_h, mosaic_v, mosaic_n))
    exit(0);
  if (hkl_info.count)
    printf("Single_crystal: %s: Read %d reflections from file '%s'\n",
      NAME_CURRENT_COMP, hkl_info.count, reflections);
  else printf("Single_crystal: %s: Using incoherent elastic scattering only\n",
      NAME_CURRENT_COMP);
%}

TRACE
%{
  double t1, t2;                /* Entry and exit times in sample */
  struct hkl_data *L;           /* Structure factor list */
  int i;                        /* Index into structure factor list */
  struct tau_data *T;           /* List of reflections close to Ewald sphere */
  int j;                        /* Index into reflection list */
  int event_counter;            /* scattering event counter */
  double kix, kiy, kiz, ki;     /* Initial wave vector [1/AA] */
  double kfx, kfy, kfz;         /* Final wave vector */
  double v;                     /* Neutron velocity */
  double tau_max;               /* Max tau allowing reflection at this ki */
  double rho_x, rho_y, rho_z;   /* the vector ki - tau */
  double rho;
  double diff;                  /* Deviation from Bragg condition */
  double ox, oy, oz;            /* Origin of Ewald sphere tangent plane */
  double b1x, b1y, b1z;         /* First vector spanning tangent plane */
  double b2x, b2y, b2z;         /* Second vector spanning tangent plane */
  double n11, n12, n22;         /* 2D Gauss description matrix N */
  double det_N;                 /* Determinant of N */
  double inv_n11, inv_n12, inv_n22; /* Inverse of N */
  double l11, l12, l22;         /* Cholesky decomposition L of 1/2*inv(N) */
  double det_L;                 /* Determinant of L */
  double Bt_D_O_x, Bt_D_O_y;    /* Temporaries */
  double y0x, y0y;              /* Center of 2D Gauss in plane coordinates */
  double alpha;                 /* Offset of 2D Gauss center from 3D center */
  int tau_count;                /* Number of reflections within cutoff */
  double V0;                    /* Volume of unit cell */
  double l_full;                /* Neutron path length for transmission */
  double l;                     /* Path length to scattering event */
  double abs_xsect, abs_xlen;   /* Absorbtion cross section and length */
  double inc_xsect, inc_xlen;   /* Incoherent cross section and length */
  double coh_xsect, coh_xlen;   /* Coherent cross section and length */
  double tot_xsect, tot_xlen;   /* Total cross section and length */
  double z1, z2, y1, y2;        /* Temporaries to choose kf from 2D Gauss */
  double adjust, coh_refl;      /* Temporaries */
  double r, sum;                /* Temporaries */
  double xsect_factor;          /* Common factor in coherent cross-section */
  double p_trans;               /* Transmission probability */
  double mc_trans, mc_interact; /* Transmission, interaction MC choices */

  if(box_intersect(&t1, &t2, x, y, z, vx, vy, vz,
                   xwidth, yheight, zthick) && t2 > 0)
  {                             /* Neutron intersects crystal */
    if(t1 > 0)
      PROP_DT(t1);                /* Move to crystal surface if not inside */
    v  = sqrt(vx*vx + vy*vy + vz*vz);
    ki = V2K*v;
    event_counter = 0;
    abs_xsect = absorbtion*2200/v;
    inc_xsect = incoherent;
    abs_xlen  = 1e2*abs_xsect/V0;
    inc_xlen  = 1e2*inc_xsect/V0;
    L = hkl_info.list;
    T = hkl_info.tau_list;
    V0= hkl_info.V0;
    for(;;)                       /* Loop over multiple scattering events */
    {
      /* Compute intersection between neutron flight path and sample. */
      if(!box_intersect(&t1, &t2, x, y, z, vx, vy, vz,
                        xwidth, yheight, zthick) || t2*v < -1e-9 || t1*v > 1e-9)
      {
        /* neutron is leaving the sample */
        fprintf(stderr,
                "Single_crystal: %s: Warning: neutron has unexpectedly left the crystal!\n"
                "                t1=%g t2=%g x=%g y=%g z=%g vx=%g vy=%g vz=%g\n",
                NAME_CURRENT_COMP, t1, t2, x, y, z, vx, vy, vz);
        break;
      }
      l_full = t2*v;

      /* (1). Compute incoming wave vector ki */
      kix = V2K*vx;
      kiy = V2K*vy;
      kiz = V2K*vz;

      /* (2). Intersection of Ewald sphere with recipprocal lattice points */

      /* Max possible tau with 5*sigma delta-d/d cutoff. */
      tau_max = 2*ki/(1 - 5*delta_d_d);

      coh_xsect = 0;
      coh_refl  = 0;
      xsect_factor = pow(2*PI, 5.0/2.0)/(V0*ki*ki);
      for(i = j = 0; i < hkl_info.count; i++)
      {
        /* Assuming reflections are sorted, stop search when max tau exceeded. */
        if(L[i].tau > tau_max)
          break;
        /* Check if this reciprocal lattice point is close enough to the
           Ewald sphere to make scattering possible. */
        rho_x = kix - L[i].tau_x;
        rho_y = kiy - L[i].tau_y;
        rho_z = kiz - L[i].tau_z;
        rho = sqrt(rho_x*rho_x + rho_y*rho_y + rho_z*rho_z);
        diff = fabs(rho - ki);
        /* Check if scattering is possible (cutoff of Gaussian tails). */
        if(diff <= L[i].cutoff)
        {
          /* Store reflection. */
          T[j].index = i;
          /* Get ki vector in local coordinates. */
          T[j].kix = kix*L[i].u1x + kiy*L[i].u1y + kiz*L[i].u1z;
          T[j].kiy = kix*L[i].u2x + kiy*L[i].u2y + kiz*L[i].u2z;
          T[j].kiz = kix*L[i].u3x + kiy*L[i].u3y + kiz*L[i].u3z;
          T[j].rho_x = T[j].kix - L[i].tau;
          T[j].rho_y = T[j].kiy;
          T[j].rho_z = T[j].kiz;
          T[j].rho = rho;
          /* Compute the tangent plane of the Ewald sphere. */
          T[j].nx = T[j].rho_x/T[j].rho;
          T[j].ny = T[j].rho_y/T[j].rho;
          T[j].nz = T[j].rho_z/T[j].rho;
          ox = (ki - T[j].rho)*T[j].nx;
          oy = (ki - T[j].rho)*T[j].ny;
          oz = (ki - T[j].rho)*T[j].nz;
          T[j].ox = ox;
          T[j].oy = oy;
          T[j].oz = oz;
          /* Compute unit vectors b1 and b2 that span the tangent plane. */
          normal_vec(&b1x, &b1y, &b1z, T[j].nx, T[j].ny, T[j].nz);
          vec_prod(b2x, b2y, b2z, T[j].nx, T[j].ny, T[j].nz, b1x, b1y, b1z);
          T[j].b1x = b1x;
          T[j].b1y = b1y;
          T[j].b1z = b1z;
          T[j].b2x = b2x;
          T[j].b2y = b2y;
          T[j].b2z = b2z;
          /* Compute the 2D projection of the 3D Gauss of the reflection. */
          /* The symmetric 2x2 matrix N describing the 2D gauss. */
          n11 = L[i].m1*b1x*b1x + L[i].m2*b1y*b1y + L[i].m3*b1z*b1z;
          n12 = L[i].m1*b1x*b2x + L[i].m2*b1y*b2y + L[i].m3*b1z*b2z;
          n22 = L[i].m1*b2x*b2x + L[i].m2*b2y*b2y + L[i].m3*b2z*b2z;
          /* The (symmetric) inverse matrix of N. */
          det_N = n11*n22 - n12*n12;
          inv_n11 = n22/det_N;
          inv_n12 = -n12/det_N;
          inv_n22 = n11/det_N;
          /* The Cholesky decomposition of 1/2*inv_n (lower triangular L). */
          l11 = sqrt(inv_n11/2);
          l12 = inv_n12/(2*l11);
          l22 = sqrt(inv_n22/2 - l12*l12);
          T[j].l11 = l11;
          T[j].l12 = l12;
          T[j].l22 = l22;
          det_L = l11*l22;
          T[j].det_L = det_L;
          /* The product B^T D o. */
          Bt_D_O_x = b1x*L[i].m1*ox + b1y*L[i].m2*oy + b1z*L[i].m3*oz;
          Bt_D_O_y = b2x*L[i].m1*ox + b2y*L[i].m2*oy + b2z*L[i].m3*oz;
          /* Center of 2D Gauss in plane coordinates. */
          y0x = -(Bt_D_O_x*inv_n11 + Bt_D_O_y*inv_n12);
          y0y = -(Bt_D_O_x*inv_n12 + Bt_D_O_y*inv_n22);
          T[j].y0x = y0x;
          T[j].y0y = y0y;
          /* Factor alpha for the distance of the 2D Gauss from the origin. */
          alpha = L[i].m1*ox*ox + L[i].m2*oy*oy + L[i].m3*oz*oz -
                       (y0x*y0x*n11 + y0y*y0y*n22 + 2*y0x*y0y*n12);
          T[j].refl = xsect_factor*det_L*exp(-alpha)/L[i].sig123;
          coh_refl += T[j].refl;
          T[j].xsect = T[j].refl*L[i].F2;
          coh_xsect += T[j].xsect;
          j++;
        }
      } /* end for */
      tau_count = j;

      /* (3). Probabilities of the different possible interactions. */
      tot_xsect = abs_xsect + inc_xsect + coh_xsect;
      /* Cross-sections are in barns = 10**-28 m**2, and unit cell volumes are
         in AA**3 = 10**-30 m**2. Hence a factor of 100 is used to convert
         scattering lengths to m**-1 */
      coh_xlen = 1e2*coh_xsect/V0;
      tot_xlen = 1e2*tot_xsect/V0;

      /* (5). Transmission */
      p_trans = exp(-tot_xlen*l_full);
      if(!event_counter && p_transmit >= 0 && p_transmit <= 1) {
        mc_trans = p_transmit; /* first event */
      } else {
        mc_trans = p_trans;
      }
      if(mc_trans >= 1 || rand01() < mc_trans)  /* Transmit */
      {
        p *= p_trans/mc_trans;
        break;
      }
      if(tot_xlen <= 0)
        ABSORB;
      mc_interact = 1 - mc_trans;
      if(mc_interact <= 0)        /* Protect against rounding errors */
        break;
      if (!event_counter) p *= fabs(1 - p_trans)/mc_interact;
      /* Select a point at which to scatter the neutron, taking
         secondary extinction into account. */
      /* dP(l) = exp(-tot_xlen*l)dl
         P(l<l_0) = [-1/tot_xlen*exp(-tot_xlen*l)]_0^l_0
                  = (1 - exp(-tot_xlen*l0))/tot_xlen
         l = -log(1 - tot_xlen*rand0max(P(l<l_full)))/tot_xlen
       */
      if(tot_xlen*l_full < 1e-6)
        /* For very weak scattering, use simple uniform sampling of scattering
           point to avoid rounding errors. */
        l = rand0max(l_full);
      else
        l = -log(1 - rand0max((1 - exp(-tot_xlen*l_full))))/tot_xlen;
      PROP_DT(l/v);
      event_counter++;

      /* (4). Account for the probability of absorbtion */
      p *= (coh_xlen + inc_xlen)/tot_xlen;
      /* Choose between coherent and incoherent scattering */
      if(coh_xlen == 0 || rand0max(coh_xlen + inc_xlen) <= inc_xlen)
      {
        /* (6). Incoherent scattering */
        randvec_target_circle(&kix, &kiy, &kiz, NULL, vx, vy, vz, 0);
        vx = kix; /* ki vector is used as tmp var with norm v */
        vy = kiy;
        vz = kiz; /* Go for next scattering event */
      } else {
        /* 7. Coherent scattering. Select reciprocal lattice point. */
        if(coh_refl <= 0)
          ABSORB;
        r = rand0max(coh_refl);
        sum = 0;
        for(j = 0; j < tau_count; j++)
        {
          sum += T[j].refl;
          if(sum > r) break;
        }
        if(j >= tau_count)
        {
          fprintf(stderr, "Single_crystal: Error: Illegal tau search "
              "(r = %g, sum = %g).\n", r, sum);
          j = tau_count - 1;
        }
        i = T[j].index;
        /* (8). Pick scattered wavevector kf from 2D Gauss distribution. */
        z1 = randnorm();
        z2 = randnorm();
        y1 = T[j].l11*z1 + T[j].y0x;
        y2 = T[j].l12*z1 + T[j].l22*z2 + T[j].y0y;
        kfx = T[j].rho_x + T[j].ox + T[j].b1x*y1 + T[j].b2x*y2;
        kfy = T[j].rho_y + T[j].oy + T[j].b1y*y1 + T[j].b2y*y2;
        kfz = T[j].rho_z + T[j].oz + T[j].b1z*y1 + T[j].b2z*y2;
        /* Normalize kf to length of ki, to account for planer
          approximation of the Ewald sphere. */
        adjust = ki/sqrt(kfx*kfx + kfy*kfy + kfz*kfz);
        kfx *= adjust;
        kfy *= adjust;
        kfz *= adjust;
        /* Adjust neutron weight (see manual for explanation). */
        p *= T[j].xsect*coh_refl/(coh_xsect*T[j].refl);
        vx = K2V*(L[i].u1x*kfx + L[i].u2x*kfy + L[i].u3x*kfz);
        vy = K2V*(L[i].u1y*kfx + L[i].u2y*kfy + L[i].u3y*kfz);
        vz = K2V*(L[i].u1z*kfx + L[i].u2z*kfy + L[i].u3z*kfz);
      }
      SCATTER;
      /* exit if multiple scattering order has been reached */
      if (order && event_counter >= order) break;
      /* Repeat loop for next scattering event. */
    } /* end for (;;) */
  } /* box intsresect */
%}

MCDISPLAY
%{
  double xmin = -0.5*xwidth;
  double xmax =  0.5*xwidth;
  double ymin = -0.5*yheight;
  double ymax =  0.5*yheight;
  double zmin = -0.5*zthick;
  double zmax =  0.5*zthick;
  magnify("xyz");
  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