/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Single_crystal
*
* %Identification
* Written by: Kristian Nielsen
* Date: December 1999
* Origin: Risoe
* Modified by: EF, 22nd Apr 2003 : now uses Read_Table library
* Modified by: EBK 2010: adapted for x-rays:
* Modified by: M. Schulz, March 2012 : allow to curve the crystal planes
* Modified by: EF, PW, May 2014: code efficiency improvement when SPLIT is used
* Modified by: EF, PW, 2015: powder and texture mode
* Modified by: PW, May 2017: Remove statement about being under validation
* Modified by: PW, June 2017: Doc updates
* Modified by: PW, Feb 2018: GPU edits
* Modified by: EF, May 2023: (re)port from n to X
* Modified by: EF, June 2023: can now read CIF files, via cif2hkl
*
* Mosaic single crystal with multiple scattering vectors, optimised for speed
* with large crystals and many reflections.
*
* %Description
* 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_a, mosaic_b, and mosaic_c
* parameters.
* The crystal lattice can be bent locally, keeping the external geometry unchanged.
* Curvature is spherical along vertical and horizontal axes.
*
* <b>Speed/stat optimisation using SPLIT</b>
* In order to dramatically improve the simulation efficiency, we recommend to
* use a SPLIT keyword on this component (or prior to it), as well as to disable
* the multiple scattering handling by setting order=1. This is especially powerful
* for large reflection lists such as with macromolecular proteins. When an incoming
* particle is identical to the preceeding, reciprocal space initialisation is 
* skipped, and a Monte Carlo choice is done on available reflections from the last
* repciprocal space calculation! To assist the user in choosing a "relevant" value
* of the SPLIT, a rolling average of the number of available reflections is
* calculated and presented in the component output.
*
* <b>Mosacitiy modes</b>
* The component features three independent ways of parametrising mosaicity:
*  a) The original algorithm where mosaicity is implemented by extending each 
*     reflection by a Gaussian "cigar" in reciprocal space, characterised by
*     the parameters mosaic and delta_d_d. 
*     (Also known as "isotropic mosaicity").
*  b) A similar mode where mosaicities can be non-isotropic and given as the
*     parameters mosaic_a, mosaic_b and mosaic_c, around the unit cell axes.
*     (Also known as "anisotropic mosaicity").
*  c) Given two "macroscopically"/experimentally measured width/mosaicities 
*     of two independent reflections, parametrised by the list 
*     mosaic_AB = {mos_a, mos_b, a_h, a_k, a_l, b_h, b_k, b_l}, a set of 
*     microscopic mosaicities as in b) are estimated (internally) and applied.
*     (Also known as "phenomenological mosaicity").
*     
* <b>Powder-mode</b>
* When these two modes are used (powder=1 ), a randomised transformation
* of the particle direction is made before and after scattering, thereby letting 
* the single crystal behave as a crystallite of either a powder (crystallite
* orientation fully randomised).
* 
* <b>Curved crystal mode</b>
* The component features a method to curve the lattice planes slightly with respect
* to the outer geometry of the crystal. The method is implemented as a transformation
* on the particle direction vector, and should be used only in cases where:
*  a) The reflection lattice vector is ~ orthogonal to the crystal surface.
*  b) The modelled curvarture is "small" with respect to the crystal surface.
* 
* <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
*   cylinder:        radius x yheight
*   sphere:          radius (yheight=0)
*   any shape:       geometry=OFF/PLY file
*
*   The complex geometry option handles any closed non-convex polyhedra.
*   It computes the intersection points of the photon 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 on the OFF/PLY file data, but its
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* <b>Crystal definition file format</b>
* 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 7-th 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.
* The reflection list should be ordered by decreasing d-spacing values.
* Lines begining by '#' are read as comments (ignored). Most sample parameters
* may be defined from the data file header, following the same mechanism as
* PowderN.
*
* Current data file header keywords include, for data format specification:
*    #column_h <index of the Bragg Qh column>
*    #column_k <index of the Bragg Qk column>
*    #column_l <index of the Bragg Ql column>
*    #column_F2 <index of the squared str. factor '|F|^2' column [b]>
*    #column_F  <index of the structure factor norm '|F|' column>
* and for material specification:
*    #sigma_inc <value of incoherent cross section [barns]>
*    #Delta_d/d <value of Delta_d/d width for all lines>
*    #lattice_a <value of the a lattice parameter [Angs]>
*    #lattice_b <value of the b lattice parameter [Angs]>
*    #lattice_c <value of the c lattice parameter [Angs]>
*    #lattice_aa <value of the alpha lattice angle [deg]>
*    #lattice_bb <value of the beta  lattice angle [deg]>
*    #lattice_cc <value of the gamma lattice angle [deg]>
*
* Last, CIF, FullProf and ShelX files can be read, and converted to F2(hkl) lists
* when 'cif2hkl' is installed. The CIF2HKL env variable can be used to point to a
* proper executable, else the McCode or the system installed versions are used.
*
* <b>Satellite Bragg peaks - surface crystal truncation rods (CTR)</b>
* It is known that scattering from a finite crystal introduces a broadening of
* Bragg peaks, seen in surface diffraction at grazing angle [Robinson and Tweet, 
* Rep. Prog. Phys. 55 (1992) 599)]. The CTR is specified as two vectors, which 
* hold the squared Fourier transform of the crystal geometry, for instance:
*   bulk:       Dirac peak      (FT of infinity)
*   half-bulk:  Dirac+1/k^2     (FT of half plane, k in rlu)
*   layer:      sinc(PI*k*d)^2  (FT of a top-hat, thickness 'd')
* These vectors are given as 'surf_k' [in 1/Angs] and 'surf_FT2', both of length
* 'surf_size'. The CTR is to be applied along vector 'surf_dir' which indicates 
* the surface normal {nx,ny,nz} in real space. When surf_size=-1, the component 
* sets the proper truncation function (only for thin box and disk shapes).
* This feature is an approximation of the real surface scattering, and does not
* handle complex geometries (clusters, wetting, multi-layers, ...). It may be 
* used as well to describe over-structure satellite peaks.
*
* Example: Single_crystal(xwidth=0.01, yheight=0.01, zdepth=0.01, mosaic = 5, 
*   reflections="Si.lau")
*
* A diamond crystal plate, cut for (002) reflections
*   Single_crystal(xwidth = 0.002, yheight = 0.1, zdepth = 0.1,
*     mosaic = 30, reflections = "C-diamond.lau",
*     ax=0,      ay=2.14,   az=-1.24,
*     bx = 0,    by = 0,    bz =  2.47,
*     cx = 6.71, cy = 0,    cz =  0)
*
* A adrenaline protein
*   Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005,
*     mosaic = 5, reflections="adrenaline.lau")
*
* Also, always use a non-zero value of delta_d_d.
*
* This sample component can advantageously benefit from the SPLIT feature, e.g.
* SPLIT COMPONENT sx = Single_crystal(...)
*
* %Parameters
* INPUT PARAMETERS
* radius:             [m] Outer radius of sample in (x,z) plane
* xwidth:             [m] Width of crystal
* yheight:            [m] Height of crystal
* zdepth:             [m] Depth of crystal (no extinction simulated)
* 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
* delta_d_d:          [1] Lattice spacing variance, gaussian RMS
* mosaic:   [arc minutes] Crystal mosaic (isotropic), gaussian RMS. Puts the crystal in the isotropic mosaic model state, thus disregarding other mosaicity parameters.
* mosaic_a: [arc minutes] Horizontal (rotation around lattice vector a) mosaic (anisotropic), gaussian RMS. Put the crystal in the anisotropic crystal vector state. I.e. model mosaicity through rotation around the crystal lattice vectors. Has precedence over in-plane mosaic model.
* mosaic_b: [arc minutes] Vertical (rotation around lattice vector b) mosaic (anisotropic), gaussian RMS.
* mosaic_c: [arc minutes] Out-of-plane (Rotation around lattice vector c) mosaic (anisotropic), gaussian RMS
* mosaic_AB: [arc_minutes, arc_minutes,1, 1, 1, 1, 1, 1]  In Plane mosaic rotation and plane vectors (anisotropic), mosaic_A, mosaic_B, A_h,A_k,A_l, B_h,B_k,B_l. Puts the crystal in the in-plane mosaic state. Vectors A and B define plane in which  the crystal roation is defined, and mosaic_A, mosaic_B, denotes the resp. mosaicities (gaussian RMS) with respect to the two reflections chosen by A and B (Miller indices).
*
* recip_cell:        [1] Choice of direct/reciprocal (0/1) unit cell definition
* ax:      [AA or AA^-1] Coordinates of first (direct/recip) unit cell vector 
* ay:      [AA or AA^-1]   a on y axis
* az:      [AA or AA^-1]   a on z axis
* bx:      [AA or AA^-1] Coordinates of second (direct/recip) unit cell vector
* by:      [AA or AA^-1]   b on y axis
* bz:      [AA or AA^-1]   b on z axis
* cx:      [AA or AA^-1] Coordinates of third (direct/recip) unit cell vector
* cy:      [AA or AA^-1]   c on y axis
* cz:      [AA or AA^-1]   c on z axis
* reflections:  [string] File name containing structure factors of reflections (LAZ LAU CIF, FullProf, ShelX). Use empty ("") or NULL for incoherent scattering only
* order:             [1] Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...)
* extra_order:       [1] When using order, allow additional multiple scattering without coherent scattering, sensible with very large unit cells (0: disable, 1: one extra, 2: two extra, ...)
*
* Optional input parameters
*
* p_transmit:        [1] Monte Carlo probability for photons to be transmitted without any scattering. Used to improve statistics from weak reflections
* sigma_inc:     [barns] Incoherent scattering cross-section per unit cell (uniform). Fully isotropic and constant. Use -1 to unactivate
* aa:              [deg] Unit cell angles alpha, beta and gamma. Then uses norms of vectors a,b and c as lattice parameters
* bb:              [deg] Beta angle
* cc:              [deg] Gamma angle
* barns:             [1] Flag to indicate if |F|^2 from 'reflections' is in barns or fm^2. barns=1 for laz and isotropic constant elastic scattering (reflections=NULL), barns=0 for lau type files
* RX:                [m] Radius of horizontal along X lattice curvature. flat for 0
* RY:                [m] Radius of vertical   along Y lattice curvature. flat for 0
* powder:            [1] Flag to indicate powder mode, for simulation of Debye-Scherrer cones via random crystallite orientation. A powder texture can be approximated with powder within 0-1
* deltak:         [AA-1] Equality-threshold for use in SPLIT settings. If difference between all ki_{x,y,z} are less than deltak from previous particle, the two are considered alike enough to jump directly to the MC choice between 'active' reflections 
* material_datafile: [Be.txt] File where the material parameters for the absorption may be found. Format is similar to what may be found off the NIST website. 
* surf_size:         [1] Length of the surf_k and surf_FT vectors. When set as -1, CTR is automatically set for the box/thin disk geometry.
* surf_k:       [1/Angs] Momentum 'k' distribution around 0 for the CTR, length 'surf_size'.
* surf_FT2:          [1] Intensity |FT(real space)|^2 distribution as CTR of a single Bragg peak, length 'surf_size'.
* surf_dir:   [nx,ny,nz] Normal vector in real space for the truncation rod spread direction. {1,0,0} is for an YZ surface, with normal along X.
*
* CALCULATED PARAMETERS:
*
* hkl_info: [structure]  Internal
* hkl_info.type: interaction type of event 't'=Transmit, 'i'=Incoherent, 'c'=Coherent [char]
* hkl_info.h:
* hkl_info.k:    wavevector indices of last coherent scattering event [Angs-1]
* hkl_info.l:
*
* %Link
* See <a href="http://icsd.ill.fr">ICSD</a> Inorganic Crystal Structure Database
* %Link
* <a href="http://www.webelements.com/">Web Elements</a>
* %Link
* <a href="http://www.ill.eu/sites/fullprof/index.html">Fullprof</a> powder refinement
* %Link
* <a href="http://www.crystallographica.com/">Crystallographica</a> software
* %Link
* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a>
* %Link
* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a>
* %Link
* <a href="http://qhull.org">qhull</a>
* %Link
* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a>
* %Link
* material datafile obtained from http://physics.nist.gov/cgi-bin/ffast/ffast.pl
* %Link
* cif2hkl https://gitlab.com/soleil-data-treatment/soleil-software-projects/cif2hkl
*
* %End
****************************************************************************/

/*
Overview of algorithm:

(1). The photon 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 sections 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 photon 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.

*/

DEFINE COMPONENT Single_crystal

SETTING PARAMETERS(string reflections=0, string geometry=0, vector mosaic_AB={0,0, 0,0,0, 0,0,0},
  xwidth=0, yheight=0, zdepth=0, radius=0, delta_d_d=1e-4,
  mosaic = -1, mosaic_a = -1, mosaic_b = -1, mosaic_c = -1,
  recip_cell=0, barns=0,
  ax = 0, ay = 0, az = 0,
  bx = 0, by = 0, bz = 0,
  cx = 0, cy = 0, cz = 0,
  p_transmit = 0.001, sigma_inc = 0,
  aa=0, bb=0, cc=0, order=1, extra_order=0, RX=0, RY=0, powder=0,
  deltak=1e-6, string material_datafile="Si.txt",
  int surf_size=0, vector surf_k=NULL, vector surf_FT2=NULL, vector surf_dir={0,0,0})

/* photon parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
SHARE
%{
/* used for reading data table from file */
%include "read_table-lib"
%include "interoff-lib"
/* Declare structures and functions only once in each instrument. */
#ifndef SINGLE_CRYSTAL_DECL
#define SINGLE_CRYSTAL_DECL

#ifndef Mosaic_AB_Undefined
#define Mosaic_AB_Undefined {0,0, 0,0,0, 0,0,0}
#endif

#ifndef MCSX_REFL_SLIST_SIZE
#define MCSX_REFL_SLIST_SIZE 128
#endif

  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 sig123;              /* The product sig1*sig2*sig3 = volume of spot */
      double m1, m2, m3;          /* Diagonal matrix representation of Gauss */
      double cutoff;              /* Cutoff value for Gaussian tails */
    };

  struct tau_data
    {
      int index;                  /* Index into reflection table */
      double refl;
      double xsect;
      /* The following vectors are in local koordinates. */
      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 b1x, b1y, b1z;       /* Spanning vectors of Ewald sphere tangent */
      double b2x, b2y, b2z;
      double l11, l12, l22;       /* Cholesky decomposition L of 2D Gauss */
      double y0x, y0y;            /* 2D Gauss center in tangent plane */
    };

  struct hkl_info_struct
    {
      int    count;               /* Number of reflections */
      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 m_a, m_b, m_c;       /* length of lattice parameter lengths */
      double m_aa, m_bb, m_cc;    /* lattice angles */
      double sigma_i;             /* inc X sect */
      double rho;                 /* density [at/Angs^3] */
      double mat_weight;          /* material atomic weight [g/mol] */
      double mat_density;         /* material atomic density [g/cm3] */
      double V0;                  /* Unit cell volume (AA**3) */
      int    column_order[5];     /* column signification [h,k,l,F,F2] */
      int    recip;               /* Flag to indicate if recip or direct cell axes given */
      int    shape;               /* 0:cylinder, 1:box, 2:sphere 3:any shape*/
      int    flag_warning;        /* number of warnings */
      char   type;                /* type of last event: t=transmit,c=coherent or i=incoherent */
      int    h,k,l;               /* last coherent scattering momentum transfer indices */
      int    tau_count;           /* Number of reflections within cutoff */
      double coh_refl, coh_xsect; /* cross section computed with last tau_list */
      double kix, kiy, kiz;       /* last incoming photon ki */
      int    nb_reuses, nb_refl, nb_refl_count;
      int    max_tau_count;
      int    ctr_size;
      double *ctr_k;
      double *ctr_FT2;
      double *ctr_dir;
    };
    
#pragma acc routine
  // sign = SX_list_compare(a,b) 
  //   used for qsort, to sort reflections
  // input:  a,b: two 'line_data' reflection pointers.
  // output: -1,0,1 for a<b, a=b, a>b
  int SX_list_compare (void const *a, void const *b)
  {
     struct hkl_data const *pa = a;
     struct hkl_data const *pb = b;
     double s = pa->tau - pb->tau;

     if (!s) return 0;
     else    return (s < 0 ? -1 : 1);
  } /* SX_list_compare */
  
#ifndef CIF2HKL
#define CIF2HKL
  // hkl_filename = cif2hkl(file, options)
  //   used to convert CIF/CFL/INS file into F2(hkl)
  //   the CIF2HKL env var can point to a cif2hkl executable
  //   else the McCode binary is attempted, then the system.
  char *cif2hkl(const char *infile, const char *options) {
    char cmd[1024];
    int  ret = 0;
    int  found = 0;
    char *OUTFILE;
    
    // get filename extension
    const char *ext = strrchr(infile, '.');
    if(!ext || ext == infile) return infile;
    else ext++;
    
    // return input when no extension or not a CIF/FullProf/ShelX file
    if ( strcasecmp(ext, "cif") 
      && strcasecmp(ext, "pcr")
      && strcasecmp(ext, "cfl")
      && strcasecmp(ext, "shx")
      && strcasecmp(ext, "ins")
      && strcasecmp(ext, "res")) return infile;
      
    OUTFILE = malloc(1024);
    if (!OUTFILE) return infile;
    
    strncpy(OUTFILE, tmpnam(NULL), 1024); // create an output temporary file name
    
    // try in order the CIF2HKL env var, then the system cif2hkl, then the McCode one
    if (!found && getenv("CIF2HKL")) {
      snprintf(cmd,  1024, "%s -o %s %s %s",
        getenv("CIF2HKL"),
        OUTFILE, options, infile);
      ret = system(cmd);
      if (ret != -1 && ret != 127) found = 1;
    }
    if (!found) {
      // try with cif2hkl command from the system PATH
      snprintf(cmd,  1024, "%s -o %s %s %s",
        "cif2hkl", OUTFILE, options, infile);
      ret = system(cmd);
      if (ret != -1 && ret != 127) found = 1;
    }
    if (!found) {
      // As a last resort, attempt with cif2hkl from $MCXTRACE/bin
      snprintf(cmd,  1024, "%s%c%s%c%s -o %s %s %s",
        getenv(FLAVOR_UPPER) ? getenv(FLAVOR_UPPER) : MCXTRACE,
        MC_PATHSEP_C, "bin", MC_PATHSEP_C, "cif2hkl",
        OUTFILE, options, infile);
      ret = system(cmd);
    }
    // ret = -1:  child process could not be created
    // ret = 127: shell could not be executed in the child process     
    if (ret == -1 || ret == 127) {
      free(OUTFILE);
      return(NULL);
    }
      
    // test if the result file has been created
    FILE *file = Open_File(OUTFILE,"r", NULL);
    if (!file) return(NULL);     
    MPI_MASTER(
    printf("%s: INFO: Converting %s into F2(HKL) list %s\n", 
      __FILE__, infile, OUTFILE);
    printf ("%s\n",cmd);
    );
    fflush(NULL);
    return(OUTFILE);
  } // cif2hkl
#endif

  /* ------------------------------------------------------------------------ */
  
  // count = read_hkl_data(filename, &info)
  //   read the given file to store F2(HKL) reflections and metadata.
  int
  read_hkl_data(char *SC_file, struct hkl_info_struct *info, struct hkl_data **hkl_list,
      double SC_mosaic, double SC_mosaic_a, double SC_mosaic_b, double SC_mosaic_c, double *SC_mosaic_AB,
      int surf_size, double *surf_k, double *surf_FT2, double *surf_dir)
  {
    struct hkl_data *list = NULL;
    int    size = 0;
    t_Table sTable; /* sample data table structure from SC_file */
    int i=0, j=0, J=0;
    double tmp_x, tmp_y, tmp_z;
    char **parsing;
    char   flag=0;
    double nb_atoms=1;
    double sum_F2=0;
    char  *filename=NULL;

    if (!SC_file || !strlen(SC_file) || !strcmp(SC_file,"NULL") || !strcmp(SC_file,"0"))
      return(0);
      
    filename = cif2hkl(SC_file, "--xtal --mode XRA");
    
    if (!filename || Table_Read(&sTable, filename, 1) < 0) return(0);
    
    if (sTable.columns < 4) {
      MPI_MASTER(
        fprintf(stderr, "Single_crystal: %s: ERROR: The number of columns in %s should be at least %d for [h,k,l,F2]\n", __FILE__, SC_file, 4);
      );
      return(0);
    }
    if (!sTable.rows) {
      MPI_MASTER(
        fprintf(stderr, "Single_crystal: %s: ERROR: The number of rows in %s should be at least %d\n", __FILE__, SC_file, 1);
      );
      return(0);
    } else size = sTable.rows;

    /* parsing of header */
    parsing = Table_ParseHeader(sTable.header,
      "sigma_abs","sigma_a ", // 0-1
      "sigma_inc","sigma_i ",
      "column_h",
      "column_k",
      "column_l",
      "column_F ",
      "column_F2",
      "Delta_d/d", // 9
      "lattice_a ",
      "lattice_b ",
      "lattice_c ",
      "lattice_aa",
      "lattice_bb",
      "lattice_cc",
      "nb_atoms","multiplicity", // 17
      "Vc","V_0","V_rho","density","weight",
      NULL);

    if (parsing) {
      //if (parsing[0] && !info->sigma_a) info->sigma_a=atof(parsing[0]);
      //if (parsing[1] && !info->sigma_a) info->sigma_a=atof(parsing[1]);
      if (parsing[2] && !info->sigma_i) info->sigma_i=atof(parsing[2]);
      if (parsing[3] && !info->sigma_i) info->sigma_i=atof(parsing[3]);
      if (parsing[4])                   info->column_order[0]=atoi(parsing[4]);
      if (parsing[5])                   info->column_order[1]=atoi(parsing[5]);
      if (parsing[6])                   info->column_order[2]=atoi(parsing[6]);
      if (parsing[7])                   info->column_order[3]=atoi(parsing[7]);
      if (parsing[8])                   info->column_order[4]=atoi(parsing[8]);
      if (parsing[9] && info->m_delta_d_d <0) info->m_delta_d_d=atof(parsing[9]);
      if (parsing[10] && !info->m_a)    info->m_a =atof(parsing[10]);
      if (parsing[11] && !info->m_b)    info->m_b =atof(parsing[11]);
      if (parsing[12] && !info->m_c)    info->m_c =atof(parsing[12]);
      if (parsing[13] && !info->m_aa)   info->m_aa=atof(parsing[13]);
      if (parsing[14] && !info->m_bb)   info->m_bb=atof(parsing[14]);
      if (parsing[15] && !info->m_cc)   info->m_cc=atof(parsing[15]);
      if (parsing[16])                  nb_atoms=atof(parsing[16]);
      if (parsing[17])                  nb_atoms=atof(parsing[17]);
      if (parsing[18] && !info->rho)    info->rho    =1/atof(parsing[18]);
      if (parsing[19] && !info->rho)    info->rho    =1/atof(parsing[19]);
      if (parsing[20] && !info->rho)    info->rho    =atof(parsing[20]);
      if (parsing[21] && !info->mat_density) info->mat_density=atof(parsing[21]);
      if (parsing[22] && !info->mat_weight)  info->mat_weight =atof(parsing[22]);
      
      for (i=0; i<=22; i++) if (parsing[i]) free(parsing[i]);
      free(parsing);
    }
    
    if (!info->rho && info->mat_density > 0 && info->mat_weight > 0 && nb_atoms > 0) {
      /* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */
      /* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */
      info->rho = info->mat_density/(info->mat_weight*nb_atoms)/1e24*NA;
    }

    if (nb_atoms > 1) { info->sigma_i *= nb_atoms; }

    /* special cases for the structure definition */
    if (info->m_ax || info->m_ay || info->m_az) {info->m_a=0; info->m_aa=0;} /* means we specify by hand the vectors */
    if (info->m_bx || info->m_by || info->m_bz) {info->m_b=0; info->m_bb=0;}
    if (info->m_cx || info->m_cy || info->m_cz) {info->m_c=0; info->m_cc=0;};

    /* compute the norm from vector a if missing */
    if (info->m_ax || info->m_ay || info->m_az) {
      double as=sqrt(info->m_ax*info->m_ax+info->m_ay*info->m_ay+info->m_az*info->m_az);
      if (!info->m_bx && !info->m_by && !info->m_bz) info->m_a=info->m_b=as;
      if (!info->m_cx && !info->m_cy && !info->m_cz) info->m_a=info->m_c=as;
    }
    if (info->m_a && !info->m_b) info->m_b=info->m_a;
    if (info->m_b && !info->m_c) info->m_c=info->m_b;

    /* compute the lattive angles if not set from data file. Not used when in vector mode. */
    if (info->m_a && !info->m_aa) info->m_aa=90;
    if (info->m_aa && !info->m_bb) info->m_bb=info->m_aa;
    if (info->m_bb && !info->m_cc) info->m_cc=info->m_bb;

    /* parameters consistency checks */
    if (!info->m_ax && !info->m_ay && !info->m_az && !info->m_a) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Wrong a lattice vector definition\n", __FILE__);
      );
      return(0);
    }
    if (!info->m_bx && !info->m_by && !info->m_bz && !info->m_b) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Wrong b lattice vector definition\n", __FILE__);
      );
      return(0);
    }
    if (!info->m_cx && !info->m_cy && !info->m_cz && !info->m_c) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Wrong c lattice vector definition\n", __FILE__);
      );
      return(0);
    }
    if (info->m_aa && info->m_bb && info->m_cc && info->recip) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Selecting reciprocal cell and angles is unmeaningful.\n", __FILE__);
      );
      return(0);
    }

    /* when lengths a,b,c + angles are given (instead of vectors a,b,c) */
    if (info->m_aa && info->m_bb && info->m_cc)
    {
      double as,bs,cs;
      if (info->m_a) as = info->m_a;
      else as = sqrt(info->m_ax*info->m_ax+info->m_ay*info->m_ay+info->m_az*info->m_az);
      if (info->m_b) bs = info->m_b;
      else bs = sqrt(info->m_bx*info->m_bx+info->m_by*info->m_by+info->m_bz*info->m_bz);
      if (info->m_c) cs = info->m_c;
      else cs =  sqrt(info->m_cx*info->m_cx+info->m_cy*info->m_cy+info->m_cz*info->m_cz);

      info->m_bz = as; info->m_by = 0; info->m_bx = 0;
      info->m_az = bs*cos(info->m_cc*DEG2RAD);
      info->m_ay = bs*sin(info->m_cc*DEG2RAD);
      info->m_ax = 0;
      info->m_cz = cs*cos(info->m_bb*DEG2RAD);
      info->m_cy = cs*(cos(info->m_aa*DEG2RAD)-cos(info->m_cc*DEG2RAD)*cos(info->m_bb*DEG2RAD))
                     /sin(info->m_cc*DEG2RAD);
      info->m_cx = sqrt(cs*cs - info->m_cz*info->m_cz - info->m_cy*info->m_cy);

      MPI_MASTER(
        printf("Single_crystal: %s structure a=%g b=%g c=%g aa=%g bb=%g cc=%g ",
        SC_file, as, bs, cs, info->m_aa, info->m_bb, info->m_cc);
      );
    } else {
      MPI_MASTER(
      if (!info->recip) {
        printf("Single_crystal: %s structure a=[%g,%g,%g] b=[%g,%g,%g] c=[%g,%g,%g] ",
	       SC_file, info->m_ax ,info->m_ay ,info->m_az,
	       info->m_bx ,info->m_by ,info->m_bz,
	       info->m_cx ,info->m_cy ,info->m_cz);
      } else {
        printf("Single_crystal: %s structure a*=[%g,%g,%g] b*=[%g,%g,%g] c*=[%g,%g,%g] ",
	       SC_file, info->m_ax ,info->m_ay ,info->m_az,
	       info->m_bx ,info->m_by ,info->m_bz,
	       info->m_cx ,info->m_cy ,info->m_cz);
      }
      );
    }
    /* Compute reciprocal or direct lattice vectors. */
    if (!info->recip) {
      vec_prod(tmp_x, tmp_y, tmp_z,
	       info->m_bx, info->m_by, info->m_bz,
	       info->m_cx, info->m_cy, info->m_cz);
      info->V0 = fabs(scalar_prod(info->m_ax, info->m_ay, info->m_az, tmp_x, tmp_y, tmp_z));
      if (!info->V0 || isnan(info->V0)) {
        MPI_MASTER(
          fprintf(stderr,
              "Single_crystal: %s: ERROR: Lattice cell volume is found invalid V0=%g. Check crystal cell parameters.\n", __FILE__, info->V0));
      }

      info->asx = 2*PI/info->V0*tmp_x;
      info->asy = 2*PI/info->V0*tmp_y;
      info->asz = 2*PI/info->V0*tmp_z;
      vec_prod(tmp_x, tmp_y, tmp_z, info->m_cx, info->m_cy, info->m_cz, info->m_ax, info->m_ay, info->m_az);
      info->bsx = 2*PI/info->V0*tmp_x;
      info->bsy = 2*PI/info->V0*tmp_y;
      info->bsz = 2*PI/info->V0*tmp_z;
      vec_prod(tmp_x, tmp_y, tmp_z, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz);
      info->csx = 2*PI/info->V0*tmp_x;
      info->csy = 2*PI/info->V0*tmp_y;
      info->csz = 2*PI/info->V0*tmp_z;
    } else {
      info->asx = info->m_ax;
      info->asy = info->m_ay;
      info->asz = info->m_az;
      info->bsx = info->m_bx;
      info->bsy = info->m_by;
      info->bsz = info->m_bz;
      info->csx = info->m_cx;
      info->csy = info->m_cy;
      info->csz = info->m_cz;

      vec_prod(tmp_x, tmp_y, tmp_z,
	       info->bsx/(2*PI), info->bsy/(2*PI), info->bsz/(2*PI),
	       info->csx/(2*PI), info->csy/(2*PI), info->csz/(2*PI));
	    info->V0 = scalar_prod(info->asx/(2*PI), info->asy/(2*PI), info->asz/(2*PI), tmp_x, tmp_y, tmp_z);
	    if (!info->V0 || isnan(info->V0)) {
        MPI_MASTER(
          fprintf(stderr,
              "Single_crystal: %s: ERROR: Lattice reciprocal cell volume is found invalid V0=%g. Check crystal cell parameters.\n", __FILE__, info->V0));
      }
      info->V0 = 1/fabs(info->V0);

      /*compute the direct cell parameters, for completeness*/
      info->m_ax = tmp_x*info->V0;
      info->m_ay = tmp_y*info->V0;
      info->m_az = tmp_z*info->V0;
      vec_prod(tmp_x, tmp_y, tmp_z,info->csx/(2*PI), info->csy/(2*PI), info->csz/(2*PI),info->asx/(2*PI), info->asy/(2*PI), info->asz/(2*PI));
      info->m_bx = tmp_x*info->V0;
      info->m_by = tmp_y*info->V0;
      info->m_bz = tmp_z*info->V0;
      vec_prod(tmp_x, tmp_y, tmp_z,info->asx/(2*PI), info->asy/(2*PI), info->asz/(2*PI),info->bsx/(2*PI), info->bsy/(2*PI), info->bsz/(2*PI));
      info->m_cx = tmp_x*info->V0;
      info->m_cy = tmp_y*info->V0;
      info->m_cz = tmp_z*info->V0;
    }
    MPI_MASTER(printf("V0=%g\n", info->V0););

    if (!info->column_order[0] || !info->column_order[1] || !info->column_order[2]) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Wrong h,k,l column definition\n", SC_file);
      );
      return(0);
    }
    if (!info->column_order[3] && !info->column_order[4]) {
      MPI_MASTER(
        fprintf(stderr,
              "Single_crystal: %s: ERROR: Wrong F,F2 column definition\n", SC_file);
      );
      return(0);
    }

    /* allocate hkl_data array */
    
    if (surf_size <= 0 || !surf_k || !surf_FT2 || !surf_dir) surf_size = 1;
    list = (struct hkl_data*) malloc(size*surf_size*sizeof(struct hkl_data));

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

      /* get data from table */
      h = Table_Index(sTable, i, info->column_order[0]-1);
      k = Table_Index(sTable, i, info->column_order[1]-1);
      l = Table_Index(sTable, i, info->column_order[2]-1);
      if (info->column_order[3])
      { F2= Table_Index(sTable, i, info->column_order[3]-1); F2 *= F2; }
      else if (info->column_order[4])
        F2= Table_Index(sTable, i, info->column_order[4]-1);

      list[J].h = h;
      list[J].k = k;
      list[J].l = l;
      list[J].F2 = F2;
      sum_F2    += F2;

      /* Precompute some values */
      list[J].tau_x = h*info->asx + k*info->bsx + l*info->csx;
      list[J].tau_y = h*info->asy + k*info->bsy + l*info->csy;
      list[J].tau_z = h*info->asz + k*info->bsz + l*info->csz;
      
      // create satellite peaks from e.g. surface truncation (CTR)
      if (surf_size > 1 && surf_FT2 && surf_k && surf_FT2[j] > 0 && surf_k[j]) {
        list[J].tau_x += surf_dir[0]*surf_k[j];
        list[J].tau_y += surf_dir[1]*surf_k[j];
        list[J].tau_z += surf_dir[2]*surf_k[j];
        list[J].F2    *= surf_FT2[j];
      }
      
      list[J].tau = sqrt(list[J].tau_x*list[J].tau_x +
                         list[J].tau_y*list[J].tau_y +
                         list[J].tau_z*list[J].tau_z);
      list[J].u1x = list[J].tau_x/list[J].tau;
      list[J].u1y = list[J].tau_y/list[J].tau;
      list[J].u1z = list[J].tau_z/list[J].tau;
      sig1 = FWHM2RMS*info->m_delta_d_d*list[J].tau;

      /* Find two arbitrary axes perpendicular to tau and each other. */
      normal_vec(&b1[0], &b1[1], &b1[2],
                 list[J].u1x, list[J].u1y, list[J].u1z);
      vec_prod(b2[0], b2[1], b2[2],
               list[J].u1x, list[J].u1y, list[J].u1z,
               b1[0], b1[1], b1[2]);

      /* Find the two mosaic axes perpendicular to tau. */
      if(SC_mosaic > 0) {
        /* Use isotropic mosaic. */
        list[J].u2x = b1[0];
        list[J].u2y = b1[1];
        list[J].u2z = b1[2];
        sig2 = FWHM2RMS*list[J].tau*MIN2RAD*SC_mosaic;
        list[J].u3x = b2[0];
        list[J].u3y = b2[1];
        list[J].u3z = b2[2];
        sig3 = FWHM2RMS*list[J].tau*MIN2RAD*SC_mosaic;
      } else if(SC_mosaic_a > 0 && SC_mosaic_b > 0 && SC_mosaic_c > 0) {
        /* Use anisotropic mosaic. */
        MPI_MASTER(
        fprintf(stderr,"Single_crystal: Warning: you are using an experimental feature:\n"
          "  anistropic mosaicity. Please examine your data carefully.\n");
        );
        /* compute the jacobian of (tau_v,tau_n) from rotations around the unit cell vectors. */
        struct hkl_data *l =&(list[J]);
        double xia_x,xia_y,xia_z,xib_x,xib_y,xib_z,xic_x,xic_y,xic_z;
        /*input parameters are in arc minutes*/
        double sig_fi_a=SC_mosaic_a*MIN2RAD;
        double sig_fi_b=SC_mosaic_b*MIN2RAD;
        double sig_fi_c=SC_mosaic_c*MIN2RAD;
        if(info->m_a==0) info->m_a=sqrt(scalar_prod( info->m_ax,info->m_ay,info->m_az,info->m_ax,info->m_ay,info->m_az));
        if(info->m_b==0) info->m_b=sqrt(scalar_prod( info->m_bx,info->m_by,info->m_bz,info->m_bx,info->m_by,info->m_bz));
        if(info->m_c==0) info->m_c=sqrt(scalar_prod( info->m_cx,info->m_cy,info->m_cz,info->m_cx,info->m_cy,info->m_cz));

        l->u2x = b1[0];
        l->u2y = b1[1];
        l->u2z = b1[2];
        l->u3x = b2[0];
        l->u3y = b2[1];
        l->u3z = b2[2];

        xia_x=l->tau_x-(M_2_PI*h/info->m_a)*info->asx;
        xia_y=l->tau_y-(M_2_PI*h/info->m_a)*info->asy;
        xia_z=l->tau_z-(M_2_PI*h/info->m_a)*info->asz;
        xib_x=l->tau_x-(M_2_PI*h/info->m_b)*info->bsx;
        xib_y=l->tau_y-(M_2_PI*h/info->m_b)*info->bsy;
        xib_z=l->tau_z-(M_2_PI*h/info->m_b)*info->bsz;
        xic_x=l->tau_x-(M_2_PI*h/info->m_c)*info->csx;
        xic_y=l->tau_y-(M_2_PI*h/info->m_c)*info->csy;
        xic_z=l->tau_z-(M_2_PI*h/info->m_c)*info->csz;

        double xia=sqrt(xia_x*xia_x + xia_y*xia_y + xia_z*xia_z);
        double xib=sqrt(xib_x*xib_x + xib_y*xib_y + xib_z*xib_z);
        double xic=sqrt(xic_x*xic_x + xic_y*xic_y + xic_z*xic_z);

        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u2x,l->u2y,l->u2z);
        double J_n_fia= xia/info->m_a/l->tau*scalar_prod(info->asx,info->asy,info->asz,tmp_x,tmp_y,tmp_z);
        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u2x,l->u2y,l->u2z);
        double J_n_fib= xib/info->m_b/l->tau*scalar_prod(info->bsx,info->bsy,info->bsz,tmp_x,tmp_y,tmp_z);
        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u2x,l->u2y,l->u2z);
        double J_n_fic= xic/info->m_c/l->tau*scalar_prod(info->csx,info->csy,info->csz,tmp_x,tmp_y,tmp_z);

        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u3x,l->u3y,l->u3z);
        double J_v_fia= xia/info->m_a/l->tau*scalar_prod(info->asx,info->asy,info->asz,tmp_x,tmp_y,tmp_z);
        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u3x,l->u3y,l->u3z);
        double J_v_fib= xib/info->m_b/l->tau*scalar_prod(info->bsx,info->bsy,info->bsz,tmp_x,tmp_y,tmp_z);
        vec_prod(tmp_x,tmp_y,tmp_z,l->tau_x,l->tau_y,l->tau_z, l->u3x,l->u3y,l->u3z);
        double J_v_fic= xic/info->m_c/l->tau*scalar_prod(info->csx,info->csy,info->csz,tmp_x,tmp_y,tmp_z);

        /*with the jacobian we can compute the sigmas in terms of the orthogonal vectors u2 and u3*/
        sig2=sig_fi_a*fabs(J_v_fia) + sig_fi_b*fabs(J_v_fib) + sig_fi_c*fabs(J_v_fic);
        sig3=sig_fi_a*fabs(J_n_fia) + sig_fi_b*fabs(J_n_fib) + sig_fi_c*fabs(J_n_fic);
      } else if (SC_mosaic_AB[0]!=0 && SC_mosaic_AB[1]!=0){
        if ( (SC_mosaic_AB[2]==0 && SC_mosaic_AB[3]==0 && SC_mosaic_AB[4]==0) || (SC_mosaic_AB[5]==0 && SC_mosaic_AB[6]==0 && SC_mosaic_AB[7]==0) ){
          MPI_MASTER(
          fprintf(stderr,"Single_crystal: %s: ERROR: in-plane mosaics are specified but one (or both)\n"
              "  in-plane reciprocal vector is the zero vector.\n", __FILE__);
          );
          return(0);
        }
        MPI_MASTER(
        fprintf(stderr,"Single_crystal: %s: Warning: you are using an experimental feature: \n"
              "  \"in-plane\" anistropic mosaicity. Please examine your data carefully.\n", __FILE__);
        );

        /*for given reflection in list - compute linear comb of tau_a and tau_b*/
        /*check for not in plane - f.i. check if (tau_a X tau_b).tau_i)==0*/
        struct hkl_data *l =&(list[J]);
        double det,c1,c2,sig_tau_c;
        double em_x,em_y,em_z, tmp_x,tmp_y,tmp_z;
        double tau_a[3],tau_b[3];
        /*convert Miller indices to taus*/
        if(info->m_a==0) info->m_a=sqrt(scalar_prod( info->m_ax,info->m_ay,info->m_az,info->m_ax,info->m_ay,info->m_az));
        if(info->m_b==0) info->m_b=sqrt(scalar_prod( info->m_bx,info->m_by,info->m_bz,info->m_bx,info->m_by,info->m_bz));
        if(info->m_c==0) info->m_c=sqrt(scalar_prod( info->m_cx,info->m_cy,info->m_cz,info->m_cx,info->m_cy,info->m_cz));
        tau_a[0]=M_2_PI*( (SC_mosaic_AB[2]/info->m_a)*info->asx + (SC_mosaic_AB[3]/info->m_b)*info->bsx + (SC_mosaic_AB[4]/info->m_c)*info->csx );
        tau_a[1]=M_2_PI*( (SC_mosaic_AB[2]/info->m_a)*info->asy + (SC_mosaic_AB[3]/info->m_b)*info->bsy + (SC_mosaic_AB[4]/info->m_c)*info->csy );
        tau_a[2]=M_2_PI*( (SC_mosaic_AB[2]/info->m_a)*info->asz + (SC_mosaic_AB[3]/info->m_b)*info->bsz + (SC_mosaic_AB[4]/info->m_c)*info->csz );
        tau_b[0]=M_2_PI*( (SC_mosaic_AB[5]/info->m_a)*info->asx + (SC_mosaic_AB[6]/info->m_b)*info->bsx + (SC_mosaic_AB[7]/info->m_c)*info->csx );
        tau_b[1]=M_2_PI*( (SC_mosaic_AB[5]/info->m_a)*info->asy + (SC_mosaic_AB[6]/info->m_b)*info->bsy + (SC_mosaic_AB[7]/info->m_c)*info->csy );
        tau_b[2]=M_2_PI*( (SC_mosaic_AB[5]/info->m_a)*info->asz + (SC_mosaic_AB[6]/info->m_b)*info->bsz + (SC_mosaic_AB[7]/info->m_c)*info->csz );

        /*check determinants to see how we should compute the linear combination of a and b (to match c)*/
        if ((det=tau_a[0]*tau_b[1]-tau_a[1]*tau_b[0])!=0){
          c1= (l->tau_x*tau_b[1] - l->tau_y*tau_b[0])/det;
          c2= (tau_a[0]*l->tau_y - tau_a[1]*l->tau_x)/det;
        }else if ((det=tau_a[1]*tau_b[2]-tau_a[2]*tau_b[1])!=0){
          c1= (l->tau_y*tau_b[2] - l->tau_z*tau_b[1])/det;
          c2= (tau_a[1]*l->tau_z - tau_a[2]*l->tau_y)/det;
        }else if ((det=tau_a[0]*tau_b[2]-tau_a[2]*tau_b[0])!=0){
          c1= (l->tau_x*tau_b[2] - l->tau_z*tau_b[0])/det;
          c2= (tau_a[0]*l->tau_z - tau_a[2]*l->tau_x)/det;
        }
        if ((c1==0) && (c2==0)){
          MPI_MASTER(
          fprintf(stderr,"Single_crystal: %s: Warning: reflection tau[%i]=(%g %g %g) "
          "has no component in defined mosaic plane\n", __FILE__,
          i, l->tau_x,l->tau_y,l->tau_z);
          );
        }
        /*compute linear combination => sig_tau_i = | c1*sig_tau_a + c2*sig_tau_b |  - also add in the minute to radian scaling factor*/;
        sig_tau_c = MIN2RAD*sqrt(c1*SC_mosaic_AB[0]*c1*SC_mosaic_AB[0] + c2*SC_mosaic_AB[1]*c2*SC_mosaic_AB[1]);
        l->u2x = b1[0]; l->u2y = b1[1]; l->u2z = b1[2];
        l->u3x = b2[0]; l->u3y = b2[1]; l->u3z = b2[2];

        /*so now let's compute the rotation around planenormal tau_a X tau_b*/
        /*g_bar (unit normal of rotation plane) = tau_a X tau_b / norm(tau_a X tau_b)*/
        vec_prod(tmp_x,tmp_y,tmp_z, tau_a[0],tau_a[1],tau_a[2],tau_b[0],tau_b[1],tau_b[2]);
        vec_prod(em_x,em_y,em_z, l->tau_x, l->tau_y, l->tau_z, tmp_x,tmp_y,tmp_z);
        NORM(em_x,em_y,em_z);
        sig2 = l->tau*sig_tau_c*fabs(scalar_prod(em_x,em_y,em_z, l->u2x,l->u2y,l->u2z));
        sig3 = l->tau*sig_tau_c*fabs(scalar_prod(em_x,em_y,em_z, l->u3x,l->u3y,l->u3z));
        /*protect against collapsing gaussians. These seem to be sensible values.*/
        if (sig2<1e-5) sig2=1e-5;
        if (sig3<1e-5) sig3=1e-5;
      }
      else {
        MPI_MASTER(
        fprintf(stderr,
                "Single_crystal: %s: ERROR: EITHER mosaic OR (mosaic_a, mosaic_b, mosaic_c)\n"
                "  must be given and be >0.\n", __FILE__);
        );
        return(0);
      }
      list[J].sig123 = sig1*sig2*sig3;
      list[J].m1 = 1/(2*sig1*sig1);
      list[J].m2 = 1/(2*sig2*sig2);
      list[J].m3 = 1/(2*sig3*sig3);
      /* Set Gauss cutoff to 5 times the maximal sigma. */
      if(sig1 > sig2)
        if(sig1 > sig3)
          list[J].cutoff = 5*sig1;
        else
          list[J].cutoff = 5*sig3;
      else
        if(sig2 > sig3)
          list[J].cutoff = 5*sig2;
        else
          list[J].cutoff = 5*sig3;
      if (!isnan(list[J].F2) && !isnan(list[J].cutoff) && !isnan(list[J].tau))
        J++; // switch to next HKL Bragg element in hkl_data
    } // for i, j
    Table_Free(&sTable);
    
    if (!sum_F2) {
      MPI_MASTER(
      exit(fprintf(stderr, "Single_crystal: %s: ERROR: all %i structure factors in file '%s' are null. Check the reflection list.\n",
        __FILE__, i, SC_file));
      );
    }
    
    // remove temporary F2(hkl) file when giving CFL/CIF/ShelX file
    if (filename != SC_file)
      unlink(filename);

    /* sort the list with increasing tau */
    qsort(list, J, sizeof(struct hkl_data),  SX_list_compare);

    *hkl_list = list;
    info->count = J;
    
    return(info->count);
  } /* read_hkl_data */

  /* ------------------------------------------------------------------------ */
  /* hkl_search
    search the HKL reflections which are on the Ewald sphere
    input:
      L,T,count,V0: constants for all calls
      kix,kiy,kiz: may be different for each call
    this function returns:
      tau_count (return), coh_refl, coh_xsect, T (updated elements in the array up to [j])
   */
#pragma acc routine
int hkl_search(struct hkl_data *L, void *TT, int count, double V0,
    double kix, double kiy, double kiz, double tau_max,
    double *coh_refl, double *coh_xsect)
  {
    double rho, rho_x, rho_y, rho_z;
    double diff;
    int    i,j;
    double ox,oy,oz;
    double b1x,b1y,b1z, b2x,b2y,b2z, kx, ky, kz, nx, ny, nz;
    double n11, n22, n12, det_N, inv_n11, inv_n22, inv_n12, l11, l22, l12,  det_L;
    double Bt_D_O_x, Bt_D_O_y, y0x, y0y, alpha;

    double ki = sqrt(kix*kix+kiy*kiy+kiz*kiz);
    int jglobal=-1;
    double coherent_refl,coherent_xsect;

    struct tau_data *T=(struct tau_data *)TT;

    //coherent_refl = *coh_refl;
    //coherent_xsect = *coh_xsect;
    coherent_refl = 0;
    coherent_xsect = 0;

    /* Common factor in coherent cross-section */
    double xsect_factor = pow(2*PI, 5.0/2.0)/(V0*ki*ki);
    j=0;
    for(i = 0; i < 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 = |ki - tau| = |kf| = |ki| within cutoff on the Ewald sphere */
        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. */
          kx = kix*L[i].u1x + kiy*L[i].u1y + kiz*L[i].u1z;
          ky = kix*L[i].u2x + kiy*L[i].u2y + kiz*L[i].u2z;
          kz = kix*L[i].u3x + kiy*L[i].u3y + kiz*L[i].u3z;
          T[j].rho_x = kx - L[i].tau;
          T[j].rho_y = ky;
          T[j].rho_z = kz;
          T[j].rho = rho;
          /* Compute the tangent plane of the Ewald sphere. */
          nx = T[j].rho_x/T[j].rho;
          ny = T[j].rho_y/T[j].rho;
          nz = T[j].rho_z/T[j].rho;
          ox = (ki - T[j].rho)*nx;
          oy = (ki - T[j].rho)*ny;
          oz = (ki - T[j].rho)*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, nx, ny, nz);
          vec_prod(b2x, b2y, b2z, nx, ny, 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;
          /* 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; /* intensity of that Bragg */
          if (isnan(T[j].refl)) continue;
          *coh_refl += T[j].refl;                                 /* total scatterable intensity */
          T[j].xsect = T[j].refl*L[i].F2;
          *coh_xsect += T[j].xsect;
          j++;
        }
        /*protect against tau shortlist buffer overrrun*/
        if (j==MCSX_REFL_SLIST_SIZE){
          break;
        }
      } /* end for */
        return (j); // this is 'tau_count', i.e. number of reachable reflections
    } /* end hkl_search */

#pragma acc routine
  int hkl_select(struct tau_data *T, int tau_count, double coh_refl, double *sum,_class_particle *_particle) {
      int j;
      double r = rand0max(coh_refl);
      *sum = 0;
      for(j = 0; j < tau_count; j++)
      {
        *sum += T[j].refl;
        if(*sum > r) break;
      }
      return j;
    }

    /* Functions for "reorientation", powder mode */
    /* Powder, forward */
#pragma acc routine
       void randrotate(double *nx, double *ny, double *nz, double a, double b, double c) {
      double x1, y1, z1, x2, y2, z2;
      rotate(x1, y1, z1, *nx,*ny,*nz, a, 1, 0, 0); /* <1> = rot(<n>,a) */
      rotate(x2, y2, z2,  x1, y1, z1, b, 0, 1, 0); /* <2> = rot(<1>,b) */
      rotate(*nx,*ny,*nz, x2, y2, z2, c, 0, 0, 1); /* <n> = rot(<2>,c) */
    }
    /* Powder, back */
#pragma acc routine
void randderotate(double *nx, double *ny, double *nz, double a, double b, double c) {
      double x1, y1, z1, x2, y2, z2;
      rotate(x1, y1, z1, *nx,*ny,*nz, -c, 0,0,1);
      rotate(x2, y2, z2,  x1, y1, z1, -b, 0,1,0);
      rotate(*nx,*ny,*nz, x2, y2, z2, -a, 1,0,0);
    }

#pragma acc routine
    /* rotate vector counterclockwise */
    void vec_rotate_2d(double* x, double* y, double angle) {
        double c, s;
        double newx, newy;

        c = cos(angle);
        s = sin(angle);

        newx = *x*c - *y*s;
        newy = *x*s + *y*c;

        *x = newx;
        *y = newy;
    }
  
  struct sx_abs_data
  {
    int     muc; // column where mu is to be found
    double  atomic_weight;
    double  delta_prefactor;
    t_Table table;
  };
  

  int sx_read_abs_data(char *ABS_file, struct sx_abs_data *abs, struct hkl_info_struct info ){
    int    status;
    char **parsing;
    double Z;
    
    if (ABS_file && strlen(ABS_file) && strcmp(ABS_file, "NULL") && strcmp(ABS_file, "0")) {
      if ( (status=Table_Read(&(abs->table),ABS_file,0))==-1){
        MPI_MASTER(
        fprintf(stderr,"Single_crystal: %s: ERROR: Could not parse file \"%s\"\n", __FILE__, ABS_file);
        );
        exit(-1);
      }
      abs->muc = (abs->table.columns==3 ? 1 : 5);
      
      parsing=Table_ParseHeader(abs->table.header,"Z","A[r]","rho",NULL);
      if (parsing) {
        if (parsing[0])                     Z               =strtol(parsing[0],NULL,10);
        if (parsing[1] && !info.mat_weight) info.mat_weight =strtod(parsing[1],NULL); // aka A[r] in g/mol
        if (parsing[2] && !info.rho)        info.rho        =strtod(parsing[2],NULL);
      }
  
      if (info.rho && info.mat_weight)
        abs->delta_prefactor= NA*(info.rho*1e-24)/info.mat_weight * 2.0*M_PI*RE;
      
      abs->atomic_weight=info.mat_weight;     
      
      return 1;
    } else {
      // init an empty table with absorption 0
      Table_Init(&(abs->table),2,2);
      abs->table.data[0]=0;abs->table.data[1]=0;
      abs->table.data[2]=FLT_MAX;abs->table.data[3]=0;
      fprintf(stderr,"Single_crystal: %s: Warning: material file %s not found. Absorption set to 0\n",
        __FILE__, ABS_file);
        return 0;
    }
  } // sx_read_abs_data
  
  
#endif /* !SINGLE_CRYSTAL_DECL */
%}


DECLARE
%{
  struct hkl_info_struct hkl_info;
  off_struct             offdata;
  struct hkl_data       *hkl_list;
  struct tau_data        tau_list[MCSX_REFL_SLIST_SIZE];
  struct sx_abs_data     abs_info;
%}

INITIALIZE
%{

  double as, bs, cs;
  int i=0;

  /* transfer input parameters */
  hkl_info.m_delta_d_d = delta_d_d;
  hkl_info.m_a  = 0;
  hkl_info.m_b  = 0;
  hkl_info.m_c  = 0;
  hkl_info.m_aa = aa;
  hkl_info.m_bb = bb;
  hkl_info.m_cc = cc;
  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;
 
  hkl_info.sigma_i = sigma_inc;
  hkl_info.recip   = recip_cell;

  /* default format h,k,l,F,F2  */
  hkl_info.column_order[0]=1;
  hkl_info.column_order[1]=2;
  hkl_info.column_order[2]=3;
  hkl_info.column_order[3]=0;
  hkl_info.column_order[4]=7;
  hkl_info.kix = hkl_info.kiy = hkl_info.kiz = 0;
  hkl_info.nb_reuses = hkl_info.nb_refl = hkl_info.nb_refl_count = 0;
  hkl_info.tau_count = 0;
  
  /*this should not be in hkl_info*/
  hkl_info.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
    MPI_MASTER(
    fprintf(stderr,"Single_crystal: %s: ERROR: You are attempting to use an OFF/PLY geometry without -DUSE_OFF. You will need to recompile with that define set!\n", NAME_CURRENT_COMP);
    );
    exit(-1);
    #else
    if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
      hkl_info.shape=3;
    }
    #endif
  }
  else if (xwidth && yheight && zdepth)  hkl_info.shape=1; /* box */
  else if (radius > 0 && yheight)        hkl_info.shape=0; /* cylinder */
  else if (radius > 0 && !yheight)       hkl_info.shape=2; /* sphere */

  if (hkl_info.shape < 0)
    exit(fprintf(stderr,"Single_crystal: %s: sample has invalid dimensions.\n"
                        "ERROR           Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
  
  /* determine if we have a surface broadening/truncation rods (CTR) */
  hkl_info.ctr_size=surf_size;
  hkl_info.ctr_k   =surf_k;
  hkl_info.ctr_FT2 =surf_FT2;
  hkl_info.ctr_dir =surf_dir;
  if (hkl_info.ctr_size < 0) { /* CTR automatic mode */
    double surf_dir_sum  =0;
    double surf_thickness=0;
    for (i=0; i<3; hkl_info.ctr_dir[i++]=0);
    if (hkl_info.shape == 1) {
      if      (xwidth <= yheight && xwidth <= zdepth)    { hkl_info.ctr_dir[0]=1; surf_thickness=xwidth;  }
      else if (yheight<= xwidth  && yheight<= zdepth)    { hkl_info.ctr_dir[1]=1; surf_thickness=yheight; }
      else                                               { hkl_info.ctr_dir[2]=1; surf_thickness=zdepth;  }
    } else if (hkl_info.shape == 0 && yheight <= radius) { hkl_info.ctr_dir[1]=1; surf_thickness=yheight; }
    // surf_size=-1: set automatically the k and FT2 vectors
    if (0 < surf_thickness) {
      hkl_info.ctr_size=21;                       // nb of satellite peaks
      double range_k = 1.0;               // k range [1/Angs], k=[-range/2 : +range/2]
      double dk = range_k/(hkl_info.ctr_size-1);  // delta-k for btw satellites
      DArray1d s_k = create_darr1d(hkl_info.ctr_size);
      DArray1d s_F = create_darr1d(hkl_info.ctr_size);
      hkl_info.ctr_k       = (double *)s_k;
      hkl_info.ctr_FT2     = (double *)s_F;
      for (i=0; i<hkl_info.ctr_size; i++) {
        hkl_info.ctr_k[i]   = -range_k/2+i*dk;    // 1/Angs, symmetric around 0 [1/Angs]
        if (surf_thickness < 1e-7) {      // thin layer thickness < 0.1 um -> sinc^2
          double x=PI*hkl_info.ctr_k[i]*surf_thickness*1e10;  // 1e10 [m -> Angs]
          hkl_info.ctr_FT2[i] = hkl_info.ctr_k[i] ? sin(x)/x : 1; 
        } else {                          // assume a half-bulk -> 1/k^2 [in rlu]
          hkl_info.ctr_FT2[i] = hkl_info.ctr_k[i] ? 1/(hkl_info.ctr_k[i]*2*PI) : 4*PI; // we assume a=2PI
        }
        hkl_info.ctr_FT2[i]*= hkl_info.ctr_FT2[i];        // squared
      }
      MPI_MASTER(
        printf("Single_crystal: %s: Info: Setting CTR surface normal XYZ=[%g,%g,%g], as a %s.\n", 
          NAME_CURRENT_COMP,
          hkl_info.ctr_dir[0],hkl_info.ctr_dir[1],hkl_info.ctr_dir[2], 
          surf_thickness < 1e-7 ? "thin layer (sinc^2)" : "half-bulk (1/k^2)");
      );
    }
  }
  if (hkl_info.ctr_size > 1 && (!hkl_info.ctr_dir || !hkl_info.ctr_dir[0] && !hkl_info.ctr_dir[1] && !hkl_info.ctr_dir[2])) { // no CTR normal ?
    MPI_MASTER(
    printf("Single_crystal: %s: Warning: No surf_dir surface normal set. Ignoring CTR.\n", i, NAME_CURRENT_COMP);
    );
    hkl_info.ctr_size = 1;
  }
  if (hkl_info.ctr_size > 1 && hkl_info.ctr_FT2 && hkl_info.ctr_k) { // do some CTR tests
    double sum_FT2=0, min_k=FLT_MAX, max_k=-FLT_MAX;
    for (i=0; i<hkl_info.ctr_size; i++) {
      if (hkl_info.ctr_FT2[i] < 0) {
        MPI_MASTER(
        printf("Single_crystal: %s: Warning: surf_FT2[%i] < 0. Must all be >= 0. Ignoring CTR.\n", i, NAME_CURRENT_COMP);
        );
        hkl_info.ctr_size=1;
        break;
      }
      sum_FT2 += hkl_info.ctr_FT2[i];
      if (hkl_info.ctr_k[i] < min_k) min_k = hkl_info.ctr_k[i];
      if (max_k < hkl_info.ctr_k[i]) max_k = hkl_info.ctr_k[i];
    }
    if (hkl_info.ctr_size>1 && min_k*max_k >= 0) { // k vector not around 0 ?
      MPI_MASTER(
        printf("Single_crystal: %s: Warning: surf_k = [%g:%g]. Must be around 0. Ignoring CTR.\n", min_k, max_k, NAME_CURRENT_COMP);
      );
      hkl_info.ctr_size=1;
    }
    if (hkl_info.ctr_size>1 && sum_FT2 <= 0) { // no positive FT2 value ?
      MPI_MASTER(
        printf("Single_crystal: %s: Warning: sum(surf_FT2) <= 0. Must be > 0. Ignoring CTR.\n", min_k, max_k, NAME_CURRENT_COMP);
      );
      hkl_info.ctr_size=1;
    }
    if (hkl_info.ctr_size>1 && sum_FT2) { // all is fine. Normalize CTR.
      NORM(hkl_info.ctr_dir[0],hkl_info.ctr_dir[1],hkl_info.ctr_dir[2]); // normal vector |u|=1
      // normalize surf_FT2 (distribute Bragg scattering along all satellites)
      for (i=0; i<hkl_info.ctr_size; hkl_info.ctr_FT2[i++] /= sum_FT2);
    }
  } // if (CTR)

  /* ought to be cleaned up as mosaic_AB now is a proper vector/array and not a define */
  double* mosaic_ABin = mosaic_AB;
  /* Read in structure factors, and do some pre-calculations. */
  if (!read_hkl_data(reflections, &hkl_info, &hkl_list, mosaic, mosaic_a, mosaic_b, mosaic_c, mosaic_ABin,
    hkl_info.ctr_size,hkl_info.ctr_k,hkl_info.ctr_FT2,hkl_info.ctr_dir)) {
    MPI_MASTER(
    printf("Single_crystal: %s: ERROR: can not set crystal structure. Check file and mosaic. Aborting.\n", NAME_CURRENT_COMP);
    );
    exit(-1);
  }

  if (hkl_info.sigma_i<0) hkl_info.sigma_i=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 uniform incoherent elastic scattering only sigma=%g.\n",
      NAME_CURRENT_COMP, hkl_info.sigma_i);
      
  if (!sx_read_abs_data(material_datafile, &abs_info, hkl_info )) abs_info.muc=0;

  MPI_MASTER(
  printf("Single_crystal: %s: Vc=%g [Angs] sigma_inc=%g [barn] reflections=%s\n",
      NAME_CURRENT_COMP, hkl_info.V0, hkl_info.sigma_i,
      reflections && strlen(reflections) ? reflections : "NULL");
  );

  if (powder && !(order==1)) {
    MPI_MASTER(
    fprintf(stderr,"Single_crystal: %s: powder mode means implicit choice of no multiple scattering!\n"
            "WARNING setting order=1\n", NAME_CURRENT_COMP);
    );
    order=1;
  }
  if (powder && (RX || RY)) {
    exit(fprintf(stderr,"Single_crystal: %s: powder mode can not be used together with crystal curvature.\n", NAME_CURRENT_COMP));
  }
  MPI_MASTER(
  printf("Direct space lattice orientation:\n");
  printf("  a = [%g %g %g]\n", hkl_info.m_ax, hkl_info.m_ay, hkl_info.m_az);
  printf("  b = [%g %g %g]\n", hkl_info.m_bx, hkl_info.m_by, hkl_info.m_bz);
  printf("  c = [%g %g %g]\n", hkl_info.m_cx, hkl_info.m_cy, hkl_info.m_cz);
  printf("Reciprocal space lattice orientation:\n");
  printf("  a* = [%g %g %g]\n", hkl_info.asx, hkl_info.asy, hkl_info.asz);
  printf("  b* = [%g %g %g]\n", hkl_info.bsx, hkl_info.bsy, hkl_info.bsz);
  printf("  c* = [%g %g %g]\n", hkl_info.csx, hkl_info.csy, hkl_info.csz);
  );

%}

TRACE
%{
  double l1, l2=0;              /* Entry and exit lengths in sample */
  struct hkl_data *L;           /* Structure factor list */
  int i;                        /* Index into structure factor list */
#ifndef OPENACC
  struct tau_data *T;           /* List of reflections close to Ewald sphere */
#else
  struct tau_data T[MCSX_REFL_SLIST_SIZE];
#endif
  int tau_count;                /* Number 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 k;                     /* photon wave vector */
  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 */
  double V0;                    /* Volume of unit cell */
  double l_full;                /* photon path length for transmission */
  double l;                     /* Path length to scattering event */
  double abs_xsect, abs_xlen;   /* Absorption cross section and length */
  double inc_xsect, inc_xlen;   /* Incoherent scattering 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, sum;           /* Temporaries */

  double p_trans;               /* Transmission probability */
  double mc_trans, mc_interact; /* Transmission, interaction MC choices */
  int    intersect=0;

  double curv_xangle;
  double curv_yangle;

  double _kx;
  double _ky;
  double _kz;

  char   type;      /* type of last event: t=transmit,c=coherent or i=incoherent */
  int    itype;     /* type of last event: t=1,c=2 or i=3 */

  #ifdef OPENACC
  #ifdef USE_OFF
  off_struct thread_offdata = offdata;
  #endif
  #else
  #define thread_offdata offdata
  #endif
  
  /* Intersection photon trajectory / sample (sample surface) */
  if (hkl_info.shape == 0)
    intersect = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius, yheight);
  else if (hkl_info.shape == 1)
    intersect = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
  else if (hkl_info.shape == 2)
    intersect = sphere_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius);
  #ifdef USE_OFF
  else if (hkl_info.shape == 3)
      intersect = off_x_intersect(&l1, &l2, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata );
  #endif
  if (l2 < 0) intersect=0;  /* we passed sample volume already */

  if(intersect)
  {                         /* photon intersects crystal */
    if(l1 > 0)
      PROP_DL(l1);          /* Move to crystal surface if not inside */
    ki  = sqrt(kx*kx + ky*ky + kz*kz);
    event_counter = 0;

    /*absorption cross-section */
    abs_xsect=0;
    if (abs_info.muc) {
      double mu=Table_Value(abs_info.table, ki*K2E, abs_info.muc);
      abs_xsect = (abs_info.atomic_weight/(6.02214076*10E23))*mu; //this is in cm^2
      //atomic weight = [g . mol^-1]
      //Avogadro number = [mol^-1]
      //mu = [cm squared . g⁻1]
      //therefore convert to barns if barns set to 1, otherwise to fm^2
      if (barns) { //convert cm^2 into barns
        abs_xsect *= 10E24; 
      } else {  //convert cm^2 into fm^2, else we assume fm^2
        abs_xsect *= 10E26;
      }
    }
    V0= hkl_info.V0;
    abs_xlen  = abs_xsect/V0;

    /*look into Compton scattering*/
    inc_xsect = hkl_info.sigma_i;
    inc_xlen  = inc_xsect/V0;

    if (barns) {
      /*If cross sections are given in barns, we need a scaling factor of 100 
       to get scattering lengths in m, since V0 is assumed to be in AA*/
    abs_xlen *= 100; inc_xlen *= 100;
    } /* else assume fm^2 */
    L = hkl_list;
    
    type = '\0';
    itype = 0;
    
#ifndef OPENACC
    T = tau_list;
    hkl_info.type = type;
#endif
    do {  /* Loop over multiple scattering events */
      /* Angles for powder randomization */
      double Alpha, Beta, Gamma;

      if (hkl_info.shape == 0)
        intersect = cylinder_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius, yheight);
      else if (hkl_info.shape == 1)
        intersect = box_intersect(&l1, &l2, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
      else if (hkl_info.shape == 2)
        intersect = sphere_intersect(&l1, &l2, x, y, z, kx, ky, kz, radius);
      #ifdef USE_OFF
      else if (hkl_info.shape == 3)
	      intersect = off_x_intersect(&l1, &l2, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata );
      #endif
      if(!intersect || l2 < -1e-7 || l1 > 1e-7)
      {
        /* photon is leaving the sample */
        hkl_info.flag_warning++;
        if (hkl_info.flag_warning < 10)
#ifndef OPENACC
          MPI_MASTER(
          fprintf(stderr,
                "Single_crystal: %s: Warning: photon has unexpectedly left the crystal!\n"
                "                l1=%g l2=%g x=%g y=%g z=%g kx=%g ky=%g kz=%g order=%i\n",
                NAME_CURRENT_COMP, l1, l2, x, y, z, kx, ky, kz, event_counter);
          );
#endif
        break;
      }

      l_full = l2;
	  
      if (   (order && !(extra_order) && event_counter >= order)
  		  || (order && extra_order && event_counter >= order + extra_order)) {
        // Exit due to truncated order, weight with relevant cross-sections to distance l_full
        p*=exp(-abs_xlen*l_full);
        intersect=0; 
        break;
      }

      /* (1). Compute incoming wave vector ki */
      if (powder) { /* orientation of crystallite is random */
        Alpha = randpm1()*PI*powder;
        Beta  = randpm1()*PI/2;
        Gamma = randpm1()*PI;
        randrotate(&kx, &ky, &kz, Alpha, Beta, Gamma);
      }

      /* ------------------------------------------------------------------------- */
      /* lattice curvature option: rotate photon velocity */
      curv_xangle = 0;
      curv_yangle = 0;

      _kx = kx;
      _ky = ky;
      _kz = kz;

      if(RY) { /* rotate v around x axis based on y pos, for vertical focus */
          curv_yangle = atan2(y, RY);
          vec_rotate_2d(&ky,&kz, curv_yangle);
          vec_rotate_2d(&Ey,&Ez, curv_yangle);

          /*changing y,z actually curves the crystal, not only the planes*/
          /*comment out if only curvature of the lattice planes is needed*/
          vec_rotate_2d(&y,&z, curv_yangle);
      }
      if(RX) { /* rotate v around y axis based on x pos, for horizontal focus */
          curv_xangle = atan2(x, RX);
          vec_rotate_2d(&kx,&kz, curv_xangle);
          vec_rotate_2d(&Ex,&Ez, curv_xangle);

          /*changing x,z actually curves the crystal, not only the planes*/
          /*comment out if only curvature of the lattice planes is needed*/
          vec_rotate_2d(&x,&z, curv_xangle);
      }

      kix = kx;
      kiy = ky;
      kiz = kz;
      kx = _kx;
      ky = _ky;
      kz = _kz;
      /* ------------------------------------------------------------------------- */

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

      double coh_xsect = 0, coh_refl = 0;
	  
	  // Condition to skip calculation of coherent cross section when, needed for extra_order feature	  
	  if (order==0 || extra_order==0 || event_counter < order) {	  
      /* in case we use 'SPLIT' then consecutive photons can be identical when entering here
         and we may skip the hkl_search call */
#ifndef OPENACC
      if (order==1 && fabs(kix - hkl_info.kix) < deltak
        && fabs(kiy - hkl_info.kiy) < deltak
        && fabs(kiz - hkl_info.kiz) < deltak) {
        hkl_info.nb_reuses++;

        /* Restore in case of matching event (e.g. SPLIT) */
        coh_refl  = hkl_info.coh_refl;
        coh_xsect = hkl_info.coh_xsect;
        tau_count = hkl_info.tau_count;

      } else {
#endif
        /* Max possible tau for this ki with 5*sigma delta-d/d cutoff. */
        double tau_max   = 2*ki/(1 - 5*hkl_info.m_delta_d_d);

        /* call hkl_search */
        
        tau_count = hkl_search(L, T, hkl_info.count, hkl_info.V0, 
                kix, kiy, kiz, tau_max,
                &coh_refl, &coh_xsect);

        /* store ki so that we can check for further SPLIT iterations */
#ifndef OPENACC
        if (tau_count>hkl_info.max_tau_count){
          hkl_info.max_tau_count=tau_count;
        }
        if (event_counter == 0 ) { /* only for incoming photon */
          hkl_info.kix = kix;
          hkl_info.kiy = kiy;
          hkl_info.kiz = kiz;

	        /* Store for potential re-use (e.g. SPLIT) */
	        hkl_info.coh_refl  = coh_refl;
	        hkl_info.coh_xsect = coh_xsect;
	        hkl_info.tau_count = tau_count;
	        hkl_info.nb_refl += tau_count;
	        hkl_info.nb_refl_count++;
	      }
      }
#endif
      } else {
        // When extra_order used, disable coherent scattering after order reached, but continue
  	    // Set coherent cross section to zero to ignore coherent part
        coh_refl = 0;
        coh_xsect = 0;
        tau_count = 0;
      }
	  
      /* (3). Probabilities of the different possible interactions. */
      /* 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 = coh_xsect/V0;
      if (barns) {
        coh_xlen *= 100;
      } /* else assume fm^2 */
      tot_xlen = abs_xlen + inc_xlen + coh_xlen;

      if(tot_xlen <= 0){
        ABSORB; // Should we really absorb here? If "nothing" can happen we perhaps ought to "pass" instead? 
      }
      
      /* (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;
      }
      mc_interact = 1 - mc_trans;
      if(mc_trans > 0 && (mc_trans >= 1 || rand01() < mc_trans))  /* Transmit */
      {
        p *= p_trans/mc_trans;
        intersect=0;
        if (powder) { /* orientation of crystallite is longer random */
          randderotate(&kx, &ky, &kz, Alpha, Beta, Gamma);
        }

        type = 't';
        if (!itype) itype = 1;
        #ifndef OPENACC
        hkl_info.type = type;
        #endif

        break; 
        /* This break means that we are leaving the while-loop, exiting the
	   crystal by "tunneling". */
      }

      /* Scattering "proper", i.e. coh or incoh */
      if(mc_interact <= 0)        /* Protect against rounding errors */
        { intersect=0;
          if (powder) { /* orientation of crystallite is no longer random */
            randderotate(&kx, &ky, &kz, Alpha, Beta, Gamma);
          }
          break;
        }

      /* First-pass considerations: */
      if (!event_counter) p *= fabs(1 - p_trans)/mc_interact;
      /* Select a point at which to scatter the photon, 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;
 
      /* Propagate to scattering point */
      PROP_DL(l);
      event_counter++;

      /* (4). Account for the probability of sigma_abs */
      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, kx, ky, kz, 0);
        kx = kix; /* ki vector is used as tmp var with norm v */
        ky = kiy;
        kz = kiz; /* Go for next scattering event */
	
        type = 'i';
        if (!itype) itype = 2;
#ifndef OPENACC
        hkl_info.type = type;
#endif
      } else {
        /* 7. Coherent scattering. Select reciprocal lattice point. */
        if(coh_refl <= 0){
          ABSORB;
        }
        sum = 0;
        j = hkl_select(T, tau_count, coh_refl, &sum,_particle);
        if(j >= tau_count)
        {
#ifndef OPENACC
          if (hkl_info.flag_warning < 10)
            MPI_MASTER(
            fprintf(stderr, "Single_crystal: %s: Warning: failed tau search "
              "(sum=%g, coh_refl=%g, j=%i, tau_count=%i). Using last reflection.\n", NAME_CURRENT_COMP, sum, coh_refl, j , tau_count);
            );
#endif
          hkl_info.flag_warning++;
          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 photon weight (see manual for explanation). */
        double pmul = T[j].xsect*coh_refl/(coh_xsect*T[j].refl);
        if (!isnan(pmul)) p *= pmul;
        kx = L[i].u1x*kfx + L[i].u2x*kfy + L[i].u3x*kfz;
        ky = L[i].u1y*kfx + L[i].u2y*kfy + L[i].u3y*kfz;
        kz = L[i].u1z*kfx + L[i].u2z*kfy + L[i].u3z*kfz;
        
        /* add thin layer k-broadening, only with box/cylinder shape */
        if (hkl_info.shape == 1) { /* box */
          double thickness_threshold = 1000*2*PI/ki/1e10; // 1000*lambda in Angs, /1e10 -> m
          if (xwidth  < thickness_threshold) kx += 2*PI/xwidth /1e10*randnorm();
          if (yheight < thickness_threshold) ky += 2*PI/yheight/1e10*randnorm();
          if (zdepth  < thickness_threshold) kz += 2*PI/zdepth /1e10*randnorm();
        } else if (hkl_info.shape == 0) { /* cylinder */
          double thickness_threshold = 1000*2*PI/ki/1e10;
          if (yheight < thickness_threshold) ky += 2*PI/yheight/1e10*randnorm();
        }
        
	
        type = 'c';
        if (!itype) itype = 3;
#ifndef OPENACC
        hkl_info.type = type;
        hkl_info.h    = L[i].h;
        hkl_info.k    = L[i].k;
        hkl_info.l    = L[i].l;
#endif

      }
      /* ------------------------------------------------------------------------- */
      /* lattice curvature option: rotate back photon velocity */
      if(RX) {
          vec_rotate_2d(&kx,&kz, -curv_xangle);
          vec_rotate_2d(&Ex,&Ez, -curv_xangle);

          /*changing x,z actually curves the crystal, not only the planes*/
          /*comment out if only curvature of the lattice planes is needed*/
          vec_rotate_2d(&x,&z, -curv_xangle);
      }
      if(RY) {
          vec_rotate_2d(&ky,&kz, -curv_yangle);
          vec_rotate_2d(&Ey,&Ez, -curv_yangle);

          /*changing y,z actually curves the crystal, not only the planes*/
          /*comment out if only curvature of the lattice planes is needed*/
          vec_rotate_2d(&y,&z, -curv_yangle);
      }
      /* ------------------------------------------------------------------------- */
      SCATTER;
      if (powder) { /* orientation of crystallite is no longer random */
        randderotate(&kx, &ky, &kz, Alpha, Beta, Gamma);
      }
      /* Repeat loop for next scattering event. */
    } while (intersect); /* end do (intersect) (multiple scattering loop) */
  } /* if intersect */
%}

FINALLY
%{
  MPI_MASTER(
  if (hkl_info.flag_warning)
    fprintf(stderr, "Single_crystal: %s: Error message was repeated %i times with absorbed photons/illegal tau search.\n",
      NAME_CURRENT_COMP, hkl_info.flag_warning);

  /* in case this instance is used in a SPLIT, we can recommend the
     optimal iteration value */
  if (hkl_info.max_tau_count>=MCSX_REFL_SLIST_SIZE){
    fprintf(stderr,"Single_crystal: %s: Warning: The reflection short list buffer was exhausted at least once. Please consider redefining MCSX_REFL_SLIST_SIZE > %d\n", NAME_CURRENT_COMP, MCSX_REFL_SLIST_SIZE);
  }

  if (hkl_info.nb_refl_count) {
    double split_iterations = (double)hkl_info.nb_reuses/hkl_info.nb_refl_count + 1;
    double split_optimal    = (double)hkl_info.nb_refl/hkl_info.nb_refl_count;
    if (split_optimal > split_iterations + 5) {
      printf("Single_crystal: %s: Info: you may improve the computation efficiency by using\n"
        "    SPLIT %i COMPONENT %s=Single_crystal(order=1, ...)\n"
        "  in the instrument description %s.\n",
        NAME_CURRENT_COMP, (int)split_optimal, NAME_CURRENT_COMP, instrument_source);
    }
  }
  );
%}

MCDISPLAY
%{
  if (hkl_info.shape == 0) {	/* cylinder */
    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);
  }
  else if (hkl_info.shape == 1) { 	/* box */
    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);
  }
  else if (hkl_info.shape == 2) {	/* sphere */
    circle("xy", 0,  0.0, 0, radius);
    circle("xz", 0,  0.0, 0, radius);
    circle("yz", 0,  0.0, 0, radius);
  }
  else if (hkl_info.shape == 3) {	/* OFF file */
    off_display(offdata);
  }
%}
END