/*****************************************************************************
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2012 Risoe National Laboratory, Roskilde, Denmark

* Component: sample_nxs
*
* %IDENTIFICATION
* Author: Mirko Boin
* Date: Jan 2019
* Origin: Helmholtz-Zentrum Berlin fuer Materialien und Energie (Germany)
* Version: $Revision: 1.3$
*
* General powder/polycrystalline sample with neutron-matter interaction based
*    on neutron cross section calculations of a unit cell
*
* %DESCRIPTION
* Features:
*      - coherent and incoherent scattering, absorption and transmission
*      - multiple scattering 
*      - wavelength-dependent neutron cross section calculation
*      - unit cell definition (via input file)
*      - focussing option for scattered neutrons (via d_phi)
*      - strain gradient simulation for axial symmetric samples (...to be extended)
*
* Geometry is a powder filled cylinder or a box defined by radius or xwidth, 
* yheight, zthick respectively. The component handles coherent and incoherent
* scattering (also multiple scattering), absorption and transmission. Hence, it
* is suitable for diffraction (scattering) and imaging (transmission) instruments
* at the same time. In order to improve the neutron statistics at the detector
* these individual features can be enabled/disabled using TransOnly, IncohScat,
* MultiScat. 
* For example, if scattering shall not be monitored (e.g. imaging), then TransOnly
* can be set to 1 to let the component calculate the neutron transmission through
* the sample only. Likewise, incoherent scattering can be switched on and off.
* 
* The decision of whether a neuron shall be scattered, absorbed or transmitted is
* based on the wavelength-dependent neutron cross sections (nxs), i.e. the neutron
* velocity is used here. These cross sections represent the individual scattering 
* and absorption probabilities for each individual neutron. However, the calculation
* of neutron cross sections depends on the material defined for this sample component.
* Therefore, several input parameters such as the crystal structure and the atoms
* involved (isotope mass, coherent scattering length, ...). This information is
* provided by an input file with the following format:
*
* <table border="1"><tr align="center"><td><B>Al.nxs</B></td></tr><tr><td><B># define the unit cell parameters:<BR>
* #   space_group                      - the space group number or Hermann or Hall symbol [string]<BR>
* #   lattice_a, ...b, ...c            - the lattice spacings a,b,c [angstrom]<BR>
* #   lattice_alpha, ...beta, ...gamma - the lattice angles alpha,beta,gamma [degree]<BR>
* #   debye_temp                       - the Debye temperature [K]<BR>
* <BR>
* space_group    = -F 4 2 3 # space group number is also allowed (= 225)<BR>
* lattice_a = 4.049<BR>
* lattice_b = 4.049<BR>
* lattice_c = 4.049<BR>
* lattice_alpha = 90<BR>
* lattice_beta = 90<BR>
* lattice_gamma = 90<BR>
* debye_temp = 429.0<BR>
* <BR>
* # add atoms to the unit cell:<BR>
* # notation is "atom_number = name b_coh sigma_inc sigma_abs_2200 molar_mass x y z"<BR>
* #   name           - labels the current atom/isotope  [string]<BR>
* #   b_coh          - the coherent scattering length [fm]<BR>
* #   sigma_inc      - the incoherent scattering cross section [barns]<BR>
* #   sigma_abs_2200 - the absorption cross sect. at 2200 m/s [barns]<BR>
* #   molar_mass     - the Molar mass [g/mol]<BR>
* #   debye_temp     - the Debye temperature [K]<BR>
* #   x y z          - the Wyckoff postion of the atom inside the unit cell<BR>
* <BR>
* [atoms]<BR>
* add_atom = Al 3.449 0.008 0.23 26.98 0.0 0.0 0.0<BR>
* </B></td></tr></table>
* The above example defines a pure (fcc) aluminium. Other example files are provided.
* Alternatively, a fallback mechanism exists to simulate an alpha-Fe (bcc) sample
* if an input file is missing. But have a look for warning messages during simulation to
* make sure you are not running the fallback alpha-Fe simulation accidently! 
*
* CHANGELOG:
*    v1.0 - first official release with McStas 2.0
*    v1.1 - enabled strain distribution
*    v1.2 - use updated libnxs (global Debye temp -> parameter files changed)
*         - corrected sigma incoherent  
*         - updated towards McStas 2.1 conventions
*    v1.3 - bugfixes


* For details, please contact the author and visit the nxs website and/or the repository.
* NXS Website: <a href="http://www.helmholtz-berlin.de/people/boin/nxs_en.html">http://www.helmholtz-berlin.de/people/boin/nxs_en.html</a>
* NXS Repository: <a href="http://bitbucket.org/mirkoboin/nxs">http://bitbucket.org/mirkoboin/nxs</a>
*
*
* %PARAMETERS
*
* INPUT PARAMETERS
*
* TransOnly:   enable/disable (=1/=0) the transmission only option
* IncohScat:   enable/disable (=1/=0) incoherent scattering
* MultiScat:   enable/disable (=1/=0) multiple scattering
* xwidth:      horizontal dimension of sample, as a width [m]
* yheight:     vertical dimension of sample, as a height [m]
* zthick:      thickness of sample [m]
* radius:      radius of sample in (x,z) plane [m]
* nxsFileName: filename of the unit cell definition
* max_hkl:     maximum (hkl) indices to consider for the nxs calculation, e.g. max_hkl=3 means (hkl)=333
* d_phi=0:     Angle corresponding to the vertical angular range
*              to focus to, e.g. detector height. 0 for no focusing [deg,0-180]
* strainFileName:   filename for simulating strain gradients
* sample_temp: sample temperature [K]
*
* %EXAMPLES
* COMPONENT sample = sample_nxs( yheight = 0.005, radius = 0.0025, nxsFileName = "Al.nxs",
*                                max_hkl=8, TransOnly = 0, IncohScat = 1, MultiScat = 0, d_phi = 110.0 )
*
* %L
* M. Boin (2012), <i>J. Appl. Cryst.</i> <b>45</b>, 603-607, <a href="http://dx.doi.org/10.1107/S0021889812016056">doi:10.1107/S0021889812016056</a>
* %L
* M. Boin, R.C. Wimpory, A. Hilger, N. Kardjilov, S.Y. Zhang, M. Strobl (2012), <i>J. Phys.: Conf. Ser.</i> <b>340</b>, 012022, <a href="http://dx.doi.org/10.1088/1742-6596/340/1/012022">doi:10.1088/1742-6596/340/1/012022</a>
* %L
* M. Boin, A. Hilger, N. Kardjilov, S.Y. Zhang, E.C. Oliver, J.A. James, C. Randau, R.C. Wimpory (2011), <i>J. Appl. Cryst.</i> <b>44</b>, 1040-1046, <a href="http://dx.doi.org/10.1107/S0021889811025970">doi:10.1107/S0021889811025970</a>
*
*
* %END
*******************************************************************************/


DEFINE COMPONENT Sample_nxs
DEFINITION PARAMETERS ( int TransOnly=0, int IncohScat=1, int MultiScat=1 )
SETTING PARAMETERS ( xwidth=0, yheight=0, zthick=0, radius=0, thickness=0, string geometry=0, string nxsFileName="", int max_hkl=8, d_phi=0, string strainFileName="", sample_temp=293.0 )
OUTPUT PARAMETERS (isrect,intersect,nxs_init_success,p_transmit,xsect_total,xsect_coherent,
xsect_incoherent,xsect_absorption,
lambda,velocity,fullpath,t1,t2,uc,A,mu_factor,sd,V2L)
DEPENDENCY "-I@MCCODE_LIB@/libs/libnxs/ @MCCODE_LIB@/libs/libnxs/libnxs.a"



SHARE
%{


%include "nxs.h"
%include "read_table-lib"
%include "interoff-lib"
  
enum StrainDistribution {
  AxialSym1D,     /* 1D axial symmetric, i.e. cylinder sample with strain gradient along radius */
  AxialSym2D      /* 2D axial symmetric, i.e. cylinder sample with strain gradient along radius and along height */
  /*NoSym1D*/
};

typedef struct {
  /** values from FILE **/
  /* strains in the sample represented by lattice spacing */
  int n_entries;          /* number of entries for hoop, radial, axial */
  double *pos_x;          /* lists of all positions */
  double *pos_y;
  double *pos_z;
  double *strains_hoop;   /* lists of all strain values */
  double *strains_radial;
  double *strains_axial;
  double d0_hoop;         /* reference values */
  double d0_radial;
  double d0_axial;
  enum StrainDistribution strain_distr;  /* strain distribution mode */
} StrainData;

int readStrainFile( char *fileName, StrainData *sd )
{
  /* read the file */
  t_Table dataTable;
  Table_Read(&dataTable, fileName, 1); /* read 1st block data from file into table */
  int i;
  char **parsing;

  /* some default values to check if input file is consistent */
  int column_radius = 0;
  int column_height = 0;
  int column_hoop = 0;
  int column_radial = 0;
  int column_axial = 0;
  
  sd->d0_hoop = sd->d0_radial = sd->d0_axial = 0.0;
  
  /* parsing of header for sample parameters */
  parsing = Table_ParseHeader( dataTable.header,"d0_hoop","d0_radial","d0_axial",NULL );  
  if (parsing)
  {
    if( parsing[0] )
      sd->d0_hoop = atof(parsing[0]);
    else
    {
      fprintf(stderr, "%s: Warning: No value for 'd0_hoop' found in file %s\n", NAME_CURRENT_COMP, fileName);
    }
    if( parsing[1] )
      sd->d0_radial = atof(parsing[1]);
    else
    {
      fprintf(stderr, "%s: Warning: No value for 'd0_radial' found in file %s\n", NAME_CURRENT_COMP, fileName);
    }
    if( parsing[2] )
      sd->d0_axial = atof(parsing[2]);
    else
    {
      fprintf(stderr, "%s: Warning: No value for 'd0_axial' found in file %s\n", NAME_CURRENT_COMP, fileName);
    }
  }
  else
  {
    fprintf(stderr, "%s: Warning: parsing file %s\n", NAME_CURRENT_COMP, fileName);
    return -1;
  }

  /* parsing of header for column assignment */
  parsing = Table_ParseHeader(dataTable.header, "column_radius", "column_height","column_hoop","column_radial","column_axial",NULL);

  /* assign columns */
  if (parsing)
  {
    if (parsing[0]) column_radius = atoi(parsing[0]);
    if (parsing[1]) column_height = atoi(parsing[1]);
    if (parsing[2]) column_hoop = atoi(parsing[2]);
    if (parsing[3]) column_radial = atoi(parsing[3]);
    if (parsing[4]) column_axial = atoi(parsing[4]);
    for (i=0; i<5; i++)
      if (parsing[i])
        free(parsing[i]);
    free(parsing);
    
    /* decide which strain distribution is stored in the file */
    int minimum_columns = 4;
    if( column_radius ) 
    {
      sd->strain_distr = AxialSym1D;
      minimum_columns = 4;
    }
    if( column_radius && column_height )
    {
      sd->strain_distr = AxialSym2D;
      minimum_columns = 5;
    }
    
    /* check if data exists (at least minimum number of columns) */
    if( dataTable.columns < minimum_columns )
    {
      fprintf(stderr, "%s: Error: There are too few table columns in '%s' to process strain data - using default values now\n", NAME_CURRENT_COMP, fileName);
      sd->n_entries = 0;
      return -1;
    }
    

    /* allocate data array */
    sd->pos_x = (double*)malloc(dataTable.rows*sizeof(double));
    if( sd->strain_distr == AxialSym2D )
      sd->pos_y = (double*)malloc(dataTable.rows*sizeof(double));

    sd->strains_hoop   = (double*)malloc(dataTable.rows*sizeof(double));
    sd->strains_radial = (double*)malloc(dataTable.rows*sizeof(double));
    sd->strains_axial  = (double*)malloc(dataTable.rows*sizeof(double));
    sd->n_entries = dataTable.rows;

    /* get data from table */
    struct list_struct{
      double pos_x;
      double pos_y;
      double hoop;
      double radial;
      double axial;
    };
    struct list_struct* dataList = (struct list_struct*)malloc( sizeof(struct list_struct)*sd->n_entries );

    
    for (i=0; i<sd->n_entries; i++)
    {
      dataList[i].pos_x = Table_Index(dataTable, i, column_radius-1);
      if( sd->strain_distr == AxialSym2D )
        dataList[i].pos_y = Table_Index(dataTable, i, column_height-1);
      dataList[i].hoop = Table_Index(dataTable, i, column_hoop-1);
      dataList[i].radial = Table_Index(dataTable, i, column_radial-1);
      dataList[i].axial = Table_Index(dataTable, i, column_axial-1);
    }

    /* sort list by position */
    int pos_x_compare( const void *par1, const void *par2)
    {
      return ( ((struct list_struct*)par1)->pos_x<((struct list_struct*)par2)->pos_x ) ? 0 : 1;
    }
    /* (re-)sort the data set to evaluate it in the TRACE section */
    qsort( dataList, dataTable.rows, sizeof(struct list_struct), pos_x_compare );
    
    if( sd->strain_distr == AxialSym2D )
    {
      int pos_y_compare( const void *par1, const void *par2)
      {
        return ( ((struct list_struct*)par1)->pos_y<((struct list_struct*)par2)->pos_y ) ? 0 : 1;
      }
      qsort( dataList, dataTable.rows, sizeof(struct list_struct), pos_y_compare );
    }
    
    for (i=0; i<sd->n_entries; i++)
    {
      sd->pos_x[i] = dataList[i].pos_x;
      
      if( sd->strain_distr == AxialSym2D )
        sd->pos_y[i] = dataList[i].pos_y;
        
      sd->strains_hoop[i] = dataList[i].hoop;
      sd->strains_radial[i] = dataList[i].radial;
      sd->strains_axial[i] = dataList[i].axial;
    }
    
    /* data set is now sorted */
  }
  return 0;
}

%}




DECLARE
%{

int shape;
off_struct offdata;
int nxs_init_success;

/* cross sections */
double xsect_coherent;
double xsect_incoherent;
double xsect_absorption;

/* sample */
NXS_UnitCell uc;

StrainData sd;

/* velocity to lambda conversion */
double V2L;
%}




INITIALIZE
%{
fprintf(stderr, "%s: Using libnxs version %s\n", NAME_CURRENT_COMP, nxs_version());

shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape  */
if( geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0") )
{
  if( off_init(geometry, xwidth, yheight, zthick, 0, &offdata) )
  {
    shape = 3; /* arbitrary */
    thickness=0; //concentric=0;
  }
}
else if( xwidth && yheight && zthick )
  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,"Sample_nxs %s has invalid dimensions.\n"
    "ERROR: Please check parameter values (xwidth, yheight, zthick, radius).\n", NAME_CURRENT_COMP) );

  
if( readStrainFile( strainFileName, &sd ) )
{
  fprintf(stderr, "%s: No strain data found! Recalculation of d-values not used.\n", NAME_CURRENT_COMP);
  sd.n_entries = 0;
}

if (TransOnly)
{
  fprintf(stderr, "%s: Performing neutron transmission only (no scattering).\n", NAME_CURRENT_COMP);
  if( sd.n_entries>0 )
    fprintf(stderr, "%s: WARNING: Cannot use strain data in TransOnly mode!\n", NAME_CURRENT_COMP);

}
else
{
  if (IncohScat)
    fprintf(stderr, "%s: Enabling incoherent scattering.\n", NAME_CURRENT_COMP);
  if (MultiScat)
  {
    fprintf(stderr, "%s: Enabling multiple scattering.\n", NAME_CURRENT_COMP);
    if (d_phi)
    {
      fprintf(stderr, "%s: WARNING: No focussing possible in multiple scattering mode. Setting d_phi=0.\n", NAME_CURRENT_COMP);
      d_phi = 0;
    }
  }
}
  

/* read unit cell parameters from file and initialise hkl */
uc = nxs_newUnitCell();
NXS_AtomInfo *atomInfoList;
int numAtoms = nxs_readParameterFile( nxsFileName, &uc, &atomInfoList); 
if( numAtoms < 1 )
{
  /* fallback solution: if no file exists, use alpha_iron */
  fprintf(stderr, "%s: WARNING: nxs parameter file %s NOT found! Using default values...\n", NAME_CURRENT_COMP, nxsFileName);
  strncpy(uc.spaceGroup,"229",MAX_CHARS_SPACEGROUP);
  uc.a = 2.866; uc.alpha = 90.0; uc.debyeTemp = 464.0;
  NXS_AtomInfo ai;
  strncpy(ai.label,"Fe",MAX_CHARS_ATOMLABEL); ai.b_coherent = 9.45;
  ai.sigmaIncoherent = 0.4; ai.sigmaAbsorption = 2.56;
  ai.molarMass = 55.85;
  ai.x[0] = ai.y[0] = ai.z[0] = 0.0;
  atomInfoList = (NXS_AtomInfo*)realloc( atomInfoList, sizeof(NXS_AtomInfo) );
  atomInfoList[0] = ai;
  numAtoms = 1;
}

unsigned int i;
if( NXS_ERROR_OK != nxs_initUnitCell(&uc) )
{
  fprintf(stderr, "%s: WARNING: No nxs parameters set! Sample will be transparent!\n", NAME_CURRENT_COMP);
  nxs_init_success = 0;
}
else
{
  nxs_init_success = 1;
  uc.temperature = sample_temp;
  for( i=0; i<numAtoms; i++ )
    nxs_addAtomInfo( &uc, atomInfoList[i] );

  uc.maxHKL_index = max_hkl;
  nxs_initHKL( &uc );

  if( sd.n_entries>0 )
  {
    if( sd.d0_hoop<=1E-6 ) sd.d0_hoop = uc.a;
    if( sd.d0_radial<=1E-6 ) sd.d0_radial = uc.a;
    if( sd.d0_axial<=1E-6 ) sd.d0_axial = uc.a;
  }

  V2L = 2.0 * PI / V2K;
}

%}



TRACE
%{

int i;
double strain_axial;
double strain_hoop;
double strain_radial;
int intersect = 0;
double p_transmit;
double lambda;
double velocity;
double fullpath;
/* intersection times */
double t1;
double t2;

double xsect_total;


/* box or cylinder? */
if ( shape == 0 )
  intersect = cylinder_intersect(&t1,&t2,x,y,z,vx,vy,vz,radius, yheight);
else if( shape == 1 )
  intersect = box_intersect(&t1,&t2,x,y,z,vx,vy,vz,xwidth, yheight, zthick);
//else if( shape == 2 )
//  intersect = sphere_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius);
else if( shape == 3 )
  intersect = off_intersect(&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, offdata );

/* neutron intersects? */
if( intersect && (t2>0) && nxs_init_success )
{
  /* get current velocity and wavelength */
  velocity  = sqrt( vx*vx + vy*vy + vz*vz );
  lambda = V2L / velocity;

  xsect_coherent = nxs_CoherentElastic(lambda, &uc );
  xsect_incoherent = nxs_IncoherentElastic(lambda, &uc ) + nxs_IncoherentInelastic(lambda, &uc ) + nxs_CoherentInelastic(lambda, &uc );
  xsect_absorption = nxs_Absorption(lambda, &uc );
  xsect_total = xsect_coherent + xsect_incoherent + xsect_absorption;

  /* Handle transmission only (imaging mode) */
  if (TransOnly)
  {
    /* path through the sample */
    fullpath = velocity * (t2-t1);

    /* change the neutron weight */
    p_transmit = exp( -xsect_total * fullpath*1E2 / uc.volume );
    p *= p_transmit;
  }
  /* ...also handle scattering events */
  else
  {
    int ms_loop = 1;
    if( MultiScat ) ms_loop = 10;
    double path;

    while( ms_loop-- )
    {
      /* go to an event point (randomly) */
      path = rand01() * (t2-t1);
      PROP_DT( path + t1 );
      // printf("%i: t1=%f  t2=%f\n",ms_loop, t1*1000,t2*1000);
      path = path * velocity;
      
      /******** consider strain data... ********/
      if( sd.n_entries>0)
      {
        /* take care of neutron position to find the correct (micro)strain value */
        double pos = sqrt( x*x + z*z );
        i = 0;
        while( pos>sd.pos_x[i] && sd.n_entries>i )
        {
          i++;
        }
        if( i==0 )
        {
          strain_hoop = sd.strains_hoop[i];
          strain_radial = sd.strains_radial[i];
          strain_axial = sd.strains_axial[i];
          pos = sd.pos_x[i];
        }
        else if( i==sd.n_entries )
        {
          strain_hoop = sd.strains_hoop[i-1];
          strain_radial = sd.strains_radial[i-1];
          strain_axial = sd.strains_axial[i-1];
          pos = sd.pos_x[i-1];
        }
        else
        {
          /* get lattice strains from current neutron position */
          /* use linear interpolation between two given positions */
          strain_axial = (sd.strains_axial[i] - sd.strains_axial[i-1]) / (sd.pos_x[i] - sd.pos_x[i-1]) * (pos - sd.pos_x[i-1]) + sd.strains_axial[i-1];
          strain_hoop = (sd.strains_hoop[i] - sd.strains_hoop[i-1]) / (sd.pos_x[i] - sd.pos_x[i-1]) * (pos - sd.pos_x[i-1]) + sd.strains_hoop[i-1];
          strain_radial = (sd.strains_radial[i] - sd.strains_radial[i-1]) / (sd.pos_x[i] - sd.pos_x[i-1]) * (pos - sd.pos_x[i-1]) + sd.strains_radial[i-1];
        }
      
        /* now take care of the neutron flight direction */
        /* use sin2psi calculation to determine current dvalue */

        double sin2psi;
        double dvalue;
        sin2psi = sin(acos( (x*vx + z*vz) / sqrt(x*x + z*z) / sqrt(vx*vx + vz*vz) ));
        sin2psi = sin2psi * sin2psi;
        double d_hoop = strain_hoop*sd.d0_hoop*1E-6+sd.d0_hoop;
        double d_radial = strain_radial*sd.d0_radial*1E-6+sd.d0_radial;
        dvalue = (d_hoop-d_radial) * sin2psi + d_radial;

        /* and now axial also... */
        sin2psi = sin(acos( vy / sqrt(vx*vx + vy*vy + vz*vz) ));
        sin2psi = sin2psi * sin2psi;
        double d_axial = strain_axial*sd.d0_axial*1E-6+sd.d0_axial;
        dvalue = (dvalue-d_axial) * sin2psi + d_axial;

        /* re-init nxs calculation & adjust lattice spacings for all hkl...*/ 
        for( i=0; i<uc.nHKL; i++ )
        {
          uc.a = dvalue;
          NXS_HKL *hkl = &(uc.hklList[i]);
          hkl->dhkl = nxs_calcDhkl( hkl->h, hkl->k, hkl->l, &uc );
          hkl->FSquare = nxs_calcFSquare( hkl, &uc );
        }
        xsect_coherent = nxs_CoherentElastic(lambda, &uc );
        xsect_incoherent = nxs_IncoherentElastic(lambda, &uc ) + nxs_IncoherentInelastic(lambda, &uc ) + nxs_CoherentInelastic(lambda, &uc );
        xsect_absorption = nxs_Absorption(lambda, &uc );
        xsect_total = xsect_coherent + xsect_incoherent + xsect_absorption;
      }
      
      /* path through the sample */
      fullpath = velocity * (t2-t1);
      p_transmit = exp( -xsect_total * fullpath*1E2 / uc.volume );

      /* check if neutron interacts with or transmits through the sample */
      if( p_transmit < rand01() )
      {
        double roulette_ball = rand01() *  (xsect_total);
        double norm = lambda*lambda*1E-2 / 2.0 / uc.volume;
      
        /* ******************** */
        /*  SCATTER coherently  */
        /* ******************** */
        if (roulette_ball <= xsect_coherent)
        {
          int j;
          int max_hkl = -1;
          double contrib = 0.0;
          while( max_hkl<(int)uc.nHKL-1 && 2.0*uc.hklList[max_hkl+1].dhkl-lambda > 1E-6 )
          {
            max_hkl++;
            contrib += uc.hklList[max_hkl].FSquare * uc.hklList[max_hkl].multiplicity * uc.hklList[max_hkl].dhkl;
          }
          
          /* determine lattice plane (for scattering) */
          roulette_ball = rand01() * contrib;
          contrib = 0.0;
          for( j=0; j<max_hkl; j++ )
          {
            contrib += uc.hklList[j].FSquare * uc.hklList[j].multiplicity * uc.hklList[j].dhkl;
            if( roulette_ball < contrib )
              break;
          }
        
          /* get scattering angle */
          double theta = asin( lambda / 2.0 / uc.hklList[j].dhkl );
          if( isnan(theta) )
          {
            /* if rounding errors occur */
            theta = PI/2.0;
          }

          /* select random point (or within a smaller range given by d_phi) on Debye Scherrer cone */
          /* maximum d_phi = 180 */
          double tmp_vx, tmp_vy, tmp_vz;
          double vout_x, vout_y, vout_z;
          double arg, d_phi0;
          
          if (d_phi)
          { 
            arg = sin(d_phi*DEG2RAD/2)/sin(2*theta);
            if (arg < -1 || arg > 1)
              d_phi = 0;
            else
              d_phi = 2*asin(arg);
          }
          if (d_phi)
          {
            d_phi0 = 2*rand01()*fabs(d_phi);
            if (d_phi0 > d_phi)
            {
              d_phi0 = PI+(d_phi0-1.5*d_phi);
            }
            else
            {
              d_phi0=d_phi0-0.5*d_phi;
            }
            p *= d_phi/PI;
          }
          else
            d_phi0 = PI*randpm1();
          
          
          vec_prod(tmp_vx,tmp_vy,tmp_vz, vx,vy,vz, 1,0,0);
          if (!tmp_vx && !tmp_vy && !tmp_vz) { tmp_vx=tmp_vz=0; tmp_vy=1; }
          /* v_out = rotate 'v' by 2*theta around tmp_v: Bragg angle */
          rotate(vout_x,vout_y,vout_z, vx,vy,vz, 2.0*theta, tmp_vx,tmp_vy,tmp_vz);
          /* tmp_v = rotate v_out by d_phi around 'v' (Debye-Scherrer cone) */
          rotate(tmp_vx,tmp_vy,tmp_vz, vout_x,vout_y,vout_z, d_phi0, vx, vy, vz);
          vx = tmp_vx;
          vy = tmp_vy;
          vz = tmp_vz;
        }
        
        /* ******************** */
        /* SCATTER incoherently */
        /* ******************** */
        else if (roulette_ball <= xsect_coherent+xsect_incoherent)
        {
          /* check the incoherent switch */
          if (IncohScat)
          {
            double solid_angle;
            randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, 0, 0, 1, 2.0*PI, d_phi*DEG2RAD, ROT_A_CURRENT_COMP);
            
            vx *= velocity;
            vy *= velocity;
            vz *= velocity;
            if (d_phi)
              p *= d_phi/PI;

          } /* end of if (bIncohScat) */
        }
        else
        {
          /* neutron absorption -> remove neutron from trajectory */
          ms_loop = 0;
          ABSORB;
        } /* end of if-else( roulette_ball <= xsect_coherent ) */
        
        
        int err = 0;
        if( shape==0 && !cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius, yheight) || t2<0 )
          err=1;
        else if( shape==1 && !box_intersect(&t1, &t2, x, y, z, vx, vy, vz, xwidth, yheight, zthick) || t2<0 )
          err=1;
        else if( shape==2 && !sphere_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius) || t2<0 )
          err=1;
        else if( shape==3 && !off_intersect(&t1, &t2, NULL, NULL, x, y, z, vx, vy, vz, offdata) || t2<0 )
          err=1;

        if( err )
        {
          /* Strange error */
          fprintf(stderr, "sample_nxs: FATAL ERROR: Did not hit sample from inside.\n, t1=%f  t2=%f\n", t1, t2);
          ABSORB;
        }
        t1 = 0.0;
        SCATTER;
        
      } /* end of if( p_transmit < rand01() ) */
      else
      {
        /* else let the neutron simply transmit through the sample */
        /* without any interaction or neutron weight change */
        ms_loop = 0;
      }
      
    } /* end of while( ms_loop-- ) */
      
  } /* end of if-else(bTransOnly) */

} /* end of if(box_intersect) */



%}



MCDISPLAY
%{

magnify("xyz");
if( geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0") )
{	/* OFF file */
  off_display(offdata);
}
else if( radius > 0 &&  yheight )
{	/* cylinder */
  circle("xz", 0,  yheight/2.0, 0, radius);
  circle("xz", 0, -yheight/2.0, 0, radius);
  line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
  line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
  line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
  line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
  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( xwidth && yheight )
{ /* box/rectangle XY */
  double xmin = -0.5*xwidth;
  double xmax =  0.5*xwidth;
  double ymin = -0.5*yheight;
  double ymax =  0.5*yheight;
  double zmin = -0.5*zthick;
  double zmax =  0.5*zthick;
  multiline(5, xmin, ymin, zmin,
               xmax, ymin, zmin,
               xmax, ymax, zmin,
               xmin, ymax, zmin,
               xmin, ymin, zmin);
  multiline(5, xmin, ymin, zmax,
               xmax, ymin, zmax,
               xmax, ymax, zmax,
               xmin, ymax, zmax,
               xmin, ymin, zmax);
  line(xmin, ymin, zmin, xmin, ymin, zmax);
  line(xmax, ymin, zmin, xmax, ymin, zmax);
  line(xmin, ymax, zmin, xmin, ymax, zmax);
  line(xmax, ymax, zmin, xmax, ymax, zmax);
}
else if ( radius > 0 && !yheight )
{	/* sphere */
  circle("xy", 0,  0.0, 0, radius);
  circle("xz", 0,  0.0, 0, radius);
  circle("yz", 0,  0.0, 0, radius);        
}
  
%}

END