/*****************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Fluorescence
*
* %Identification
* Written by: E. Farhi
* Date:       March 2022
* Origin:     Synchrotron SOLEIL
*
* Sample model handling absorption, fluorescence, Compton and Rayleigh scattering.
*
* %Description
* Sample that models many photon-matter interactions:
* - absorption          (photon excites an electron and creates a hole)
* - fluorescence        (excited electrons emit light while falling into lower states)
* - Compton scattering  (inelastic, transfer some energy to an electron,    incoherent)
* - Rayleigh scattering (elastic,   dipolar radiation of excitated electrons, coherent)
*
* Use the <b>FluoPowder</b> sample component to properly handle powder diffraction with fluorescence.
*
* The 'material' specification is given as a chemical formulae, e.g. "LaB6". It
* may also be given as a file name (CIF/LAU/LAZ/FullProf format) in which case 
* the formulae is guessed (but may be approximative). Atoms from Z=5 to Z=90 are
* handled.
*
* By setting the 'order' to 1, the absorption along the scattered path is handled.
* A higher 'order' will handle multiple scattering events, and final absorption.
* For instance, a value order>=2 handles e.g. fluorescence iterative cascades 
* in the material. Leaving 'order=0' handles the single scattering only.
*
* Example: Fluorescence(material="LaB6",
*  xwidth=0.001,yheight=0.001,zdepth=0.0001, p_interact=0.99,
*  target_index=1, focus_xw=0.0005, focus_yh=0.0005)
*
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or any other shape
*   box/plate:       xwidth x yheight x zdepth (thickness=0)
*   hollow box/plate:xwidth x yheight x zdepth and thickness>0
*   cylinder:        radius x yheight (thickness=0)
*   hollow cylinder: radius x yheight and thickness>0
*   sphere:          radius (yheight=0 thickness=0)
*   hollow sphere:   radius and thickness>0 (yheight=0)
*   any shape:       geometry=OFF file
*
*   The complex geometry option handles any closed non-convex polyhedra.
*   It computes the intersection points of the photon ray with the object
*   transparently, so that it can be used like a regular sample object.
*   It supports the OFF, PLY and NOFF file format but not COFF (colored faces).
*   Such files may be generated from XYZ data using:
*     qhull < coordinates.xyz Qx Qv Tv o > geomview.off
*   or
*     powercrust coordinates.xyz
*   and viewed with geomview or java -jar jroff.jar (see below).
*   The default size of the object depends of the OFF file data, but its
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* <b>Concentric components:</b>
* This component has the ability to contain other components when used in
* hollow cylinder geometry (namely sample environment, e.g. cryostat and
* furnace structure). Such component 'shells' should be split into input and
* output side surrounding the 'inside' components. First part must then use
* 'concentric=1' flag to enter the inside part. The component itself must be
* repeated to mark the end of the concentric zone. The number of concentric
* shells and number of components inside is not limited.
*
* COMPONENT F_in = Fluorescence(material="Al", concentric=1, ...)
* AT (0,0,0) RELATIVE sample_position
*
* COMPONENT something_inside ... // e.g. the sample itself or other materials
*
* COMPONENT F_out = COPY(F_in)(concentric=0)
* AT (0,0,0) RELATIVE sample_position
*
* <b>Enhancing computation efficiency:</b>
* * An important option to enhance statistics is to set 'p_interact' to, say,
* 30 percent (0.3) in order to force a fraction of the beam to scatter. This
* will result on a larger number of scattered events, retaining intensity.
*
* In addition, it may be desirable to define a 'target' for the fluorescence
* processes via e.g. the 'target_index' and the 'focus_xw / focus_yh' options.
* This target should e.g. be the SDD area.
*
* This sample component can advantageously benefit from the SPLIT feature, e.g.
* SPLIT COMPONENT sample = Fluorescence(...)
*
* The computation is made via the XRayLib (apt install libxrl-dev) https://github.com/tschoonj/xraylib.
*
* %Parameters
* material:  [str]    Chemical formulae, e.g. "LaB6", "Pb2SnO4". If may also be a CIF/LAZ/LAU file.
* weight:     [g/mol] Atomic/molecular weight of material.
* density:   [g/cm^3] Density of material. V_rho=density/weight/1e24*N_A.
* rho:         [AA-3] Density of scattering elements (nb atoms/unit cell V_0).
* packing_factor: [1] How dense is the material compared to bulk 0-1.
* radius:         [m] Outer radius of sample in (x,z) plane. cylinder/sphere.
* xwidth:         [m] Width for a box sample shape.
* yheight:        [m] Height of sample in vertical direction for box/cylinder shapes.
* zdepth:         [m] Depth for a box sample shape.
* thickness:      [m] Thickness of hollow sample Negative value extends the hollow volume outside of the box/cylinder.
* concentric:     [1] Indicate that this component has a hollow geometry and may contain other components. It should then be duplicated after the inside part (only for box, cylinder, sphere).
* 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.
* p_interact:     [1] Force a given fraction of the beam to scatter, keeping intensity right, to enhance small signals (-1 inactivate).
* focus_xw:       [m] Horiz. dimension of a rectangular area.
* focus_yh:       [m] Vert.  dimension of a rectangular area.
* focus_aw:     [deg] Horiz. angular dimension of a rectangular area.
* focus_ah:     [deg] Vert.  angular dimension of a rectangular area.
* focus_r:        [m] Radius of disk containing target. Use 0 for full space.
* target_index:   [1] Relative index of component to focus at, e.g. next is +1.
* target_x:       [m] Position of target to focus at, along X.
* target_y:       [m] Position of target to focus at, along Y.
* target_z:       [m] Position of target to focus at, along Z.
* flag_compton:   [1] When 0, the Compton  scattering is ignored.
* flag_rayleigh:  [1] When 0, the Rayleigh scattering is ignored.
* flag_lorentzian:[1] When 1, the line shapes are assumed to be Lorentzian, else Gaussian.
* flag_kissel:    [1] When 1 (slower), handle M-lines XRF from Kissel for Z>=52 Te (else only K and L-lines).
* flag_low_z:     [1] When 1, enhances low concentration atoms, else uses cross-sections weighting.
* order:          [1] Limit multiple fluorescence up to given order. Last iteration is absorption only.
*
* CALCULATED PARAMETERS:
* type: scattering event type 0=FLUORESCENCE fluorescence, 1=RAYLEIGH Rayleigh, 2=COMPTON Compton, 3=TRANSMISSION Transmit (absorbsion), 4=DIFFRACTION Diffraction (not handled here)
*
* %Link
* The XRayLib https://github.com/tschoonj/xraylib http://dx.doi.org/10.1016/j.sab.2011.09.011
* %Link
* Fluorescence https://en.wikipedia.org/wiki/Fluorescence
* %Link
* Rayleigh https://en.wikipedia.org/wiki/Rayleigh_scattering
* %Link
* Compton https://en.wikipedia.org/wiki/Compton_scattering
* %Link
* X-ray absorption edges http://skuld.bmsc.washington.edu/scatter/AS_periodic.html
* %Link
* X-ray fluorescence spectra http://www.xrfresearch.com/xrf-spectra/
* %Link
* X-ray edges and fluo lines https://physics.nist.gov/PhysRefData/XrayTrans/Html/search.html
*
* %E
***********************************************************/

DEFINE COMPONENT Fluorescence
SETTING PARAMETERS(
  string geometry=0,
  radius=0, thickness=0,
  xwidth=0, yheight=0, zdepth=0,
  int concentric=0,
  string material="LaB6", packing_factor=0, rho=0, density=0, weight=0, 
  p_interact=0, 
  target_x = 0, target_y = 0, target_z = 0, focus_r = 0,
  focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, int target_index=0,
  int flag_compton=1, int flag_rayleigh=1, int flag_lorentzian=0,
  int flag_kissel=0,  int flag_low_z=0,
  int order=1)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
DEPENDENCY " @XRLFLAGS@ -DUSE_OFF "
NOACC

/* ========================================================================== */

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

  #ifndef FLUORESCENCE_DECL
  #define FLUORESCENCE_DECL

  // WARNING: mush NOT use rand01 calls here, as they depend on _particle

  #define XRAYLIB_LINES_MAX 383
  #define FLUORESCENCE 0  // Fluo
  #define RAYLEIGH     1  // Coherent
  #define COMPTON      2  // Incoherent
  #define DIFFRACTION  3  // diffraction
  #define TRANSMISSION 4

  #include <xraylib/xraylib.h>
  
/* See XrayLib code (c) T. Schoonjans
 * /usr/include/xraylib/xraylib-parser.h      for compoundData
 * /usr/include/xraylib/xraylib.h             for XS
 */

/* inspired from:
 * https://github.com/golosio/xrmc src/photon/photon.cpp (c) Bruno Golosio    
 */
  
  /* XRMC_CrossSections: Compute interaction cross sections in [barn/atom]
   * When kissel=1, M-lines are computed, but computation is 10x slower.
   * Computes xs[FLUORESCENCE,RAYLEIGH,COMPTON] that must have been allocated.
   * Return total cross section, given Z and E0:
   *   total_xs = XRMC_CrossSections(Z, E0, xs[3], flag_kissel);
   */
  double XRMC_CrossSections(int Z, double E0, double *xs, int kissel) {
    int    i_line;
    
    if (xs == NULL) return 0;
    for(i_line=0; i_line<=COMPTON; i_line++) xs[i_line]=0;
    
    // loop on possible fluorescence lines
    for (i_line=0; i_line<XRAYLIB_LINES_MAX; i_line++) { 
      // cumulative sum of the line cross sections
      xs[FLUORESCENCE] += kissel && Z >= 50 ?
        CSb_FluorLine_Kissel(Z, -i_line, E0, NULL) /* XRayLib lines are negative integers */
      : CSb_FluorLine(Z, -i_line, E0, NULL);
    }

    // coherent and incoherent cross sections
    xs[RAYLEIGH] = CSb_Rayl( Z, E0, NULL);
    xs[COMPTON]  = CSb_Compt(Z, E0, NULL);
    
    // total interaction cross section, should converge to CSb_Total(Z, E0, NULL)
    return xs[FLUORESCENCE] + xs[RAYLEIGH] + xs[COMPTON];
  } // XRMC_CrossSections

  /* XRMC_SelectFromDistribution: select a random element from a distribution
   *   index = XRMC_SelectFromDistribution(cum_sum[N], N);
   * index is returned within 0 and N-1
   * The x_arr must be a continuously increasing cumulated sum, which last element is the max.
   */
  int XRMC_SelectFromDistribution(double x_arr[], int N, double r)
  {
    if (x_arr == NULL || N <=0) return (0);
    double x = r*x_arr[N-1]; // the threshold cumulated value to search in distribution
    if (x<x_arr[0]) {    // x is smaller than lower limit
      return 0;
    }
    if (x>=x_arr[N-1]) { // x is greater or equal to upper limit
      return N-1;
    }
    int id=0, iu=N-1; // lower and upper index of the subarray to search
    while (iu-id>1) { // search until the size of the subarray to search is >1
      int im = (id + iu)/2; // use the midpoint for equal partition
      // decide which subarray to search
      if (x>=x_arr[im]) id=im; // change min index to search upper subarray
      else iu=im; // change max index to search lower subarray
    }

    return id;
  } // XRMC_SelectFromDistribution

  /* XRMC_SelectInteraction: select interaction type Fluo/Compton/Rayleigh
   * Return the interaction type from a random choice within cross sections 'xs'
   *   type = XRMC_SelectInteraction(xs[3], 3);
   * 'xs' is computed with XRMC_CrossSections. Can be unsorted.
   * type is e.g. one of FLUORESCENCE | RAYLEIGH | COMPTON
   */
  int XRMC_SelectInteraction(double *xs, int nb_processes, double r)
  {
    double sum_xs, *cum_xs;
    int    i;

    if (xs == NULL || nb_processes < 1) return 0;
    cum_xs = malloc((nb_processes+1)*sizeof(double));
    if (!cum_xs)          return 0;
    cum_xs[0]=sum_xs=0;
    for (i=0; i< nb_processes; i++) {
      sum_xs += xs[i];
      cum_xs[i+1]= sum_xs;
    }
    i = XRMC_SelectFromDistribution(cum_xs, nb_processes+1, r);
    free(cum_xs);
    return i;
  } // XRMC_SelectInteraction

  /* XRMC_SelectFluorescenceEnergy: select outgoing fluo photon energy, when incoming with 'E0'
   * When kissel=1, M-lines are included, but computation is 10x slower
   *   Ef = XRMC_SelectFluorescenceEnergy(Z, E0, &dE, kissel);
   */
  double XRMC_SelectFluorescenceEnergy(int Z, double E0, double *dE, int kissel, double r)
  {
    int i_line;
    double sum_xs, cum_xs_lines[XRAYLIB_LINES_MAX+1];

    // compute cumulated XS for all fluo lines
    cum_xs_lines[0] = sum_xs = 0;
    for (i_line=0; i_line<XRAYLIB_LINES_MAX; i_line++) { // loop on fluorescent lines
      double xs = kissel && Z >= 50 ?
        CSb_FluorLine_Kissel(Z, -i_line, E0, NULL) /* XRayLib lines are negative integers */
      : CSb_FluorLine(Z, -i_line, E0, NULL);
      // when a line is inactive: E=xs=0
      sum_xs += xs;
      cum_xs_lines[i_line+1] = sum_xs; // cumulative sum of their cross sections
    }
    // select randomly one of these lines
    i_line = XRMC_SelectFromDistribution(cum_xs_lines, XRAYLIB_LINES_MAX, r); // extract a line
    // get the K shell line width as approximation of fluorescence line width
    if (dE) *dE = AtomicLevelWidth(Z, K_SHELL, NULL); // keV, full-width in keV

    return LineEnergy(Z, -i_line, NULL); // fluorescent line energy
  } // XRMC_SelectFluorescenceEnergy

  // Function removing spaces from string
  char * removeSpacesFromStr(char *string)
  {
      // non_space_count to keep the frequency of non space characters
      int non_space_count = 0;
      if (string == NULL) return NULL;
   
      //Traverse a string and if it is non space character then, place it at index non_space_count
      for (int i = 0; string[i] != '\0'; i++)
      {
          if (isalpha(string[i]) || isdigit(string[i]))
          {
              string[non_space_count] = string[i];
              non_space_count++; //non_space_count incremented
          }    
      }
      
      //Finally placing final character at the string end
      string[non_space_count] = '\0';
      return string;
  }

  // ok = fluo_get_material(material, formula)
  // extracts material atoms from file header
  // the result is concatenated into 'formula'
  int fluo_get_material(char *filename, char *formula) {
  
    int  ret = 0;
    char Line[65535];
    int  flag_found_cif=0;
    if (filename == NULL || formula == NULL) return (0);
    FILE *file = Open_File(filename, "r", NULL);
    if (!file)
      exit(fprintf(stderr, "%s: ERROR: can not open file %s\n", 
      __FILE__, filename));
    
    // Read the file, and search tokens in rows
    formula[0]=0;
    while (fgets(Line, sizeof(Line), file) != NULL) {
      char *token = NULL;
      int   flag_exit=0;
      char *first_non_space=NULL;
      char *next_non_space =NULL;
      // CIF: _chemical_formula_structural 'chemical_formulae'
      // CIF: _chemical_formula_sum 'chemical_formulae'
      // LAZ/LAU: # ATOM <at> <trailing>
      // LAZ/LAU: # Atom <at> <trailing>
      // LAZ/LAU: # TITLE <at> <at> ... [ trailing...]
      // CFL: Title <chemical_formulae>
      // CFL: Atom <at> <trailing>
      
      // search for CIF token
      // single line: search " '\'" delimiter after CIF token, marks reading the formula
      if (!strncasecmp(Line, "_chemical_formula_structural", 28) 
       || !strncasecmp(Line, "_chemical_formula_sum", 21) 
       || !strncasecmp(Line, "_chemical_formula_moiety", 24)
       || flag_found_cif) {
        if  (flag_found_cif) { flag_found_cif=0; /* can not span on more that 2 lines */} 
        else flag_found_cif=1;
        // search for delimiter after the CIF token
        char *first_space_after_token=strpbrk(Line, " \'\n");
        if (first_space_after_token) {
          // search for the characters that may compose the formula
          first_non_space = strpbrk(first_space_after_token, "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789()[]");
          if (strchr(first_space_after_token,'?')) {
            // we have an invalid/unknown formula in the CIF. Skip CIF token/line
            flag_found_cif=0;
            continue;
          }
        }
        if (first_non_space) 
          next_non_space  = strpbrk(first_non_space+1, "\'\n\r"); // position of formula end
        if (next_non_space) { 
          flag_exit = 1;
          token = first_non_space;
        }
      } else if (!strncasecmp(Line, "# TITLE", 7) && strchr(Line, '['))
        token = Line+7;
      else if (!strncasecmp(Line, "# Atom", 6))
        token = Line+6;
      else if (!strncasecmp(Line, "Atom", 4))
        token = Line+4;
      else if (!strncasecmp(Line, "Title", 5))
        token = Line+7;
      if (!token) continue;
      if (!strncasecmp(Line, "# TITLE", 7)) {
        first_non_space = Line+7;
        next_non_space  = strchr(Line+7, '[');
        if (next_non_space) flag_exit = 1;
      }
      if (!first_non_space) first_non_space = strtok(token, " \t\n");
      if (!first_non_space) continue;
      if (!next_non_space)  next_non_space  = strtok(NULL,  " \t\n"); // end of formulae
      if (!next_non_space)  next_non_space = Line+strlen(Line);
      // remove spaces
      strncat(formula, first_non_space, next_non_space-first_non_space);
      ret++;
      if (flag_exit) break;
    }
    fclose(file);
    removeSpacesFromStr(formula);
    return(ret);
  } // fluo_get_material
  
  #endif // FLUORESCENCE_DECL
%}

/* ========================================================================== */

DECLARE %{
  struct   compoundData *compound;
  DArray1d cum_massFractions;
  DArray1d cum_xs_fluo;
  DArray1d cum_xs_Compton;
  DArray1d cum_xs_Rayleigh;
  DArray1d xs_fluo;
  DArray1d xs_Compton;
  DArray1d xs_Rayleigh;
  int  shape;
  off_struct offdata;
  int  n_fluo;
  int  n_Compton;
  int  n_Rayleigh;
  double p_fluo;
  double p_Compton;
  double p_Rayleigh;
  
  // cached values for SPLIT
  int    reuse_nb;
  int    reuse_intersect;
  double reuse_ki;
  double reuse_l0;
  double reuse_l1;
  double reuse_l2;
  double reuse_l3;
  double reuse_sigma_barn;
  DArray1d reuse_cum_xs_fluo;
  DArray1d reuse_cum_xs_Compton;
  DArray1d reuse_cum_xs_Rayleigh;
  DArray1d reuse_xs_fluo;
  DArray1d reuse_xs_Compton;
  DArray1d reuse_xs_Rayleigh;
%}

INITIALIZE %{

  /* energies en [keV], angles in [radians], XRL CSb cross sections are in [barn/atom] */
  xrl_error *error = NULL;
  int        i;
  
  XRayInit();
  
  shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape  */
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    #ifndef USE_OFF
    fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n");
    exit(-1);
    #else
    if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
      shape=3; thickness=0; concentric=0;
    }
    #endif
  }
  else if (xwidth && yheight && zdepth)  shape=1; /* box */
  else if (radius > 0 &&  yheight)       shape=0; /* cylinder */
  else if (radius > 0 && !yheight)       shape=2; /* sphere */

  if (shape < 0)
    exit(fprintf(stderr,"Fluorescence: %s: sample has invalid dimensions.\n"
                        "ERROR       Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
  
  if (!material || !strlen(material) || !strcmp(material, "NULL") || !strcmp(material, "0")) 
    exit(fprintf(stderr, "Fluorescence: ERROR: %s: Null material specification.\n", NAME_CURRENT_COMP));

  if (order > 1 && flag_lorentzian)
    fprintf(stderr, "Fluorescence: WARNING: %s: Using flag_lorentzian=1 Lorentzian line-shape with order>1 may create artifacts.\n", NAME_CURRENT_COMP);
    
  // test if the material is given as a file
  char path[1024];
  char formula[65536];
  formula[0]='\0';
  FILE *file = Open_File(material, "r", path);
  if (file != NULL) {
    fclose(file);
    // open the material structure file (laz/lau/cif...)
    // search (case sensitive) along lines
    if (!fluo_get_material(material, formula))
      exit(fprintf(stderr, "ERROR: %s: file %s does not contain material formulae.\n", NAME_CURRENT_COMP, material));
    
    MPI_MASTER(
    fprintf(stderr, "Fluorescence: INFO: %s: found material %s from file %s\n", NAME_CURRENT_COMP, formula, material);
    );
    strcpy(material, formula);
    // CIF: _chemical_formula_structural 'chemical_formulae'
    // CIF: _chemical_formula_sum 'chemical_formulae'
    // LAZ/LAU: # ATOM <at> <trailing>
    // LAZ/LAU: # Atom <at> <trailing>
    // LAZ/LAU: # TITLE <at> <at> ... [ trailing...]
    // CFL: Title <chemical_formulae>
    // CFL: Atom <at> <trailing>
    
  }
    
  compound = CompoundParser(material, &error); /* XRayLib */
  if (error != NULL) {
    exit(fprintf(stderr, "ERROR: %s: Invalid material %s: %s\n",
      NAME_CURRENT_COMP, material, error->message));
  }
  xrl_error_free(error);
  
  /* compute total density for raw material and display information ========= */  
  if (weight <= 0) weight = compound->molarMass; /* g/mol */
  MPI_MASTER(
    printf("%s: Material %s mass fractions:\n", 
      NAME_CURRENT_COMP, material);
  )
  
  double mat_density       = 0; /* g/cm3 */
  double sum_massFractions = 0;
  cum_massFractions        = create_darr1d(compound->nElements+1);
  cum_xs_fluo              = create_darr1d(compound->nElements+1);
  cum_xs_Compton           = create_darr1d(compound->nElements+1);
  cum_xs_Rayleigh          = create_darr1d(compound->nElements+1);
  xs_fluo                  = create_darr1d(compound->nElements+1);
  xs_Compton               = create_darr1d(compound->nElements+1);
  xs_Rayleigh              = create_darr1d(compound->nElements+1);
  cum_massFractions[0]     = 0;
  
  /* print material information, and check for elements */
  for (i=0; i< compound->nElements; i++) {
    int    Z      = compound->Elements[i];
    error = NULL;
    double Z_dens = ElementDensity(Z, &error);
    if (error != NULL)
      exit(fprintf(stderr, "ERROR: %s: Z=%i %s\n", NAME_CURRENT_COMP, Z, error->message));
    mat_density           += compound->massFractions[i]*Z_dens;
    sum_massFractions     += compound->massFractions[i];
    cum_massFractions[i+1] = sum_massFractions;
    MPI_MASTER(
      printf("  | %6.2g %%: Z=%3i %3s %8.3g [g/mol] %8.3g [g/cm3]\n",
        compound->massFractions[i]*100, Z, AtomicNumberToSymbol(Z,NULL), AtomicWeight(Z, NULL),
        Z_dens);
    )
  }
  xrl_error_free(error);
  if (density <= 0)        density        = mat_density;   /* g/cm3 */
  if (packing_factor <= 0) packing_factor = density/mat_density;
  
  /* 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 */
  if (!rho) rho = density/weight/1e24*NA; // atom density [at/Angs-3]
  MPI_MASTER(
    printf("%s: Material %s M=%g [g/mol] density=%g [g/cm3] rho=%g [at/Angs-3]",
      NAME_CURRENT_COMP, material, weight, density, rho);
    if (fabs(packing_factor-1) > 1e-2)
      printf(" packing_factor=%g", packing_factor);
    printf("\n");
  )
  if (0 < packing_factor && packing_factor < 1) rho *= packing_factor;
  
  /* target for scattering ================================================== */
  if (!target_index && !target_x && !target_y && !target_z) target_index=1;
  if (target_index)
  {
    Coords ToTarget;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
    coords_get(ToTarget, &target_x, &target_y, &target_z);
  }
  if (!(target_x || target_y || target_z)) {
    MPI_MASTER(
    printf("Fluorescence: %s: The target is not defined. Using 4PI.\n",
      NAME_CURRENT_COMP);
    );
  }
  
  n_fluo = n_Compton = n_Rayleigh = 0;
  p_fluo = p_Compton = p_Rayleigh = 0;
  
  // cached variables set to 0 (for SPLIT)
  reuse_nb=0;
  reuse_ki=reuse_l0=reuse_l1=reuse_l2=reuse_l3=reuse_sigma_barn=0;
  reuse_cum_xs_fluo              = create_darr1d(compound->nElements+1);
  reuse_cum_xs_Compton           = create_darr1d(compound->nElements+1);
  reuse_cum_xs_Rayleigh          = create_darr1d(compound->nElements+1);
  reuse_xs_fluo                  = create_darr1d(compound->nElements+1);
  reuse_xs_Compton               = create_darr1d(compound->nElements+1);
  reuse_xs_Rayleigh              = create_darr1d(compound->nElements+1);
  reuse_cum_xs_fluo[0]=reuse_cum_xs_Compton[0]=reuse_cum_xs_Rayleigh[0]=0;
%}

TRACE %{

int    intersect=0;       /* flag to continue/stop */
int    reuse    =0;       /* flag raised when reusing (SPLIT) */
int    type     =-1;      /* type of interaction */
double l0,  l1,  l2,  l3; /* distances for intersections */
double dl0, dl1, dl2, dl; /* distance intervals */
int    flag_concentric = 0;
int    flag_ishollow   = 0;
int    event_counter   = 0;         /* scattering event counter (multiple fluorescence) */
int    force_transmit  = 0;         /* Flag to handle cross-section weighting in case of finite order */
double sigma_barn=0, xs[5];         /* cross sections [barn/atom] fluo/Compton/Rayleigh/diff/transmit */
double aim_x=0, aim_y=0, aim_z=1;   /* Position of target relative to scattering point */


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

double k, E;
double ki,pi;

/* Store Initial photon state */
k  = ki = sqrt(kx*kx+ky*ky+kz*kz); //  Angs-1
E  = K2E*k; // keV
pi = p;     // used to test for multiple fluo weighting and order cutoff

do { /* while (intersect) Loop over multiple scattering events */

  // test for a SPLIT event (same particle comes in)
  l0=l1=l2=l3=dl0=dl1=dl2=0;
  if (!event_counter && fabs(reuse_ki - ki) < 1e-15) {
    reuse     = 1;
    reuse_nb++;
    // use cached values and skip actual computation
    intersect = reuse_intersect;
    l0        = reuse_l0;
    l1        = reuse_l1;
    l2        = reuse_l2;
    l3        = reuse_l3;

  } else {
    // we have a different event: compute intersection lengths
    
    /* ========================================================================== */
    /*                                   GEOMETRY                                 */
    /* ========================================================================== */

    /* Intersection photon trajectory / sample (sample surface) */
    if (thickness >= 0) {
      if (shape==0)
        intersect=cylinder_intersect(&l0,&l3, x,y,z,kx,ky,kz, radius,yheight);
      else if (shape==1)
        intersect=box_intersect     (&l0,&l3, x,y,z,kx,ky,kz, xwidth,yheight,zdepth);
      else if (shape==2)
        intersect=sphere_intersect  (&l0,&l3, x,y,z,kx,ky,kz, radius);
      #ifdef USE_OFF
      else if (shape == 3)
        intersect=off_x_intersect(&l0, &l3, NULL, NULL, x, y, z, kx,ky,kz, thread_offdata );
      #endif
    } else {
      if (shape==0)
        intersect=cylinder_intersect(&l0,&l3, x,y,z,kx,ky,kz, radius-thickness,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight);
      else if (shape==1)
        intersect=box_intersect     (&l0,&l3, x,y,z,kx,ky,kz,
          xwidth-2*thickness > 0 ?  xwidth-2*thickness : xwidth,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight,
          zdepth-2*thickness > 0 ?  zdepth-2*thickness : zdepth);
      else if (shape==2)
        intersect=sphere_intersect  (&l0,&l3, x,y,z,kx,ky,kz, radius-thickness);
      #ifdef USE_OFF
      else if (shape == 3)
        intersect=off_x_intersect(&l0, &l3, NULL, NULL, x, y, z, kx,ky,kz, thread_offdata );
      #endif
    }


    /* Computing the intermediate lengths */
    if (intersect && p_interact >= 0) {
      flag_ishollow = 0;
      if (thickness > 0) {
        if (shape==0 && cylinder_intersect(&l1,&l2, x,y,z,kx,ky,kz, radius-thickness,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight))
          flag_ishollow=1;
        else if (shape==2 && sphere_intersect   (&l1,&l2, x,y,z,kx,ky,kz, radius-thickness))
          flag_ishollow=1;
        else if (shape==1 && box_intersect(&l1,&l2, x,y,z,kx,ky,kz,
          xwidth-2*thickness > 0 ? xwidth-2*thickness : xwidth,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight,
          zdepth-2*thickness > 0 ? zdepth-2*thickness : zdepth))
          flag_ishollow=1;
      } else if (thickness<0) {
        if (shape==0 && cylinder_intersect(&l1,&l2, x,y,z,kx,ky,kz, radius,yheight))
          flag_ishollow=1;
        else if (shape==2 && sphere_intersect   (&l1,&l2, x,y,z,kx,ky,kz, radius))
          flag_ishollow=1;
        else if (shape==1 && box_intersect(&l1,&l2, x,y,z,kx,ky,kz, xwidth, yheight, zdepth))
          flag_ishollow=1;
      }
      if (l3 > 1000) l3=0;              // we passed the sample, far away
      if (!flag_ishollow) l1 = l2 = l3; /* no empty space inside */
    } /* if intersect */


    // store values for potential next SPLIT
    if (!event_counter) {
      reuse_intersect = intersect;
      reuse_l0        = l0;
      reuse_l1        = l1;
      reuse_l2        = l2;
      reuse_l3        = l3;
      reuse_ki        = ki;
    }
    reuse           = 0;
  } // if !reuse (SPLIT)

  if (intersect) { /* the photon hits the sample */

    if (l0 > 0) {  /* we are before the sample */
      PROP_DL(l0); /* propagates photon to the entry of the sample */
    } else if (l1 > 0 && l1 > l0) { /* we are inside first part of the sample */
      /* no propagation, stay inside */
    } else if (l2 > 0 && l2 > l1) { /* we are in the hole */
      PROP_DL(l2); /* propagate to inner surface of 2nd part of sample */
    } else if (l3 > 0 && l3 > l2) { /* we are in the 2nd part of sample */
      /* no propagation, stay inside */
    }

    dl0=l1-(l0 > 0 ? l0 : 0); /* Distance in first part of hollow/cylinder/box */
    dl1=l2-(l1 > 0 ? l1 : 0); /* Distance in hole */
    dl2=l3-(l2 > 0 ? l2 : 0); /* Distance in 2nd part of hollow cylinder */

    if (dl0 < 0) dl0 = 0;
    if (dl1 < 0) dl1 = 0;
    if (dl2 < 0) dl2 = 0;

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

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

    if (!dl0 && !dl2) {
      intersect = 0; /* the sample was passed entirely */
    }
  } // if intersect (geometry)

    
  /* ========================================================================== */
  /*                             INTERACTION PROCESS                            */
  /* ========================================================================== */
  xs[FLUORESCENCE]=xs[COMPTON]=xs[RAYLEIGH]=xs[TRANSMISSION]=xs[DIFFRACTION]=sigma_barn=0;
  cum_xs_fluo[0] = cum_xs_Compton[0] = cum_xs_Rayleigh[0] = 0;
  
  if (intersect) {
    double my_s;
    int    i_Z,i;
    int    flag=0;
    double d_path, p_trans, p_scatt, mc_trans, mc_scatt;
    
    /* actual fluorescence calculation */
    
    /* compute total scattering cross section for incoming photon energy Ei */
    /* compute each contribution XS */
    if (!reuse) {
      for (i_Z=0; i_Z< compound->nElements; i_Z++) {
        int    Z   = compound->Elements[i_Z];
        double frac= compound->massFractions[i_Z];
        double xs_Z[3];

        // get Fluorescence xs
        XRMC_CrossSections(Z, E, xs_Z, flag_kissel); // [barn/atom]

        // local XS for atom (Z)
        xs_fluo[i_Z+1]     = frac*xs_Z[FLUORESCENCE];
        xs_Compton[i_Z+1]  = frac*xs_Z[COMPTON];
        xs_Rayleigh[i_Z+1] = frac*xs_Z[RAYLEIGH];
        // selection on cumulated distribution xs*frac
        cum_xs_fluo[i_Z+1]     = cum_xs_fluo[i_Z]    +xs_fluo[i_Z+1];
        cum_xs_Compton[i_Z+1]  = cum_xs_Compton[i_Z] +xs_Compton[i_Z+1];
        cum_xs_Rayleigh[i_Z+1] = cum_xs_Rayleigh[i_Z]+xs_Rayleigh[i_Z+1];
        for (i=0; i<3; i++) { xs[i] += frac*xs_Z[i]; }
      } // for Z in compound
      for (i=0; i<3; i++) sigma_barn += xs[i]; // total XS
      
      // store values into cache for SPLIT
      if (!event_counter) {
        for (i_Z=0; i_Z< compound->nElements+1; i_Z++) {
          reuse_cum_xs_fluo[i_Z]     = cum_xs_fluo[i_Z];
          reuse_cum_xs_Compton[i_Z]  = cum_xs_Compton[i_Z];
          reuse_cum_xs_Rayleigh[i_Z] = cum_xs_Rayleigh[i_Z];
          reuse_xs_fluo[i_Z]         = xs_fluo[i_Z];
          reuse_xs_Compton[i_Z]      = xs_Compton[i_Z];
          reuse_xs_Rayleigh[i_Z]     = xs_Rayleigh[i_Z];
        }
        reuse_sigma_barn = sigma_barn;
      }

    } else {
      // reuse cached values (SPLIT), event_counter=0
      for (i_Z=0; i_Z< compound->nElements+1; i_Z++) {
        cum_xs_fluo[i_Z]      = reuse_cum_xs_fluo[i_Z];
        cum_xs_Compton[i_Z]   = reuse_cum_xs_Compton[i_Z];
        cum_xs_Rayleigh[i_Z]  = reuse_cum_xs_Rayleigh[i_Z];
        xs[FLUORESCENCE]     += reuse_cum_xs_fluo[i_Z];
        xs[COMPTON]          += reuse_cum_xs_Compton[i_Z];
        xs[RAYLEIGH]         += reuse_cum_xs_Rayleigh[i_Z];
        xs_fluo[i_Z]          = reuse_xs_fluo[i_Z];
        xs_Compton[i_Z]       = reuse_xs_Compton[i_Z];
        xs_Rayleigh[i_Z]      = reuse_xs_Rayleigh[i_Z];
      }
      sigma_barn = reuse_sigma_barn;
    }
//    printf("DEBUG: pi=%g p=%g ; ki=%g k=%g ; xs=[%g %g %g] type=FOLLOW sigma_barn=%g\n",
//        pi, p, ki, k, xs[0], xs[1], xs[2], sigma_barn);
    
    /* probability to absorb/scatter */
    my_s   = rho*100*sigma_barn; /* mu, 100: convert from barns to fm^2. my_s in [1/m] */
    d_path = ( dl0 +dl2 );  /* total path lenght in sample */
    
    /* Proba of transmission/interaction along length d_path */    
    p_trans = exp(-my_s*d_path); /* probability to not-interact (transmit) */
    p_scatt = 1 - p_trans;       /* portion of beam which scatters */
    
    /* force a given fraction of the beam to scatter */
    if (!event_counter && p_interact>0 && p_interact<=1) {
      /* we force a portion of the beam to interact */
      /* This is used to improve statistics */
      mc_trans = 1-p_interact;
    } else {
      mc_trans = p_trans; /* 1 - p_scatt */
    }
    mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */
    if (mc_scatt <= 0) ABSORB;
    
    if (!force_transmit && mc_scatt > 0 && (mc_scatt >= 1 || rand01() < mc_scatt)) {
      /* we "scatter" with one of the interaction processes */
        
      dl = -log(1 - rand0max((1 - exp(-my_s*d_path)))) / my_s; /* length */

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

      /* photon propagation to the scattering point */
      PROP_DL(dl);
      p *= fabs(p_scatt/mc_scatt); /* account for p_interact, lower than 1 */
      
    } else { // force_transmit
      /* we go through the material without interaction, and exit */
      if (type <0) type = TRANSMISSION; // transmission
      intersect = 0;
      PROP_DL(d_path);
      /* attenuate beam by portion which is scattered (and left along) */
      p *= p_trans;
      if (!event_counter && p_interact>0 && p_interact<=1) p /= mc_trans;
//      printf("DEBUG: pi=%g p=%g ; ki=%g k=%g ; xs=[%g %g %g] type=TRANSMIT\n",
//        pi, p, ki, k, xs[0], xs[1], xs[2]);
      break; // end while (intersect)
    }
    
  } /* if intersect (propagate) */

  if (intersect) { /* scattering event */
    int    i_Z, Z;
    double solid_angle=0;
    double theta, dsigma;
    double Ef, dE;
    double kf=k, kf_x, kf_y, kf_z;  // next particle direction

    /* MC choose process from cross sections 'xs': fluo, Rayleigh, Compton  */
    type = XRMC_SelectInteraction(xs, COMPTON+1, rand01()); // up to COMPTON
    
    /* choose Z (element) on associated XS, taking into account mass-fractions */
    if (flag_low_z) {
      i_Z = (int)floor(compound->nElements * rand01());
      switch (type) {
        case FLUORESCENCE:
          p *= xs_fluo[i_Z+1]/xs[type];
          break;
        case RAYLEIGH:
          p *= xs_Rayleigh[i_Z+1]/xs[type];
          break;
        case COMPTON:
          p *= xs_Compton[i_Z+1]/xs[type];
          break;
        default:
          // printf("%s: WARNING: process %i unknown. Absorb.\n", NAME_CURRENT_COMP, type);
          ABSORB;
      }
    } else {
      switch (type) {
        case FLUORESCENCE:
          i_Z = XRMC_SelectFromDistribution(cum_xs_fluo,     compound->nElements+1, rand01());
          break;
        case RAYLEIGH:
          if (!flag_rayleigh) ABSORB;
          i_Z = XRMC_SelectFromDistribution(cum_xs_Rayleigh, compound->nElements+1, rand01());
          break;
        case COMPTON:
          if (!flag_compton) ABSORB;
          i_Z = XRMC_SelectFromDistribution(cum_xs_Compton,  compound->nElements+1, rand01());
          break;
        default:
          // printf("%s: WARNING: process %i unknown. Absorb.\n", NAME_CURRENT_COMP, type);
          ABSORB;
      }
    }
    Z   = compound->Elements[i_Z];
    
    /* select outgoing vector (kf_x, kf_y, kf_z) and |kf| */
    if ((target_x || target_y || target_z)) {
      aim_x = target_x-x;       /* Vector pointing at target (anal./det.) */
      aim_y = target_y-y;
      aim_z = target_z-z;
    }
    if(focus_aw && focus_ah) {
      randvec_target_rect_angular(&kf_x, &kf_y, &kf_z, &solid_angle,
        aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
    } else if(focus_xw && focus_yh) {
      randvec_target_rect(&kf_x, &kf_y, &kf_z, &solid_angle,
        aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP);
    } else {
      randvec_target_circle(&kf_x, &kf_y, &kf_z, &solid_angle, aim_x, aim_y, aim_z, focus_r);
    }
    p *= solid_angle/(4*PI); // correct for selected solid-angle

    NORM(kf_x, kf_y, kf_z);  // normalize the outout direction |kf|=1
    
    // determine final energy (kf)
    switch (type) {
      case FLUORESCENCE: /* 0 Fluo: choose line */
        Ef      = XRMC_SelectFluorescenceEnergy(Z, E, &dE, flag_kissel, rand01()); // dE full-width in keV
        if (dE) {
          dE /= 2;  // half-width
          if (flag_lorentzian) dE  *= tan(PI/2*randpm1()); // Lorentzian distribution
          else                 dE  *= randnorm();          // Gaussian distribution
          Ef = Ef + dE;
        }
        kf      = Ef*E2K;
        theta   = acos(scalar_prod(kf_x,kf_y, kf_z, kx,ky,kz)/k);
        /* add polarisation factor */
        if (Ex!=0 || Ey!=0 || Ez!=0){
          double EE=sqrt(Ex*Ex+Ey*Ey+Ez*Ez);
          double s=scalar_prod(kf_x,kf_y,kf_z,Ex,Ey,Ez)/EE;
          p *= (1-s)*(1-s);
        } else {
          /*unpolarized light in - means an effective reduction according to only theta*/
          p *= (1+cos(theta)*cos(theta))*0.5;
        }
        n_fluo++;
        p_fluo += p;
        break;
        
      case RAYLEIGH:     /* 1 Rayleigh: Coherent, elastic */
        theta      = acos(scalar_prod(kf_x,kf_y, kf_z, kx,ky,kz)/k);
        dsigma     = DCSb_Rayl(Z,  E, theta, NULL); // [barn/at/st]
        p         *= 4*PI*dsigma/xs[RAYLEIGH];
        n_Rayleigh++;
        p_Rayleigh += p;
        break;
      
      case COMPTON:      /* 2 Compton: Incoherent: choose final energy */
        theta      = acos(scalar_prod(kf_x,kf_y, kf_z, kx,ky,kz)/k);
        dsigma     = DCS_Compt(Z,  E, theta, NULL); // [barn/at/st]
        kf         = ComptonEnergy(E, theta, NULL)*E2K; /* XRayLib, fraction of emc2=511 keV */
        p         *= 4*PI*dsigma/xs[COMPTON];
        n_Compton++;
        p_Compton += p;
        break;
    }
    
    kx = kf*kf_x;
    ky = kf*kf_y;
    kz = kf*kf_z;
    k  = kf;
    E  = k*K2E;
    SCATTER;
    event_counter++;
    
    /* exit if multiple scattering order has been reached */
    if (!order) break; // skip final absorption
    // stop when order has been reached, or weighting is very low
    if (order && (event_counter >= order || p/pi < 1e-15)) { force_transmit=1; }
 
  } // if intersect (scatter)
  
//  if (type>=0)
//  printf("DEBUG: pi=%g p=%g ; ki=%g k=%g ; xs=[%g %g %g] type=%i sigma_barn=%g %s\n",
//    pi, p, ki, k, xs[0], xs[1], xs[2], type, sigma_barn, reuse ? "SPLIT" : "");

} while(intersect); /* end do (intersect) (multiple scattering loop) */



%}

FINALLY %{
  FreeCompoundData(compound);

  if (n_fluo || n_Compton || n_Rayleigh)
    MPI_MASTER(
      printf("INFO: %s: scattered intensity: fluo=%g Compton=%g Rayleigh=%g\n",
             NAME_CURRENT_COMP, p_fluo, p_Compton, p_Rayleigh);
    );
%}

DISPLAY %{
    if (shape == 0) { /* cyl */
      circle("xz", 0,  yheight/2.0, 0, radius);
      circle("xz", 0, -yheight/2.0, 0, radius);
      line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
      line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
      line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
      line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
      if (thickness) {
        double radius_i=radius-thickness;
        circle("xz", 0,  yheight/2.0, 0, radius_i);
        circle("xz", 0, -yheight/2.0, 0, radius_i);
        line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0);
        line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0);
        line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i);
        line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i);
      }
    } else if (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);
      if (zdepth) {
        xmin = -0.5*xwidth;
        xmax =  0.5*xwidth;
        ymin = -0.5*yheight;
        ymax =  0.5*yheight;
        zmin = -0.5*zdepth;
        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);
      }
    } if (shape == 2) { /* sphere */
      circle("xy",0,0,0,radius);
      circle("xz",0,0,0,radius);
      circle("yz",0,0,0,radius);
    } else if (shape == 3) {        /* OFF file */
      off_display(offdata);
    }
%}

END