/*****************************************************************************
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: PowderN
*
* %Identification
* Written by: P. Willendrup, L. Chapon, K. Lefmann, A.B.Abrahamsen, N.B.Christensen, E.M.Lauridsen.
* Date: 4.2.98
* Origin: McXtrace release 1.2
* Modified by: KL, KN 20.03.98   (rewrite)
* Modified by: KL,    28.09.01   (two lines)
* Modified by: KL,    22.05.03   (background)
* Modified by: KL, PW 01.05.05   (N lines)
* Modified by: PW, LC 04.10.05   (Merge with Chapon Powder_multi)
* Modified by: PW, KL 05.06.07   (Concentricity)
* Modified by: EF,    17.10.08   (added any shape sample geometry)
* Modified by: EK,    28.10.10   (for X-ray use)
* Modified by: EK,               (Added options for incoherent scattering)
* Modified by: EF, June 2023     (can now read CIF files)
*
* General powder sample (N lines, single scattering, incoherent scattering)
*
* %Description
* General powder sample with
*         many scattering vectors
*         possibility for intrinsic line broadening
*         incoherent background ratio is computed from material datafile.
*         No multiple scattering. No secondary extinction.
*
* Based on Powder1/Powder2/Single_crystal.
* Geometry is a powder filled cylinder, sphere, box or any shape from an OFF file.
* Incoherent scattering is only provided here to account for a background.
* The efficient is highly improved when restricting the vertical scattering
*   range on the Debye-Scherrer cone (with 'd_phi' and 'focus_flip').
* The unit cell volume Vc may also be computed when giving the density,
*   the atomic/molecular weight and the number of atoms per unit cell.
*
* <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 xray with the object
*   transparently, so that it can be used like a regular sample object.
*   It supports the PLY, OFF and NOFF file format but not COFF (colored faces).
*   Such files may be generated from XYZ data using:
*   qhull < coordinates.xyz Qx Qv Tv o > geomview.off
*   or
*     powercrust coordinates.xyz
*   and viewed with geomview or java -jar jroff.jar (see below).
*   The default size of the object depends of the OFF file data, but its
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* If you use this component and produce valuable scientific results, please
* cite authors with references bellow (in <a href="#links">Links</a>).
*
* Example: PowderN(reflections = "c60.lau", d_phi = 15 , radius = 0.01,
*   yheight = 0.05, Vc = 1076.89, delta_d_d=0, DW=1)
*
* <b>Powder definition file format</b>
* Powder structure is specified with an ascii data file 'reflections'.
* The powder data are free-text column based files.
* The reflection list should be ordered by decreasing d-spacing values.
*       ... d ... F2
* Lines begining by '#' are read as comments (ignored) but they may contain
* the following keywords (in the header):
*   #Vc           <value of unit cell volume Vc [Angs^3]>
*   #Debye_Waller <value of Debye-Waller factor DW>
*   #delta_d_d/d    <value of delta_d_d/d width for all lines>
* These values are not read if entered as component parameters (Vc=...)
*
* The signification of the columns in the numerical block may be
* set using the 'format' parameter, by defining signification of the
* columns as a vector of indexes in the order
*   format={j,d,F2,DW,delta_d_d/d,1/2d,q,F}
* Signification of the symbols is given below. Indices start at 1.
* Indices with zero means that the column are not present, so that:
*   Crystallographica={ 4,5,7,0,0,0,0,0 }
*   Fullprof         ={ 4,0,8,0,0,5,0,0 }
*   Lazy             ={17,6,0,0,0,0,0,13}
*
* At last, the format may be overridden by direct definition of the
* column indexes in the file itself by using the following keywords
* in the header (e.g. '#column_j 4'):
*   #column_j     <index of the multiplicity 'j' column>
*   #column_d     <index of the d-spacing 'd' column [Angs]>
*   #column_F2    <index of the squared str. factor '|F|^2' column [b]>
*   #column_F     <index of the structure factor norm '|F|' column>
*   #column_DW    <index of the Debye-Waller factor 'DW' column>
*   #column_Dd    <index of the relative line width delta_d_d/d broadening 'Dd' column>
*   #column_inv2d <index of the 1/2d=sin(theta)/lambda 'inv2d' column>
*   #column_q     <index of the scattering wavevector 'q' column [Angs-1]>
*
* Last, CIF, FullProf and ShelX files can be read, and converted to F2(hkl) lists
* if 'cif2hkl' is installed. The CIF2HKL env variable can be used to point to a
* proper executable, else the McCode, then the system installed versions are used.
*
* <b>Concentricity</b>
*
* PowderN assumes 'concentric' shape, i.e. can contain other components inside its
* optional inner hollow. Example, Sample in Al cryostat:
*
*
* COMPONENT Cryo = PowderN(reflections="Al.laz", radius = 0.01, thickness = 0.001,
*                          concentric = 1, p_interact=0.1)
* AT (0,0,0) RELATIVE Somewhere
*
* COMPONENT Sample = some_other_component(with geometry FULLY enclosed in the hollow)
* AT (0,0,0) RELATIVE Somewhere
*
* COMPONENT Cryo2 = COPY(Cryo)(concentric = 0)
* AT (0,0,0) RELATIVE Somewhere
*
*
* (The second instance of the cryostat component can also be written out completely
*  using PowderN(...). In both cases, this second instance needs concentric = 0.)
* The concentric arrangment can not be used with OFF geometry specification.
*
* This sample component can advantageously benefit from the SPLIT feature, e.g.
* SPLIT COMPONENT pow = PowderN(...)
*
* %Parameters
* INPUT PARAMETERS
* radius:        [m] Outer radius of sample in (x,z) plane.
* xwidth:        [m] Horiz. dimension of sample, as a width.
* yheight:       [m] Height of sample y direction.
* zdepth:        [m] Depth of box sample.
* thickness:     [m] Thickness of hollow sample. Negative value extends the hollow volume outside of the box/cylinder.
* reflections: [str] Input file for reflections (LAZ LAU CIF, FullProf, ShelX). Use only incoherent scattering if NULL or "".
*
* Optional parameters:
* barns:         [1] Flag to indicate if |F|^2 from 'reflections' is in barns or fm^2, (barns=1 for laz/cif, barns=0 for lau type files).
* 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).
* d_omega:     [deg] Horizontal focus range (only for incoherent scattering), 0 for no focusing.
* d_phi:       [deg] Angle corresponding to the vertical angular range to focus to, e.g. detector height. 0 for no focusing.
* delta_d_d:     [1] Global relative Delta_d/d spreading when the 'w' column is not available. Use 0 if ideal.
* density:  [g/cm^3] Density of material. rho=density/weight/1e24*N_A.
* DW:            [1] Global Debye-Waller factor when the 'DW' column is not available. Use 1 if included in F2.
* focus_flip:    [1] Controls the sense of d_phi. If 0 d_phi is measured against the xz-plane. If !=0 d_phi is measured against zy-plane.
* format:       [{}] Name of the format, or list of column indexes (see Description). N.b. no quotes!
* 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.
* mat_format:   [{}] Format of the asorption parameter file.
* material: [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.
* nb_atoms:      [1] Number of sub-unit per unit cell, that is ratio of sigma for chemical formula to sigma per unit cell.
* p_inc:         [1] Fraction of incoherently scattered rays.
* p_interact:    [1] Fraction of events interacting with sample, e.g. 1-p_transmit-p_inc.
* p_transmit:    [1] Fraction of transmitted (only attenuated) rays.
* pack:          [1] Packing factor
* target_index:  [1] Relative index of component to focus incoherent scattering at, e.g. next is +1
* tth_sign:      [1] Sign of the scattering angle. If 0, the sign is chosen randomly (left and right). ONLY functional in combination with d_phi and ONLY applies to bragg lines.
* Vc:         [AA^3] Volume of unit cell=nb atoms per cell/density of atoms.
* weight:    [g/mol] Atomic/molecular weight of material.
*
* CALCULATED PARAMETERS:
* line_info: [struct]     Internal structure containing many members/info
* line_info.type: [char]  Interaction type of event 't'=Transmit, 'i'=Incoherent, 'c'=Coherent
* line_info.dq: [Angs^-1] Wavevector transfer of last coherent scattering event.
* 
* %Link
* See also: Single_crystal
* %Link
* See <a href="http://icsd.ill.fr">ICSD</a> Inorganic Crystal Structure Database
* <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 (free license)
* %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
* cif2hkl https://gitlab.com/soleil-data-treatment/soleil-software-projects/cif2hkl
* %Link
* material datafile obtained from http://physics.nist.gov/cgi-bin/ffast/ffast.pl
*
* %End
*****************************************************************************/
DEFINE COMPONENT PowderN

SETTING PARAMETERS (
  string reflections="NULL", string material="NULL", string geometry="NULL",
  vector format={0,0,0,0,0,0,0,0}, vector mat_format={0,0,0,0,0},
  radius=0, yheight=0, xwidth=0, zdepth=0, thickness=0,
  pack=1, Vc=0, delta_d_d=0, p_inc=0.1, p_transmit=0.1,
  DW=0, int nb_atoms=1, d_omega=0, d_phi=0, int tth_sign=0, p_interact=0,
  int concentric=0, density=0, weight=0, int barns=1, int focus_flip=0, int target_index=0)

/* X-ray 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 POWDERN_DECL
  #define POWDERN_DECL
  // single reflection data
  struct line_data {
    double F2;       /* Value of structure factor */
    double q;        /* Qvector */
    int j;           /* Multiplicity */
    double DWfactor; /* Debye-Waller factor */
    double w;        /* Intrinsic line width */
  };

  // component data
  struct line_info_struct {
    struct line_data* list; /* Reflection array */
    int count;              /* Number of reflections */
    double Dd;
    double DWfactor;
    double V_0;
    double rho;
    double at_weight;
    double at_nb;
    char compname[256];
    int flag_barns;
    int shape;           /* 0 cylinder, 1 box, 2 sphere, 3 OFF file */
    int column_order[8]; /* column signification */
    int flag_warning;
    double dq;      /* wavevector transfer [Angs-1] */
    double Epsilon; /* global strain in ppm */
    double XsectionFactor;
    double my_s_k2_sum;
    double my_a;
    double my_inc;
    double lfree; // store mean free path for the last event;
    double *w, *q, *my_s_k2;
    double radius_i, xwidth_i, yheight_i, zdepth_i;
    double k; /* last velocity (cached) */
    double Nq;
    int nb_reuses, nb_refl, nb_refl_count;
    double k_min, k_max;
    double xs_Nq[CHAR_BUF_LENGTH];
    double xs_sum[CHAR_BUF_LENGTH];
    unsigned int photon_passed;
    long xs_compute, xs_reuse, xs_calls;
    t_Table mat_table;
    int mat_column_order[5]; /*column signification for the coeff. in material data file*/
  };

  // shared functions **********************************************************

  // sign = PN_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
  PN_list_compare (const void* a, const void* b) {
    const struct line_data* pa = a;
    const struct line_data* pb = b;

    /* Sort by q */
    if (pa->q < pb->q)
      return -1;
    if (pa->q > pb->q)
      return 1;

    /* In case of tie, sort by F2 also */
    if (pa->F2 < pb->F2)
      return -1;
    if (pa->F2 > pb->F2)
      return 1;

    /* In case of tie, sort by j also */
    if (pa->j < pb->j)
      return -1;
    if (pa->j > pb->j)
      return 1;

    return 0;
  } /* PN_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 (char* infile, char* options) {
    char cmd[1024];
    int ret = 0;
    int found = 0;
    char* OUTFILE;
    char* inpath;

    // 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;
    inpath = malloc (1024);
    if (!inpath)
      return infile;

    // get input file path from read-table:Open_File
    FILE* f_infile = Open_File (infile, "r", inpath);
    if (!f_infile)
      return infile;
    fclose (f_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, inpath);
      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, inpath);
      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, inpath);
      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_line_data(filename, &info)
  //   read the given file to store F2(HKL) reflections and metadata.
  int
  read_line_data (char* SC_file, struct line_info_struct* info) {
    struct line_data* list = NULL;
    int size = 0;
    t_Table sTable; /* sample data table structure from SC_file */
    int i = 0;
    int mult_count = 0;
    char flag = 0;
    double q_count = 0, j_count = 0, F2_count = 0;
    char** parsing;
    int list_count = 0;
    double sum_F2 = 0;
    char* filename = NULL;

    if (!SC_file || !strlen (SC_file) || !strcmp (SC_file, "NULL")) {
      MPI_MASTER (printf ("PowderN: %s: Using incoherent elastic scattering only.\n", info->compname););
      info->count = 0;
      return (0);
    }
    filename = cif2hkl (SC_file, "--mode XRA");
    if (filename && filename != SC_file)
      info->flag_barns = 1; // cif2hkl returns barns

    /* read 1st block data from SC_file into sTable*/
    if (!filename || Table_Read (&sTable, filename, 1) < 0)
      exit (fprintf (stderr, "PowderN: ERROR: Could not open file %s - exiting!\n", filename));

    // remove temporary F2(hkl) file when giving CFL/CIF/ShelX file
    if (filename && filename != SC_file)
      unlink (filename);

    /* parsing of header */
    parsing = Table_ParseHeader (sTable.header, "Vc", "V_0", "column_j", "column_d", "column_F2", "column_DW", "column_Dd", "column_inv2d", "column_1/2d",
                                 "column_sintheta/lambda", "column_q", /* 10 */
                                 "DW", "Debye_Waller", "delta_d/d", "column_F ", "V_rho", "density", "weight", "nb_atoms", "multiplicity", NULL);

    if (parsing) {
      if (parsing[0] && !info->V_0)
        info->V_0 = atof (parsing[0]);
      if (parsing[1] && !info->V_0)
        info->V_0 = atof (parsing[1]);
      if (parsing[2])
        info->column_order[0] = atoi (parsing[2]);
      if (parsing[3])
        info->column_order[1] = atoi (parsing[3]);
      if (parsing[4])
        info->column_order[2] = atoi (parsing[4]);
      if (parsing[5])
        info->column_order[3] = atoi (parsing[5]);
      if (parsing[6])
        info->column_order[4] = atoi (parsing[6]);
      if (parsing[7])
        info->column_order[5] = atoi (parsing[7]);
      if (parsing[8])
        info->column_order[5] = atoi (parsing[8]);
      if (parsing[9])
        info->column_order[5] = atoi (parsing[9]);
      if (parsing[10])
        info->column_order[6] = atoi (parsing[10]);
      if (parsing[11] && info->DWfactor <= 0)
        info->DWfactor = atof (parsing[11]);
      if (parsing[12] && info->DWfactor <= 0)
        info->DWfactor = atof (parsing[12]);
      if (parsing[13] && info->Dd < 0)
        info->Dd = atof (parsing[13]);
      if (parsing[14])
        info->column_order[7] = atoi (parsing[14]);
      if (parsing[15] && !info->V_0)
        info->V_0 = 1 / atof (parsing[15]);
      if (parsing[16] && !info->rho)
        info->rho = atof (parsing[16]);
      if (parsing[17] && !info->at_weight)
        info->at_weight = atof (parsing[17]);
      if (parsing[18] && info->at_nb <= 1)
        info->at_nb = atof (parsing[18]);
      if (parsing[19] && info->at_nb <= 1)
        info->at_nb = atof (parsing[19]);
      for (i = 0; i <= 19; i++)
        if (parsing[i])
          free (parsing[i]);
      free (parsing);
    }

    if (!sTable.rows)
      exit (fprintf (stderr,
                     "PowderN: %s: Error: The number of rows in %s "
                     "should be at least %d\n",
                     info->compname, SC_file, 1));
    else
      size = sTable.rows;

    MPI_MASTER (Table_Info (sTable););

    if (filename == SC_file) { // only when not from cif2hkl
      if (info->column_order[0] == 4 && info->flag_barns != 0)
        MPI_MASTER (printf ("PowderN: %s: Powder file probably of type Crystallographica/Fullprof (lau)\n"
                            "WARNING: but F2 unit is set to barns=1 (barns). Intensity might be 100 times too high.\n",
                            info->compname););
      if (info->column_order[0] == 17 && info->flag_barns == 0)
        MPI_MASTER (printf ("PowderN: %s: Powder file probably of type Lazy Pulver (laz)\n"
                            "WARNING: but F2 unit is set to barns=0 (fm^2). Intensity might be 100 times too low.\n",
                            info->compname););
    }
    /* allocate line_data array */
    list = (struct line_data*)malloc (size * sizeof (struct line_data));

    for (i = 0; i < size; i++) {
      /*      printf("Reading in line %i\n",i);*/
      double j = 0, d = 0, w = 0, q = 0, DWfactor = 0, F2 = 0;
      int index;

      if (info->Dd >= 0)
        w = info->Dd;
      if (info->DWfactor > 0)
        DWfactor = info->DWfactor;

      /* get data from table using columns {j d F2 DW Dd inv2d q F} */
      /* column indexes start at 1, thus need to substract 1 */
      if (info->column_order[0] > 0)
        j = Table_Index (sTable, i, info->column_order[0] - 1);
      if (info->column_order[1] > 0)
        d = Table_Index (sTable, i, info->column_order[1] - 1);
      if (info->column_order[2] > 0)
        F2 = Table_Index (sTable, i, info->column_order[2] - 1);
      if (info->column_order[3] > 0)
        DWfactor = Table_Index (sTable, i, info->column_order[3] - 1);
      if (info->column_order[4] > 0)
        w = Table_Index (sTable, i, info->column_order[4] - 1);
      if (info->column_order[5] > 0 && !(info->column_order[1] > 0)) // Only use if d not read already
      {
        d = Table_Index (sTable, i, info->column_order[5] - 1);
        d = (d > 0 ? 1 / d / 2 : 0);
      }
      if (info->column_order[6] > 0 && !(info->column_order[1] > 0)) // Only use if d not read already
      {
        q = Table_Index (sTable, i, info->column_order[6] - 1);
        d = (q > 0 ? 2 * PI / q : 0);
      }
      if (info->column_order[7] > 0 && !F2) {
        F2 = Table_Index (sTable, i, info->column_order[7] - 1);
        F2 *= F2;
      }

      /* assign and check values */
      j = (j > 0 ? j : 0);
      q = (d > 0 ? 2 * PI / d : 0); /* this is q */
      DWfactor = (DWfactor > 0 ? DWfactor : 1);
      w = (w > 0 ? w : 0); /* this is q and d relative spreading */
      F2 = (F2 >= 0 ? F2 : 0);
      if (j == 0 || q == 0) {
        MPI_MASTER (printf ("PowderN: %s: line %i has invalid definition (mult=0 or q=0 or d=0)\n", info->compname, i););
        continue;
      }
      list[list_count].j = j;
      list[list_count].q = q;
      list[list_count].DWfactor = DWfactor;
      list[list_count].w = w;
      list[list_count].F2 = F2;
      sum_F2 += F2;

      /* adjust multiplicity if j-column + multiple d-spacing lines */
      /* if  d = previous d, increase line duplication index */
      if (!q_count)
        q_count = q;
      if (!j_count)
        j_count = j;
      if (!F2_count)
        F2_count = F2;
      if (fabs (q_count - q) < 0.0001 * fabs (q) && fabs (F2_count - F2) < 0.0001 * fabs (F2) && j_count == j) {
        mult_count++;
        flag = 0;
      } else
        flag = 1;
      if (i == size - 1)
        flag = 1;
      /* else if d != previous d : just passed equivalent lines */
      if (flag) {
        if (i == size - 1)
          list_count++;
        /*   if duplication index == previous multiplicity */
        /*      set back multiplicity of previous lines to 1 */
        if ((mult_count && list_count > 0)
            && (mult_count == list[list_count - 1].j || ((list_count < size) && (i == size - 1) && (mult_count == list[list_count].j)))) {
          MPI_MASTER (printf ("PowderN: %s: Set multiplicity to 1 for lines [%i:%i]\n"
                              "         (d-spacing %g is duplicated %i times)\n",
                              info->compname, list_count - mult_count, list_count - 1, list[list_count - 1].q, mult_count););
          for (index = list_count - mult_count; index < list_count; list[index++].j = 1)
            ;
          mult_count = 1;
          q_count = q;
          j_count = j;
          F2_count = F2;
        }
        if (i == size - 1)
          list_count--;
        flag = 0;
      }
      list_count++;
    } /* end for */

    Table_Free (&sTable);

    if (!sum_F2) {
      MPI_MASTER (fprintf (stderr, "PowderN: %s: ERROR: all %i structure factors in file '%s' are null. Check the reflection list.\n", info->compname, list_count,
                           SC_file););
      return (0);
    }

    /* sort the list with increasing q */
    qsort (list, list_count, sizeof (struct line_data), PN_list_compare);

    MPI_MASTER (printf ("PowderN: %s: Read %i reflections from file '%s'\n", info->compname, list_count, SC_file););

    info->list = list;
    info->count = list_count;

    return (list_count);
  } /* read_line_data */

  /* calc_xsect: computes the number of possible reflections (return value), and the total xsection 'sum' */
  /* this routine looks for a pre-computed value in the Nq and sum cache tables               */
  /* when found, the earch starts from the corresponding lower element in the table           */
  #pragma acc routine
  int
  calc_xsect (double k, double* q, double* my_s_k2, int count, double* sum, struct line_info_struct* line_info) {
    int Nq = 0, line = 0, line0 = 0;
    *sum = 0;

    /* check if a line_info element has been recorded already - not on OpenACC */
    #ifndef OPENACC
    if (k >= line_info->k_min && k <= line_info->k_max && line_info->photon_passed >= CHAR_BUF_LENGTH) {
      line = (int)floor (k - line_info->k_min) * CHAR_BUF_LENGTH / (line_info->k_max - line_info->k_min);
      Nq = line_info->xs_Nq[line];
      *sum = line_info->xs_sum[line];
      if (!Nq && *sum == 0) {
        /* not yet set: we compute the sum up to the corresponding wavevector in the table cache */
        double line_k = line_info->k_min + line * (line_info->k_max - line_info->k_min) / CHAR_BUF_LENGTH;
        for (line0 = 0; line0 < count; line0++) {
          if (q[line0] <= 2 * line_k) { /* q < 2*kf: restrict structural range */
            *sum += my_s_k2[line0];
            if (Nq < line0 + 1)
              Nq = line0 + 1; /* determine maximum line index which can scatter */
          } else
            break;
        }
        line_info->xs_Nq[line] = Nq;
        line_info->xs_sum[line] = *sum;
        line_info->xs_compute++;
      } else
        line_info->xs_reuse++;
      line0 = Nq;
    }

    line_info->xs_calls++;
    #endif

    for (line = line0; line < count; line++) {
      if (q[line] <= 2 * k) { /* q < 2*kf: restrict structural range */
        *sum += my_s_k2[line];
        if (Nq < line + 1)
          Nq = line + 1; /* determine maximum line index which can scatter */
      } else
        break;
    }

    return (Nq);
  } /* calc_xsect */

  #pragma acc routine
  int
  calc_abs_xsect (double k, double* abs, struct line_info_struct* line_info) {
    /* compute the absorption cross section.
     * amounts to a table lookup in material data_file by energy*/
    double e = k * K2E;
    *abs = Table_Value ((line_info->mat_table), e, line_info->mat_column_order[1]);
    *abs *= 100 * line_info->rho;
    return 0;
  }

  #pragma acc routine
  int
  calc_inc_xsect (double k, double* inc, struct line_info_struct* line_info) {
    /*TODO choose one of these bits*/

    /* compute the incoherent cross section.
     * amounts to a table lookup in material data_file by energy
     * In time this should be improved to actually handle Compton scattering.*/
    double e = k * K2E;

    double correction;
    int i = 2;
    /* To find the inc. scattering column in material datafiles,
     * Find the first non-zero column index (btw. 2..4).*/
    for (i = 2; i < 5; i++) {
      if (line_info->mat_column_order[i])
        break;
    }

    switch (i) {
    case 2:
      /*inc is explicitly in column i - no correction term.*/
      correction = 0;
      break;
    case 3:
      /*inc is given as inc+ćoh in column i*/
      correction = line_info->my_s_k2_sum / (k * k);
      break;
    case 4:
      /*inc is given as total attenuation = inc+coh+abs in column i*/
      correction = line_info->my_s_k2_sum / (k * k) + line_info->my_a;
      break;
    case 5:
      /*means no non-zero index has been found - set xsect to 0*/
      *inc = 0;
      return 1;
    }

    *inc = Table_Value ((line_info->mat_table), e, line_info->mat_column_order[i] - 1) - correction;
    if (*inc < 0) {
      *inc = 0;
    }
    *inc *= 100 * line_info->rho; /*scale by rho and 100 to get in units m^-1*/
    return 0;
  }

  #endif /* !POWDERN_DECL */
%}

DECLARE
%{
  struct line_info_struct line_info;
  double* columns;
  double* mat_columns;
  off_struct offdata;
  double tgt_x;
  double tgt_y;
  double tgt_z;
%}

INITIALIZE
%{
  /* We ought to clean up the columns variable as format is now a proper vector/array */
  columns = format;
  mat_columns = mat_format;
  int i = 0;
  struct line_data* L;
  line_info.Dd = delta_d_d;
  line_info.DWfactor = DW;
  line_info.V_0 = Vc;
  line_info.rho = density;
  line_info.at_weight = weight;
  line_info.at_nb = nb_atoms;
  line_info.flag_barns = barns;
  line_info.shape = 0;
  line_info.flag_warning = 0;
  line_info.radius_i = line_info.xwidth_i = line_info.yheight_i = line_info.zdepth_i = 0;
  line_info.k = 0;
  line_info.Nq = 0;
  line_info.k_min = FLT_MAX;
  line_info.k_max = 0;
  line_info.photon_passed = 0;
  line_info.nb_reuses = line_info.nb_refl = line_info.nb_refl_count = 0;
  line_info.xs_compute = line_info.xs_reuse = line_info.xs_calls = 0;
  for (i = 0; i < 8; i++) {
    line_info.column_order[i] = (int)columns[i];
  }
  for (i = 0; i < 4; i++) {
    line_info.mat_column_order[i] = (int)mat_columns[i];
  }
  strncpy (line_info.compname, NAME_CURRENT_COMP, 256);

  line_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
    fprintf (stderr, "PowderN: %s: Error: You are attempting to use an OFF 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)) {
      line_info.shape = 3;
      thickness = 0;
      concentric = 0;
    }
    #endif
  } else if (xwidth && yheight && zdepth)
    line_info.shape = 1; /* box */
  else if (radius > 0 && yheight)
    line_info.shape = 0; /* cylinder */
  else if (radius > 0 && !yheight)
    line_info.shape = 2; /* sphere */

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

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

  if (thickness > 0) {
    if (radius > thickness) {
      line_info.radius_i = radius - thickness;
    } else {
      if (xwidth > 2 * thickness)
        line_info.xwidth_i = xwidth - 2 * thickness;
      if (yheight > 2 * thickness)
        line_info.yheight_i = yheight - 2 * thickness;
      if (zdepth > 2 * thickness)
        line_info.zdepth_i = zdepth - 2 * thickness;
    }
  } else if (thickness < 0) {
    thickness = fabs (thickness);
    if (radius) {
      line_info.radius_i = radius;
      radius = line_info.radius_i + thickness;
    } else {
      line_info.xwidth_i = xwidth;
      line_info.yheight_i = yheight;
      line_info.zdepth_i = zdepth;
      xwidth = xwidth + 2 * thickness;
      yheight = yheight + 2 * thickness;
      zdepth = zdepth + 2 * thickness;
    }
  }

  if (!line_info.yheight_i) {
    line_info.yheight_i = yheight;
  }

  if (!p_interact) {
    if (!p_transmit) {
      MPI_MASTER (fprintf (stderr, "PowderN: %s: WARNING: p_interact=0, adjusting to 0.01, to avoid algorithm instability\n", NAME_CURRENT_COMP););
      p_interact = 1e-2;
    } else {
      p_interact = 1 - p_transmit - p_inc;
    }
  }
  if (!p_transmit) {
    MPI_MASTER (fprintf (stderr, "PowderN: %s: WARNING: p_transmit=0, adjusting to 0.01, to avoid algorithm instability\n", NAME_CURRENT_COMP););
    p_transmit = 1e-2;
  }
  double p_sum = p_interact + p_inc + p_transmit;
  p_interact = p_interact / p_sum;
  p_inc = p_inc / p_sum;
  p_transmit = p_transmit / p_sum;

  if (concentric) {
    MPI_MASTER (printf ("PowderN: %s: Concentric mode - remember to include the 'opposite' copy of this component !\n"
                        "WARNING  The equivalent, 'opposite' comp should have concentric=0\n",
                        NAME_CURRENT_COMP););
    if (p_transmit < 0.1) {
      MPI_MASTER (printf ("PowderN: %s: Concentric mode and p_transmit<0.1 !\n"
                          "WARNING  Consider increasing p_transmit as few particles will reach the inner hollow.\n",
                          NAME_CURRENT_COMP););
    }
  }

  if (reflections && strlen (reflections) && strcmp (reflections, "NULL") && strcmp (reflections, "0")) {
    i = read_line_data (reflections, &line_info);
    if (i == 0)
      exit (fprintf (stderr,
                     "PowderN: %s: reflection file %s is not valid.\n"
                     "ERROR    Please check file format.\n",
                     NAME_CURRENT_COMP, reflections));
  }

  /* compute the scattering unit density from material weight and density */
  /* the weight of the scattering element is the chemical formula molecular weight
   * times the nb of chemical formulae in the scattering element (nb_atoms) */
  if (!line_info.V_0 && line_info.at_nb > 0 && line_info.at_weight > 0 && line_info.rho > 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 */
    line_info.V_0 = line_info.at_nb / (line_info.rho / line_info.at_weight / 1e24 * 6.02214199e23);
  }

  if (line_info.V_0 <= 0)
    MPI_MASTER (printf ("PowderN: %s: density/unit cell volume is NULL (Vc). Unactivating component.\n", NAME_CURRENT_COMP););

  if (line_info.V_0 > 0 && p_inc > 0 && (!(strlen (material) && strcmp (material, "NULL")))) {
    p_inc = 0; // no inc scatt (no material file for abs)
  }

  /* Factor 100 to convert from barns to fm^2 */
  line_info.XsectionFactor = line_info.flag_barns ? 100 : 1;

  if (line_info.V_0 > 0 && i) {
    L = line_info.list;

    line_info.q = malloc (line_info.count * sizeof (double));
    line_info.w = malloc (line_info.count * sizeof (double));
    line_info.my_s_k2 = malloc (line_info.count * sizeof (double));
    if (!line_info.q || !line_info.w || !line_info.my_s_k2)
      exit (fprintf (stderr, "PowderN: %s: ERROR allocating memory (init)\n", NAME_CURRENT_COMP));
    for (i = 0; i < line_info.count; i++) {
      line_info.my_s_k2[i] = 4 * PI * PI * PI * 1.0 * (L[i].DWfactor ? L[i].DWfactor : 1) / (line_info.V_0 * line_info.V_0) * (L[i].j * L[i].F2 / L[i].q)
                             * line_info.XsectionFactor;
      line_info.q[i] = L[i].q;
      line_info.w[i] = L[i].w;
    }
  }

  if (material && strlen (material) && strcmp (material, "NULL")) {
    int status;
    char** parsing;
    if ((status = Table_Read (&(line_info.mat_table), material, 0)) == -1) {
      fprintf (stderr, "PowderN: %s Error reading material data from file %s.\n", NAME_CURRENT_COMP, material);
    }
    parsing = Table_ParseHeader (line_info.mat_table.header, "column_e", "column_abs", "column_inc", "column_cohinc", "column_tot", NULL);

    if (parsing) {
      int i;
      for (i = 0; i < 5; i++) {
        if (parsing[i])
          line_info.mat_column_order[i] = atoi (parsing[i]);
      }
    }
  }
%}

TRACE
%{
  double l0, l1, l2, l3, k, k1, l_full, l, l_1, dl, alpha0, alpha, theta, my_s, my_s_n, sg;
  double solid_angle, photontype, ptype;
  double arg, tmp_kx, tmp_ky, tmp_kz, kout_x, kout_y, kout_z, nx, ny, nz, pmul = 1;
  int line;
  char intersect = 0;
  char intersecti = 0;

  // Variables calculated within thread for thread purpose only
  char type = '\0';
  int itype = 0;
  double d_phi_thread = d_phi;
  // These ones are injected back to struct at the end of TRACE in non-OpenACC case
  int nb_reuses = line_info.nb_reuses;
  int nb_refl = line_info.nb_refl;
  int nb_refl_count = line_info.nb_refl_count;
  double kcache = line_info.k;
  double Nq = line_info.Nq;
  double k_min = line_info.k_min;
  double k_max = line_info.k_max;
  double lfree = line_info.lfree;
  unsigned int photon_passed = line_info.photon_passed;
  long xs_compute = line_info.xs_compute;
  long xs_reuse = line_info.xs_reuse;
  long xs_calls = line_info.xs_calls;
  // Number of warnings in this specific TRACE run,
  // added to comp struct end of TRACE
  int flag_warning = 0;
  double dq = line_info.dq;

  double my_s_k2_sum = line_info.my_s_k2_sum;
  double my_a = line_info.my_a;
  double my_inc = line_info.my_inc;

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

  if (line_info.V_0 > 0 && (line_info.count || line_info.my_inc)) {
    if (line_info.shape == 1) {
      intersect = box_intersect (&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
      intersecti = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
    } else if (line_info.shape == 0) {
      intersect = cylinder_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
      intersecti = cylinder_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
    } else if (line_info.shape == 2) {
      intersect = sphere_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius);
      intersecti = sphere_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i);
    }
    #ifdef USE_OFF
    else if (line_info.shape == 3) {
      intersect = off_x_intersect (&l0, &l3, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata);
      intersecti = 0;
    }
    #endif
  }

  if (intersect && l3 > 0) {

    if (concentric) {
      /* Set up for concentric case */
      /* 'Remove' the backside of this comp */
      if (!intersecti) {
        l1 = (l3 + l0) / 2;
      }
      l2 = l1;
      l3 = l1;
      dl = -1.0 * rand01 (); /* In case of scattering we will scatter on 'forward' part of sample */
    } else {
      if (!intersecti) {
        l1 = (l3 + l0) / 2;
        l2 = l1;
      }
      dl = randpm1 (); /* Possibility to scatter at all points in line of sight */
    }

    /* X-ray enters at t=l0. */
    if (l0 < 0)
      l0 = 0; /* already in sample */
    if (l1 < 0)
      l1 = 0; /* already in inner hollow */
    if (l2 < 0)
      l2 = 0; /* already past inner hollow */
    k = sqrt (kx * kx + ky * ky + kz * kz);
    l_full = l3 - l2 + l1 - l0;

    if (line_info.photon_passed < CHAR_BUF_LENGTH) {
      if (k < line_info.k_min)
        line_info.k_min = k;
      if (k > line_info.k_max)
        line_info.k_max = k;
      photon_passed++;
    }

    /* Calculate total scattering cross section at relevant wavevector.
       Reuse calculation if possible (i.e if photon is identical ot preceeding one and not GPU.*/
    #ifndef OPENACC
    if (fabs (k - line_info.k) < 1e-15) {
      line_info.nb_reuses++;
      my_s_k2_sum = line_info.my_s_k2_sum;
      my_inc = line_info.my_inc;
      my_a = line_info.my_a;
    } else {
      #endif
      Nq = calc_xsect (k, line_info.q, line_info.my_s_k2, line_info.count, &my_s_k2_sum, &line_info);
      calc_abs_xsect (k, &my_a, &line_info);
      calc_inc_xsect (k, &my_inc, &line_info);
      if (pack != 1) {
        /*apply packing factor correction*/
        my_a *= pack;
        my_inc *= pack;
        my_s_k2_sum *= pack;
      }
      #ifndef OPENACC
      line_info.my_s_k2_sum = my_s_k2_sum;
      line_info.my_inc = my_inc;
      line_info.my_a = my_a;
      line_info.Nq = Nq;
      line_info.k = k;
      line_info.nb_refl += line_info.Nq;
      line_info.nb_refl_count++;
    }
    #endif

    if (l3 < 0) {
      l3 = 0; /* Already past sample?! */
      if (line_info.flag_warning < 10)
        printf ("PowderN: %s: WARNING: xray has already passed us? (Skipped).\n"
                "         In concentric geometry, this may be caused by a missing concentric=0 option in 2nd enclosing instance.\n",
                NAME_CURRENT_COMP);
      flag_warning++;
    } else {
      if (dl < 0) {                 /* Calculate scattering point position */
        dl = fabs (dl) * (l1 - l0); /* 'Forward' part */
      } else {
        dl = dl * (l3 - l2) + (l2 - l0); /* Possibly also 'backside' part */
      }

      my_s = my_s_k2_sum / (k * k) + my_inc;
      /* Total attenuation from scattering*/

      ptype = rand01 ();
      /* How to handle this one? Transmit (1) / Incoherent (2) / Coherent (3) ? */
      if (ptype < p_transmit) {
        photontype = 1;
        l = l_full; /* Passing through, full length */
        PROP_DL (l3);
      } else if (p_inc > 0 && ptype >= p_transmit && ptype < (p_transmit + p_inc)) {
        photontype = 2;
        l = dl;            /* Penetration in sample */
        PROP_DL (dl + l0); /* Point of scattering */
        SCATTER;
      } else if (ptype >= p_transmit + p_inc) {
        photontype = 3;
        l = dl;
        PROP_DL (dl + l0); /* Point of scattering */
        SCATTER;
      } else {
        exit (fprintf (stderr, "PowderN %s: DEAD - this shouldn't happen!\n", NAME_CURRENT_COMP));
      }

      if (photontype == 3) { /* Make coherent scattering event */
        if (line_info.count > 0) {
          /* choose line */
          if (Nq > 1)
            line = floor (Nq * rand01 ()); /* Select between Nq powder lines */
          else
            line = 0;
          if (line_info.w[line])
            arg = line_info.q[line] * (1 + line_info.w[line] * randnorm ()) / (2.0 * k);
          else
            arg = line_info.q[line] / (2.0 * k);
          my_s_n = pack * line_info.my_s_k2[line] / (k * k);
          if (fabs (arg) > 1)
            ABSORB; /* No bragg scattering possible*/
          if (tth_sign == 0) {
            sg = randpm1 ();
            if (sg > 0)
              sg = 1;
            else
              sg = -1;
          } else {
            sg = tth_sign / fabs (tth_sign);
          }
          theta = asin (arg); /* Bragg scattering law */

          /* Choose point on Debye-Scherrer cone */
          if (d_phi_thread) { /* relate height of detector to the height on DS cone */
            arg = sin (d_phi_thread * DEG2RAD / 2) / sin (2 * theta);
            /* If full Debye-Scherrer cone is within d_phi, don't focus */
            if (arg < -1 || arg > 1)
              d_phi_thread = 0;
            /* Otherwise, determine alpha to rotate from scattering plane
               into d_phi focusing area*/
            else
              alpha = 2 * asin (arg);
          }
          if (d_phi_thread) {
            /* Focusing */
            alpha = fabs (alpha);
            alpha0 = 0.5 * randpm1 () * alpha;
            if (focus_flip) {
              alpha0 += M_PI_2;
            }
          } else
            alpha0 = M_PI * randpm1 ();

          /* now find a nearly vertical rotation axis:
           * Either
           *  (k along Z) x (X axis) -> nearly Y axis
           * Or
           *  (k along X) x (Z axis) -> nearly Y axis
           */

          /* update originating from McStas, added by JS, 1/7/2017
             If a target is defined, try to define vertical axis as a normal to the plane
                defined by the incident particle velocity and target position.
                Check that k is not ~ parallel to the target direction.
          */
          double knorm = 0.0;
          if (target_index) {
            vec_prod (tmp_kx, tmp_ky, tmp_kz, kx, ky, kz, tgt_x, tgt_y, tgt_z);
            knorm = sqrt (tmp_kx * tmp_kx + tmp_ky * tmp_ky + tmp_kz * tmp_kz) / k;
          }
          // no target or direction is nearly parallel to k:
          if (knorm < 0.01) {
            if (fabs (kx / k) < fabs (kz / k)) {
              nx = 1;
              ny = 0;
              nz = 0;
            } else {
              nx = 0;
              ny = 0;
              nz = 1;
            }
            vec_prod (tmp_kx, tmp_ky, tmp_kz, kx, ky, kz, nx, ny, nz);
          }

          /* k_out = rotate 'k' by 2*theta around tmp_k: Bragg angle */
          rotate (kout_x, kout_y, kout_z, kx, ky, kz, 2 * sg * theta, tmp_kx, tmp_ky, tmp_kz);

          /* tmp_k = rotate k_out by alpha0 around 'k' (Debye-Scherrer cone) */
          rotate (tmp_kx, tmp_ky, tmp_kz, kout_x, kout_y, kout_z, alpha0, kx, ky, kz);
          kx = tmp_kx;
          ky = tmp_ky;
          kz = tmp_kz;

          /*weight the outgoing signal according to polarization*/
          if (Ex != 0 || Ey != 0 || Ez != 0) {
            double EE = sqrt (Ex * Ex + Ey * Ey + Ez * Ez);
            double s = scalar_prod (kx, ky, kz, Ex, Ey, Ez) / k / 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;
          }

          /* Since now scattered and new direction given, calculate path to exit */
          if (line_info.shape == 1) {
            intersect = box_intersect (&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
            intersecti = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
          } else if (line_info.shape == 0) {
            intersect = cylinder_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
            intersecti = cylinder_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
          } else if (line_info.shape == 2) {
            intersect = sphere_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius);
            intersecti = sphere_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i);
          } else if (line_info.shape == 3) {
            intersect = off_x_intersect (&l0, &l3, NULL, NULL, x, y, z, kx, ky, kz, offdata);
            intersecti = 0;
          }

          if (!intersect) {
            /* Strange error: did not hit cylinder */
            if (line_info.flag_warning < 10)
              printf ("PowderN: %s: WARNING: Did not hit sample from inside (coh). ABSORB.\n", NAME_CURRENT_COMP);
            flag_warning++;
            #pragma acc atomic write
            line_info.flag_warning += flag_warning;
            ABSORB;
          }

          if (!intersecti) {
            l1 = (l3 + l0) / 2;
            l2 = l1;
          }

          if (concentric && intersecti) {
            /* In case of concentricity, 'remove' backward wall of sample */
            l2 = l1;
            l3 = l1;
          }

          if (l0 < 0)
            l0 = 0; /* already in sample */
          if (l1 < 0)
            l1 = 0; /* already in inner hollow */
          if (l2 < 0)
            l2 = 0; /* already past inner hollow */

          l_1 = l3 - l2 + l1 - l0; /* Length to exit */

          pmul *= Nq * l_full * my_s_n * exp (-(my_a + my_s) * (l + l_1)) / (1 - (p_inc + p_transmit));
          /* Correction in case of d_phi focusing - BUT only when d_phi != 0 */
          if (d_phi_thread) {
            pmul *= alpha / PI;
            if (tth_sign)
              pmul *= 0.5;
          }

          type = 'c';
          itype = 1;
          dq = line_info.q[line];
        } /* else transmit <-- No powder lines in file */
      } /* Coherent scattering event */
      else if (p_inc > 0 && photontype == 2) { /* Make incoherent scattering event */
        /*should be replaced by Compton scattering and/or TDS*/
        if (d_omega && d_phi_thread) {
          randvec_target_rect_angular (&kx, &ky, &kz, &solid_angle, tgt_x, tgt_y, tgt_z, d_omega * DEG2RAD, d_phi_thread * DEG2RAD, ROT_A_CURRENT_COMP);
        } else if (d_phi_thread) {
          randvec_target_rect_angular (&kx, &ky, &kz, &solid_angle, tgt_x, tgt_y, tgt_z, 2 * PI, d_phi_thread * DEG2RAD, ROT_A_CURRENT_COMP);
        } else {
          randvec_target_circle (&kx, &ky, &kz, &solid_angle, 0, 0, 1, 0);
        }
        k1 = sqrt (kx * kx + ky * ky + kz * kz);
        kx *= k / k1;
        ky *= k / k1;
        kz *= k / k1;

        /* Since now scattered and new direction given, calculate path to exit */
        if (line_info.shape == 1) {
          intersect = box_intersect (&l0, &l3, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
          intersecti = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.xwidth_i, line_info.yheight_i, line_info.zdepth_i);
        } else if (line_info.shape == 0) {
          intersect = cylinder_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius, yheight);
          intersecti = cylinder_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i, line_info.yheight_i);
        } else if (line_info.shape == 2) {
          intersect = sphere_intersect (&l0, &l3, x, y, z, kx, ky, kz, radius);
          intersecti = sphere_intersect (&l1, &l2, x, y, z, kx, ky, kz, line_info.radius_i);
        }
        #ifdef USE_OFF
        else if (line_info.shape == 3) {
          intersect = off_x_intersect (&l0, &l3, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata);
          intersecti = 0;
        }
        #endif

        if (!intersect) {
          /* Strange error: did not hit cylinder */
          if (line_info.flag_warning < 10)
            printf ("PowderN: %s: WARNING: Did not hit sample from inside (inc). ABSORB.\n", NAME_CURRENT_COMP);
          flag_warning++;
          #pragma acc atomic write
          line_info.flag_warning += flag_warning;
          ABSORB;
        }

        if (!intersecti) {
          l1 = (l3 + l0) / 2;
          l2 = l1;
        }

        if (concentric && intersecti) {
          /* In case of concentricity, 'remove' backward wall of sample */
          l2 = l1;
          l3 = l1;
        }

        if (l0 < 0)
          l0 = 0; /* already in sample */
        if (l1 < 0)
          l1 = 0; /* already in inner hollow */
        if (l2 < 0)
          l2 = 0; /* already past inner hollow */

        l_1 = (l3 - l2 + l1 - l0); /* Length to exit */

        pmul *= l_full * my_inc * exp (-(my_a + my_s) * (l + l_1)) / (p_inc);
        pmul *= solid_angle / (4 * PI);

        type = 'i';
        itype = 2;

      } /* Incoherent scattering event */
      else if (photontype == 1) {
        /* Make transmitted (absorption-corrected) event */
        /* No coordinate changes here, simply change xray weight */
        pmul *= exp (-(my_a + my_s) * (l)) / (p_transmit);

        type = 't';
        itype = 3;
      }
      p *= pmul;
    } /* Photon leaving since it has passed already */
  } /* else transmit non interacting xrays */

  // Inject these back to component struct in non-OpenACC case
  #ifndef OPENACC
  line_info.nb_reuses = nb_reuses;
  line_info.nb_refl = nb_refl;
  line_info.nb_refl_count = nb_refl_count;
  line_info.k = kcache;
  line_info.Nq = Nq;
  line_info.k_min = k_min;
  line_info.k_max = k_max;
  line_info.lfree = lfree;
  line_info.xs_compute = xs_compute;
  line_info.xs_reuse = xs_reuse;
  line_info.xs_calls = xs_calls;
  line_info.dq = dq;
  line_info.photon_passed = photon_passed;
  #endif

  // Add warning status to comp struct (atomic on GPU)
  #pragma acc atomic write
  line_info.flag_warning += flag_warning;
%}

FINALLY
%{
  free (line_info.list);
  free (line_info.q);
  free (line_info.w);
  free (line_info.my_s_k2);

  MPI_MASTER (if (line_info.flag_warning > 10)
                  printf ("PowderN: %s: Warning messages were repeated %i times with absorbed xrays.\n", NAME_CURRENT_COMP, line_info.flag_warning);

              /* in case this instance is used in a SPLIT, we can recommend the
                 optimal iteration value */
              if (line_info.nb_refl_count) {
                double split_iterations = (double)line_info.nb_reuses / line_info.nb_refl_count + 1;
                double split_optimal = (double)line_info.nb_refl / line_info.nb_refl_count;
                if (split_optimal > split_iterations + 5)
                  printf ("PowderN: %s: INFO: You may highly improve the computation efficiency by using\n"
                          "    SPLIT %i COMPONENT %s=PowderN(...)\n"
                          "  in the instrument description %s.\n",
                          NAME_CURRENT_COMP, (int)split_optimal, NAME_CURRENT_COMP, instrument_source);
              });
%}

MCDISPLAY
%{
  if (line_info.V_0) {

    if (line_info.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 (line_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);
      if (line_info.zdepth_i) {
        xmin = -0.5 * line_info.xwidth_i;
        xmax = 0.5 * line_info.xwidth_i;
        ymin = -0.5 * line_info.yheight_i;
        ymax = 0.5 * line_info.yheight_i;
        zmin = -0.5 * line_info.zdepth_i;
        zmax = 0.5 * line_info.zdepth_i;
        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 (line_info.shape == 2) { /* sphere */
      if (line_info.radius_i) {
        circle ("xy", 0, 0, 0, line_info.radius_i);
        circle ("xz", 0, 0, 0, line_info.radius_i);
        circle ("yz", 0, 0, 0, line_info.radius_i);
      }
      circle ("xy", 0, 0, 0, radius);
      circle ("xz", 0, 0, 0, radius);
      circle ("yz", 0, 0, 0, radius);
    } else if (line_info.shape == 3) { /* OFF file */
      off_display (offdata);
    }
  }
%}
END