/*****************************************************************************
*
*  McStas, neutron ray-tracing package
*  Copyright(C) 2005 Risoe National Laboratory.
*
* Component: Isotropic_Sqw
*
* %I
* Written by: Virginie Hugouvieux, E. Farhi
* Date: August 2003
* Version:
* Origin:  ILL
* Modified by: E. Farhi, Jul 2005: made it work, concentric mode, multiple use
*
* Isotropic sample handling multiple scattering and absorption for a general
* S(q,w) (coherent and/or incoherent)
*
* %D
* An isotropic sample handling multiple scattering and including as input the
* dynamic structure factor of the chosen sample (e.g. from Molecular
* Dynamics). Handles elastic/inelastic, coherent and incoherent scattering -
* depending on the input S(Q,w) - with multiple scattering and absorption.
* Only the norm of Q is handled (not the vector), and thus suitable for
* liquids, gazes, amorphous and powder samples.
*
* If incoherent S(Q,w) file is specified as empty (0 or "") then the
* scattering is constant isotropic (Vanadium like).
* In case you only have one S(q,w) data containing both coherent and
* incoherent contributions you should e.g. use 'file_coh' and set 'sigma_coh'
* to the total scattering cross section.
*
* The S(q,w) are used as probability tables, and the intensity is weighted to
* the S(q,w) integral value. This integral value may be set from:
*   norm_sqw=-1 compute |S| from the density (refer to the Component Manual)
*   norm_sqw>0  use     |S| = norm_sqw
*   norm_sqw=0  compute |S| = \int S(q,w) dq dw from raw data files
* The default and recommanded normalization is 'norm_sqw=-1'.
*
* Focusing on the relevant [q,w] data range corresponding to the instrument
* setting can speed-up substantially the computation and prevent some neutrons
* to be removed when energy transfert is higher than actual neutron energy.
* An automatic such process may be activated when 'auto_qw=1' (recommanded).
* Additionally, for single order scattering (order=1), you may restrict the
* vertical spreading of the scattering area using d_phi parameter.
*
* If you use this component and produce valuable scientific results, please
* cite authors with references bellow (in Links).
*
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or a hollow cylinder/sphere.
*   box/plate:       xwidth x yheight x zthick (thickness=0)
*   hollow box/plate:xwidth x yheight x zthick and thickness>0
*   cylinder:        radius_o x yheight (radius_i=0)
*   hollow cylinder: radius_o x yheight and radius_i>0 or thickness>0
*   sphere:          radius_o (yheight=0 radius_i=0)
*   hollow sphere:   radius_o and radius_i>0 or thickness>0 (yheight=0)
*
* <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 S_in = Isotropic_Sqw(Sqw_coh="Al.laz", concentric=1, ...)
* AT (0,0,0) RELATIVE sample_position
*
* COMPONENT something_inside ... // e.g. the sample itself or other walls
* ...
*
* COMPONENT S_out = Isotropic_Sqw(Sqw_coh="Al.laz", ...)
* AT (0,0,0) RELATIVE sample_position
*
* <b>Sqw file format:</b>
* File format for S(Q,w) (coherent and incoherent) should contain 3 numerical
* blocks, defining q axis values (vector), then energy axis values (vector),
* then a matrix with one line per q axis value, containing Sqw values for
* each energy axis value. Comments (starting with '#') and non numerical lines
* are ignored and used to separate blocks. Sampling must be regular.
*
* Example:
* # q axis values
* # vector of m values in Angstroem-1
* 0.001000 .... 3.591000
* # w axis values
* # vector of n values in meV
* 0.001391 ... 1.681391
* # sqw values (one line per q axis value)
* # matrix of S(q,w) values (m rows x n values), one line per q value,
* 9.721422  10.599145 ... 0.000000
* 10.054191 11.025244 ... 0.000000
* ...
* 0.000000            ... 3.860253
*
* See for instance file He4.sqw. Such files may be obtained from e.g. INX,
* Nathan, Lamp and IDA softwares, as well as Molecular Dynamics.
*
* <b>Powder file format:</b>
* Files for coherent elastic powder scattering may also be used.
* Format specification follows the same principle as in the PowderN
* component, with parameters:
*
*     powder_format=Crystallographica
* or  powder_format=Fullprof
* or  powder_format=Lazy
* or  powder_format={j,d,F2,DW,Delta_d/d,1/2d,q,F} (column indexes 1:n)
*
* or column indexes (starting from 1) given as comments in the file header
* (e.g. '#column_j 4'). Refer to the PowderN component for more details.
* Delta_d/d and Debye-Waller factor may be specified for all lines with the
* 'powder_Dd' and 'powder_DW' parameters.
*
* Additionally a special [q,Sq] format is also defined with:
*   powder_format=qSq
* for which column 1 is 'q' and column 2 is 'S(q)'.
*
* <b>Examples:</b>
* 1- Vanadium-like incoherent elastic scattering
*   Isotropic_Sqw(V_rho=1/13.827,
*     sigma_abs=5.08, sigma_inc=4.935, sigma_coh=0)
*
* 2- liq-4He parameters
*   Isotropic_Sqw(..., Sqw_coh="He4.sqw", T=10)
*
* 3- powder sample
*  Isotropic_Sqw(..., Sqw_coh="Al.laz", save_sqw=1)
*
* %BUGS:
* When used in concentric mode, multiple bouncing scattering
* (traversing the hollow part) is not taken into account.
*
* %VALIDATION
* For Vanadium incoherent scattering mode, V_sample, Single_crystal and
* Isotropic_Sqw produce equivalent results, eventhough the two later are more
* accurate (geometry, multiple scattering).
*
* %P
* INPUT PARAMETERS:
* Sqw_coh:    [str] Name of the file containing the values of Q, w and S(Q,w)
*                     Coherent part; Q in Angs-1, E in meV, S(q,w) in meV-1
*                     Use 0 or "" to disable.
* Sqw_inc:    [str] Name of the file containing the values of Q, w and S(Q,w)
*                     Incoherent (self) part.
*                     Use 0 or "" to scatter isotropically (V-like).
* sigma_coh:[barns] Coherent Scattering cross-section. 0 to disable
* sigma_inc:[barns] Incoherent Scattering cross-section. 0 to disable
* sigma_abs:[barns] Absorption cross-section at 2200 m/s. 0 to disables
* V_rho:     [AA-3] Density of atoms (nb atoms/unit cell V_0).
* T:            [K] Temperature of sample, detailed balance
*
* Geometry parameters:
* radius_o:     [m] Outer radius of sample in (x,z) plane. cylinder/sphere.
* radius_i:     [m] Inner radius of sample in (x,z) plane. cylinder/sphere.
* yheight:      [m] Height of sample in vertical direction for box/cylinder
*                     shapes
* xwidth:       [m] width for a box sample shape
* zthick:       [m] thickness for a bulk box sample shape
* thickness:    [m] Thickness of hollow sample (overrides radius_i).
*
* OPTIONAL PARAMETERS:
* norm_sqw:     [1] When positive, the value is used as the integral of S(q,w)
*                   Null value uses computed integral.
*                   Any negative value will normalize on the density.
* auto_qw:      [1] When set to 1, the [q,w] range will automatically be tuned
*                   to optimal setting whenever required.
* concentric:   [1] indicates that this component in hollow cylinder geometry
*                   may contain other components.
*                   It should be then duplicated after the inside part.
*                   See description for an example.
* order:        [1] Limit multiple scattering up to given order
*                     0:all (default), 1:single, 2:double, ...
* verbose:      [1] Verbosity level (0:silent, 1:normal, 2:verbose, 3:debug).
* d_phi:      [deg] scattering vertical angular spreading (usually the heigh
*                     of the next component/detector). Use 0 for full space.
*                     This is only relevant for single scattering (order=1).
* save_sqw:     [1] When set to 1, saves S(q), DOS and S(q,w) as monitors.
* weight:   [g/mol] atomic/molecular weight of material
* density: [g/cm^3] density of material. V_rho=density/weight/1e24*N_A
* qmin:    [Angs-1] Minimum Q value to use in S(q,w).
* qmax:    [Angs-1] Maximum Q value to use in S(q,w).
* wmin:       [meV] Minimum Energy value to use in S(q,w).
* wmax:       [meV] Maximum Energy value to use in S(q,w).
* threshold:    [1] Value under which S(Q,w) is not accounted for.
*                     to set according to the S(Q,w) values ,i.e. not too low.
*
* POWDER ELASTIC SCATTERING PARAMETERS:
* powder_Dd:             [1] global Delta_d/d spreading, or 0 if ideal.
* powder_DW:             [1] global Debey-Waller factor, if not in |F2| or 1.
* powder_format: [no quotes] name or definition of column indexes in file
* powder_Vc:          [AA^3] volume of the unit cell
*
* OUTPUT PARAMETERS:
* VarSqw : internal structure containing many members/info
* VarSqw.dq wavevector transfert [Angs-1]
* VarSqw.dw energy transfert [Angs-1]
* VarSqw.type interaction type of event
*         'c' (coherent),  't' (transmitted)
*         'i' (incoherent) 'v' (isotropic incoherent, Vanadium-like).
* SCATTERED: order of multiple scattering
*
* %Links
* Hugouvieux V, Farhi E, Johnson MR, Virtual neutron scattering experiments, Physica B, 350 (2004) 151.
* %L
* Hugouvieux V, PhD, University of Montpellier II, France (2004).
* %L
* Cross sections for single elements: http://www.ncnr.nist.gov/resources/n-lengths/
* %L
* Cross sections for compounds:       http://www.ncnr.nist.gov/resources/sldcalc.html
* %L
* Example data file <a href="../data/He4.sqw">He4.sqw</a>
* %L
* The <a href="PowderN.html">PowderN</a> component.
* %L
* Web Elements                        http://www.webelements.com/
* %E
*****************************************************************************/

DEFINE COMPONENT Isotropic_Sqw
  DEFINITION PARAMETERS (Sqw_coh=0, Sqw_inc=0, powder_format=Undefined)
  SETTING PARAMETERS(radius_i=0,radius_o=0,thickness=0,
    xwidth=0, yheight=0, zthick=0,
    qmin=0, qmax=0, wmin=0, wmax=0, auto_qw=0,
    threshold=1e-10, int order=0, T=0, verbose=1, d_phi=0, int concentric=0,
    V_rho=0, sigma_abs=0, sigma_coh=0, sigma_inc=0, save_sqw=0,
    norm_sqw=-1, powder_Dd=0, powder_DW=0, powder_Vc=0, density=0, weight=0)
  OUTPUT PARAMETERS (VarSqw, columns)
  STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

/*****************************************************************************/
/*****************************************************************************/

SHARE
%{

#ifndef ISOTROPIC_SQW
#define ISOTROPIC_SQW

/* {j d F2 DW Dd inv2d q F} + { Sq if j == -1}*/
#ifndef Crystallographica
#define Crystallographica { 4,5,7,0,0,0,0,0 }
#define Fullprof          { 4,0,8,0,0,5,0,0 }
#define Undefined         { 0,0,0,0,0,0,0,0 }
#define Lazy              {17,6,0,0,0,0,0,13}
#endif
/* special case for [q,Sq] table */
#define qSq               {-1,0,0,0,0,0,1,0 }

%include "read_table-lib"

/* For the density of states */
struct Sqw_D_struct
{
  double omega;    /* omega value for the data block */
  double I;        /* intensity for the current DOS block */
  double cumul_I;  /* cumulated intensity (between 0 and 1) */
};

/* For the Q(w) probabilities */
struct Sqw_Q_struct
{
   double Q;           /* omega value for the data block */
   double proba;       /* normalized probability for the current DOS block */
   double cumul_proba; /* normalized cumulated probability */
};

struct Sqw_Data_struct /* contains normalized Sqw data for probabilities */
{
  struct Sqw_D_struct  *DOS;           /* P(w)  = density of states */
  struct Sqw_Q_struct **Q_proba_per_w; /* P(Q|w)= probability of each Q with w */

  long  *DOS_lookup;
  long **QW_lookup;
  double*Sq;  /* S(q) = \int_w S(q,w) */
  long   nb_q;
  long   nb_w; /* length of q and w vectors/axes */
  double q_min, q_max;
  double w_min, w_max;
  long   lookup_length;
  char   filename[80];
  double intensity;
  double intensity_total;
};

struct Sqw_sample_struct {
  char   compname[256];

  struct Sqw_Data_struct Data_inc;
  struct Sqw_Data_struct Data_coh;

  double s_abs, s_coh, s_inc;
  double my_s;
  double my_a_v;
  double rho;
  double T2E;
  double sqSE2K;
  int    maxloop;
  long   neutron_removed;
  long   neutron_enter;
  long   neutron_pmult;
  long   neutron_exit;
  double Dd, DWfactor;

  double sqw_threshold;

  char   verbose_output;
  char   shape; /* 0:cylinder, 1:box, 2:sphere */

  double dq, dw; /* q/w transfert */
  char   type; /* interaction type: c(coherent), i(incoherent),
                                  V(isotropic incoherent), t(transmitted) */
  double ki_x,ki_y,ki_z,kf_x,kf_y,kf_z;
  double ti, tf;
  double vi, vf;
  double ki, kf;

  double w_max, w_min, q_max, q_min;
  int    column_order[9]; /* column signification */
  double Temperature;
  double sqw_norm;
  char   sqw_auto;
};

#include <stdio.h>
#include <math.h>

void Sqw_Data_init(struct Sqw_Data_struct *Sqw_Data)
{
  Sqw_Data->nb_q         =0;
  Sqw_Data->nb_w         =0;
  Sqw_Data->q_min        =0;
  Sqw_Data->q_max        =0;
  Sqw_Data->w_min        =0;
  Sqw_Data->w_max        =0;
  Sqw_Data->lookup_length=100; /* length of lookup tables */
  Sqw_Data->intensity    =0;
  Sqw_Data->intensity_total=0;
  strcpy(Sqw_Data->filename, "");
  Sqw_Data->DOS          =NULL;
  Sqw_Data->Q_proba_per_w=NULL;
  Sqw_Data->DOS_lookup   =NULL;
  Sqw_Data->QW_lookup    =NULL;
  Sqw_Data->Sq           =NULL;
}

double Sqw_powder_gauss(double x, double mean, double rms) {
  return (exp(-((x)-(mean))*((x)-(mean))/(2*(rms)*(rms)))/(sqrt(2*PI)*(rms)));
}

/*****************************************************************************
* Sqw_read_PowderN: Read PowderN data files
*   Returns t_Table array or NULL in case of error
* Used in : Sqw_readfile (1)
*****************************************************************************/
t_Table *Sqw_read_PowderN(struct Sqw_sample_struct *Sqw, t_Table sqwTable)
{
  struct line_data
  {
    double F2;                  /* Value of structure factor */
    double q;                   /* Q vector */
    int j;                      /* Multiplicity */
    double DWfactor;            /* Debye-Waller factor */
    double w;                   /* Intrinsic line width */
  };
  struct line_data *list = NULL;
  double q_count=0, j_count=0, F2_count=0;
  int    mult_count  =0;
  double q_step      =FLT_MAX;
  long   size        =0;
  int    i, index;
  double q_min=0, q_max=0;
  char   flag=0;
  int    list_count=0;
  double q_step_cur;
  char   flag_qSq = 0;

  t_Table *retTable;

  flag_qSq = (Sqw->column_order[8]>0 && Sqw->column_order[6]>0);

  size = sqwTable.rows;
  Table_Info(sqwTable);
  if (Sqw->verbose_output > 0)
    printf("Isotropic_sqw: Converting %d powder lines from %s into S(q,w) data\n",
        size, sqwTable.filename);
  /* 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, DWfactor=0, F2=0, Sq=-1, q=0;
      int index;

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

      /* get data from table using columns {j d F2 DW Dd inv2d q} + { Sq }*/
      /* column indexes start at 1, thus need to substract 1 */
      if (Sqw->column_order[0]>0)
        j = Table_Index(sqwTable, i, Sqw->column_order[0]-1);
      if (Sqw->column_order[1]>0)
        d = Table_Index(sqwTable, i, Sqw->column_order[1]-1);
      if (Sqw->column_order[2]>0)
        F2 = Table_Index(sqwTable, i, Sqw->column_order[2]-1);
      if (Sqw->column_order[3]>0)
        DWfactor = Table_Index(sqwTable, i, Sqw->column_order[3]-1);
      if (Sqw->column_order[4]>0)
        w = Table_Index(sqwTable, i, Sqw->column_order[4]-1);
      if (Sqw->column_order[5]>0)  {
        d = Table_Index(sqwTable, i, Sqw->column_order[5]-1); if (d) d = 1/d/2; }
      if (Sqw->column_order[6]>0)
        q = Table_Index(sqwTable, i, Sqw->column_order[6]-1);
      if (Sqw->column_order[7]>0 && !F2)
        {F2= Table_Index(sqwTable, i, Sqw->column_order[7]-1); F2 *= F2;}

      if (Sqw->column_order[8]>0)
        Sq= Table_Index(sqwTable, i, Sqw->column_order[8]-1);

      if (q > 0 && Sq >= 0) F2 = Sq;
      if (d > 0 && q <= 0)  q = 2*PI/d;

      /* assign and check values */
      j = (j > 0 ? j : 0);
      DWfactor = (DWfactor > 0 ? DWfactor : 1);
      w = (w>0 ? w : 0);
      F2 = (F2 >= 0 ? F2 : 0);
      d = (q > 0 ? 2*PI/d : 0);
      if (j == 0 || d == 0 || q == 0) {
        printf("Isotropic_sqw: %s: line %i has invalid definition\n"
               "         (mult=0 or q=0 or d=0)\n", Sqw->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; /* or S(q) if flag_qSq */

      if (q_max < d) q_max = q;
      if (q_min > d) q_min = q;
      if (list_count > 1) {
        q_step_cur = fabs(list[list_count].q - list[list_count-1].q);
        if (q_step_cur > 0 && (!q_step || q_step_cur < q_step))
         q_step = q_step_cur;
      }

      /* 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[list_count-1].j
        || (mult_count == list[list_count].j && i == size-1)) {
          printf("Isotropic_Sqw: %s: Set multiplicity to 1 for lines [%i:%i]\n"
                  "         (d-spacing %g is duplicated %i times)\n",
            Sqw->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 */

  /* now builds new Table_Array to continue with Sqw_readfile */
  if (q_max == q_min || !q_step) return(NULL);
  if (!flag_qSq)
    size = 3*(q_max - q_min)/q_step; /* set a default of 3 q values per line */
  else size = list_count;

  if (Sqw->verbose_output > 0)
    printf("Isotropic_sqw: q range [%g:%g], creating %li elements vector\n",
        q_min, q_max, size);

  retTable = (t_Table*)calloc(4, sizeof(t_Table));
  if (!retTable) printf("Isotropic_Sqw: ERROR: Cannot allocate PowderN->Sqw table.\n");
  else {
    char *header;
    if (!Table_Init(&retTable[0], size, 1))
      { printf("Isotropic_Sqw: Cannot allocate q-axis from Powder lines.\n"); return(NULL); }
    if (!Table_Init(&retTable[1], 1, 1))
      { printf("Isotropic_Sqw: Cannot allocate w-axis from Powder lines.\n"); return(NULL); }
    if (!Table_Init(&retTable[2], size, 1))
      { printf("Isotropic_Sqw: Cannot allocate Sqw from Powder lines.\n"); return(NULL); }
    Table_Init(&retTable[3], 0,0);

    header = malloc(64); if (header)
    { retTable[0].header = header; strcpy(retTable[0].header, "q"); }
    header = malloc(64); if (header)
    { retTable[1].header = header; strcpy(retTable[1].header, "w"); }
    header = malloc(64); if (header)
    { retTable[2].header = header; strcpy(retTable[2].header, "Sqw"); }
    for (i=0; i < 4; i++) {
      retTable[i].array_length = 3;
      retTable[i].block_number = i+1;
    }
    if (!flag_qSq)
      for (i=0; i<size; i++)
        retTable[0].data[i]  = q_min + i*(q_max - q_min)/size;
      for (i=0; i<list_count; i++) { /* loop on each Bragg peak */
      double peak_qmin, peak_qmax,factor,q;
      if (list[i].w > 0 && !flag_qSq) {
        peak_qmin = list[i].q*(1 - list[i].w*3);
        peak_qmax = list[i].q*(1 + list[i].w*3);
      } else { /* Dirac peak, no width */
        peak_qmin = peak_qmax = list[i].q;
      }
      factor = 2*(2*PI)*(2*PI)*(2*PI)*(2*PI)*Sqw->rho
               *list[i].j*(list[i].DWfactor ? list[i].DWfactor : 1)
               /(Sqw->type == 'c' ? Sqw->s_coh : Sqw->s_inc)*list[i].F2;
      for (q=peak_qmin; q <= peak_qmax; q += q_step) {
        index = (long)floor(size*(q - q_min)/(q_max - q_min));
        if (index < 0) index=0;
        else if (index >= size) index = size-1;
        if (flag_qSq) {
          retTable[2].data[index] += list[i].F2;
          retTable[0].data[index]  = list[i].q;
        } else if (list[i].w <=0 || list[i].w*q < q_step) /* step function */
          retTable[2].data[index] += factor/q_step;
        else /* gaussian */
          retTable[2].data[index] += factor
                  * Sqw_powder_gauss(q, list[i].q, list[i].w*list[i].q);
      }
    } /* end for i */
    Table_Stat(&retTable[0]); Table_Stat(&retTable[1]); Table_Stat(&retTable[2]);
  }

  return(retTable);
} /* Sqw_read_PowderN */

/*****************************************************************************
* Sqw_readfile: Read Sqw data files
*   Returns Sqw_Data_struct or NULL in case of error
* Used in : Sqw_init (2)
*****************************************************************************/
struct Sqw_Data_struct *Sqw_readfile(
  struct Sqw_sample_struct *Sqw, char *file, struct Sqw_Data_struct *Sqw_Data)
{

  t_Table *Table_Array= NULL;
  t_Table *newTable;
  long     nblocks     = 0;

  t_Table  Sqw_restrict;

  long     i, j;
  char     flag=1; /* exit_flag */
  double   q, w, val, sum;
  double   mat_density=0, mat_weight=0;
  double   alpha=1;
  long    *DOS_lookup;
  long   **QW_lookup;
  char   **parsing;

  /* setup default */
  Sqw_Data_init(Sqw_Data);

  if (!file || !strlen(file)) return(Sqw_Data);
  Table_Array = Table_Read_Array(file, &nblocks);
  strncpy(Sqw_Data->filename, file, 80);
  if (!Table_Array) return(NULL);

  /* parsing of header */
  parsing = Table_ParseHeader(Table_Array[0].header,
    "Vc","V_0",
    "sigma_abs","sigma_a ",
    "sigma_inc","sigma_i ",
    "column_j", /* 6 */
    "column_d",
    "column_F2",
    "column_DW",
    "column_Dd",
    "column_inv2d", "column_1/2d", "column_sintheta_lambda",
    "column_q", /* 14 */
    "sigma_coh","sigma_c ",
    "Temperature",
    "column_Sq",
    "column_F ", /* 19 */
    "V_rho",
    "density",
    "weight",
    NULL);
  if (parsing) {
    if (parsing[0] && !Sqw->rho)   Sqw->rho    =1/atof(parsing[0]);
    if (parsing[1] && !Sqw->rho)   Sqw->rho    =1/atof(parsing[1]);
    if (parsing[2] && !Sqw->s_abs) Sqw->s_abs  =  atof(parsing[2]);
    if (parsing[3] && !Sqw->s_abs) Sqw->s_abs  =  atof(parsing[3]);
    if (parsing[4] && !Sqw->s_inc) Sqw->s_inc  =  atof(parsing[4]);
    if (parsing[5] && !Sqw->s_inc) Sqw->s_inc  =  atof(parsing[5]);
    if (parsing[6])                   Sqw->column_order[0]=atoi(parsing[6]);
    if (parsing[7])                   Sqw->column_order[1]=atoi(parsing[7]);
    if (parsing[8])                   Sqw->column_order[2]=atoi(parsing[8]);
    if (parsing[9])                   Sqw->column_order[3]=atoi(parsing[9]);
    if (parsing[10])                  Sqw->column_order[4]=atoi(parsing[10]);
    if (parsing[11])                  Sqw->column_order[5]=atoi(parsing[11]);
    if (parsing[12])                  Sqw->column_order[5]=atoi(parsing[12]);
    if (parsing[13])                  Sqw->column_order[5]=atoi(parsing[13]);
    if (parsing[14])                  Sqw->column_order[6]=atoi(parsing[14]);
    if (parsing[15] && !Sqw->s_coh)   Sqw->s_coh=atof(parsing[15]);
    if (parsing[16] && !Sqw->s_coh)   Sqw->s_coh=atof(parsing[16]);
    if (parsing[17] && !Sqw->Temperature) Sqw->Temperature=atof(parsing[17]);
    if (parsing[18])                  Sqw->column_order[8]=atoi(parsing[18]);
    if (parsing[19])                  Sqw->column_order[7]=atoi(parsing[19]);
    if (parsing[20] && !Sqw->rho)     Sqw->rho    =atof(parsing[20]);
    if (parsing[21] && !Sqw->rho)     mat_density =atof(parsing[21]);
    if (parsing[22] && !Sqw->rho)     mat_weight  =atof(parsing[22]);

    for (i=0; i<=22; i++) if (parsing[i]) free(parsing[i]);
    free(parsing);
  }

  if (!Sqw->rho && mat_density > 0 && mat_weight > 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 */
    Sqw->rho = mat_density/mat_weight/1e24*6.02214199e23;
    if (Sqw->verbose_output) printf("Isotropic_Sqw: %s: Computing atom density V_rho=%g [AA^-3] \n",NAME_CURRENT_COMP, Sqw->rho);
  }

  if (nblocks) {
    if (nblocks == 1) {
      /* import Powder file */
      newTable = Sqw_read_PowderN(Sqw, Table_Array[0]);
      if (!newTable) {
        printf("Isotropic_Sqw: %s: ERROR importing powder line file %s.\n"
               "               Check format definition.\n",
              Sqw->compname, file);
        exit(-1);
      } else flag=0;
      Table_Free_Array(Table_Array);
      Table_Array = newTable;
    } else if (nblocks != 3) {
      printf("Isotropic_Sqw: %s: "
             "File %s contains %i block%s instead of 3.\n",
              Sqw->compname, file, nblocks, (nblocks == 1 ? "" : "s"));
    } else flag=0;
  }

  /* print some info about Sqw files */
  if (flag) Sqw->verbose_output = 2;
  if (nblocks && Sqw->verbose_output > 1) {
    printf("Isotropic_Sqw: %s file read, analysing...\n", file);
    Table_Info_Array(Table_Array);
  }
  if (flag) {
    if (nblocks) printf("ERROR          Wrong file format.\n"
      "               Disabling contribution.\n"
      "               File must contain 3 blocks for [q,w,sqw]\n");
    return(Sqw_Data);
  }

  /* restrict gridding range to q and omega min/max */
  for (i=0; i< 2; i++) { /* i=[q,w] */
    long    index;
    double  val_min, val_max;
    double *usr_min,*usr_max;
    t_Table Table;

    Table  = Table_Array[i];
    /* make q,w axes a single row */
    if (Table.columns == 1) { Table.columns = Table.rows; Table.rows = 1; }
    index  = Table.rows*Table.columns; /* full length of q,w */
    val_min= Table.min_x;
    val_max= Table.max_x;

    /* update min/max accroding to input parameters */
    if (i == 0) { /* q */
      val_min = (val_min < Sqw->q_min ? Sqw->q_min : val_min);
      val_max = (val_max > Sqw->q_max ? Sqw->q_max : val_max);
    } else {
      val_min = (val_min < Sqw->w_min ? Sqw->w_min : val_min);
      val_max = (val_max > Sqw->w_max ? Sqw->w_max : val_max);
    }

    if (i == 0) Sqw_Data->nb_q = index; else Sqw_Data->nb_w = index;

    /* set restricted range */
    usr_min = (i == 0 ? &(Sqw_Data->q_min) : &(Sqw_Data->w_min));
    usr_max = (i == 0 ? &(Sqw_Data->q_max) : &(Sqw_Data->w_max));
    if (!*usr_min && !*usr_max) {
      *usr_min = val_min; *usr_max = val_max;
    }
  }

  if (!Sqw_Data->nb_q  || !Sqw_Data->nb_w) {
    printf("Isotropic_Sqw: %s: Data file %s has null q or omega information (%lix%li).\n"
           "ERROR          Exiting.\n",
      Sqw->compname, file, Sqw_Data->nb_w, Sqw_Data->nb_q);
    return(NULL);
  }

  /* build the 'total' = Sqw table in restricted range */
  Table_Init(&Sqw_restrict, Sqw_Data->nb_q, Sqw_Data->nb_w);
  Sqw_restrict.block_number = 1;

  if (!Sqw_restrict.data || !Sqw_restrict.rows*Sqw_restrict.columns) {
    printf("Isotropic_Sqw: %s: Cannot allocate Sqw_restrict Table (%lix%li).\n"
           "ERROR          Exiting.\n",
      Sqw->compname, Sqw_Data->nb_q, Sqw_Data->nb_w);
    return(NULL);
  }
  sprintf(Sqw_restrict.filename, "S(q,w) from %s", file);

  /* now compute S(q) and normalize */
  Sqw_Data->Sq =
    (double *)calloc(Sqw_Data->nb_q, sizeof(double));
  if (!Sqw_Data->Sq) {
    printf("Isotropic_Sqw: %s: Cannot allocate S(Q) array (%li bytes).\n"
           "ERROR          Exiting.\n",
      Sqw->compname, Sqw_Data->nb_q*sizeof(double));
    return(NULL);
  }

  sum = val = alpha = 0;
  for (i=0; i < (Sqw_Data->nb_q) ; i++) {
    double sqw=0, sq=0;
    q  = Table_Index(Table_Array[0], 0, i);
    for (j=0; j < (Sqw_Data->nb_w) ; j++) {
      flag = 1;
      if (nblocks) sqw = Table_Index(Table_Array[2], i, j);

      if (sqw < Sqw->sqw_threshold) flag = 0; /* Sqw lower than threshold */
      /* check for q,w within min/max */
      w  = Table_Index(Table_Array[1], 0, j);
      if (q < Sqw_Data->q_min || q > Sqw_Data->q_max) flag = 0; /* within q range  ? */
      if (w < Sqw_Data->w_min || w > Sqw_Data->w_max) flag = 0; /* within w range  ? */
      if (flag) val += sqw; /* only on restricted range */
      sum += sqw;           /* on total Sqw range */
      if (!Table_SetElement(&Sqw_restrict, i, j, (flag ? sqw : 0) ))
        printf("Isotropic_Sqw: %s: "
               "Error when setting Sqw[%i,%i]=%g from file %s\n",
               Sqw->compname, i,j, sqw, file);
      sq += sqw;
    }
    Sqw_Data->Sq[i] = sq;
    alpha += q*q*sq;
  }
  if (Sqw->verbose_output > 1)
    printf("Isotropic_Sqw: %s: Generated S(Q)[%i]\n",
      Sqw->compname, Sqw_Data->nb_q);
  alpha /= Sqw_Data->nb_q; /* this is \int q^2 S(q) dq */

  if (alpha) { /* compute theoretical |S| norm */
    alpha = (Sqw_Data->q_max*Sqw_Data->q_max*Sqw_Data->q_max/3 - 2*PI*PI*Sqw->rho)/alpha;
    val *= alpha;
    sum *= alpha;
    for (i=0; i < (Sqw_Data->nb_q) ; Sqw_Data->Sq[i++] *= alpha);
  }

  val /= Sqw_Data->nb_q*Sqw_Data->nb_w;
  sum /= Sqw_Data->nb_q*Sqw_Data->nb_w;

  if (!alpha) {
    printf("Isotropic_Sqw: %s: Can not compute S(q) norm\n"
           "WARNING        Using raw norm |S|=%g (not safe)\n",
      Sqw->compname, sum);
   alpha=1;
  } else {
    if (Sqw->verbose_output > 1
        || fabs(log10(fabs(alpha))) > 1 )
      printf("Isotropic_Sqw: %s: Theoretical %s Norm |S| = %g\n",
                        Sqw->compname,
                        (Sqw->type=='c' ? "coherent" : "incoherent"),
                        sum);
    if (Sqw->sqw_norm>0 && fabs(log10(fabs(sum/Sqw->sqw_norm))) > 0.1)
      printf("WARNING        The specified norm %g is probably wrong\n"
             "               Multiply intensistes by factor %g\n",
        Sqw->sqw_norm, sum/Sqw->sqw_norm);
  }

  if (!val || !sum) {
    printf("Isotropic_Sqw: %s: selected (Q,w) range is outside\n"
           "ERROR          available Sqw data. |S|=%g\n",
      Sqw->compname, sum);
    printf("                 q=[%g:%g] w=[%g:%g] size=[%ix%i]\n",
           Sqw_Data->q_min, Sqw_Data->q_max,
           Sqw_Data->w_min, Sqw_Data->w_max,
           Sqw_Data->nb_q, Sqw_Data->nb_w);
    return(NULL);
  }

  Table_Stat(&Sqw_restrict);

  Sqw_Data->intensity       = val;
  Sqw_Data->intensity_total = sum;

  if (Sqw->verbose_output > 0) {
    printf("Isotropic_Sqw: %s: Generated %s normalized %scoherent Sqw\n"
           "                   q=[%g:%g] w=[%g:%g] |S|=%g size=[%ix%i]\n",
           Sqw->compname, file, (Sqw->type == 'i' ? "in" : ""),
           Sqw_Data->q_min, Sqw_Data->q_max,
           Sqw_Data->w_min, Sqw_Data->w_max, val,
           Sqw_Data->nb_q, Sqw_Data->nb_w);
  if (Sqw_Data->w_min == Sqw_Data->w_max && Sqw_Data->nb_w == 1)
    printf("               Elastic scattering only.\n");
  }

  /* set up density of states and index_wconst */
  /* uses: Sqw_restrict and w axes */

  Sqw_Data->DOS =
    (struct Sqw_D_struct*)calloc(Sqw_Data->nb_w, sizeof(struct Sqw_D_struct));

  if (!Sqw_Data->DOS) {
    printf("Isotropic_Sqw: %s: Cannot allocate DOS (%li bytes).\n"
           "ERROR          Exiting.\n",
      Sqw->compname, Sqw_Data->nb_w*sizeof(struct Sqw_D_struct));
    return(NULL);
  }
  val = 0;
  for (j=0; j < (Sqw_Data->nb_w) ; j++) {
    double local_val=0;
    w  = Table_Index(Table_Array[1], 0, j);
    for (i=0; i < (Sqw_Data->nb_q) ; i++) { /* integrate on all q values */
      local_val += Table_Index(Sqw_restrict, i, j)*alpha;
    }
    Sqw_Data->DOS[j].omega = w;
    Sqw_Data->DOS[j].I     = local_val;
    val                   += local_val; /* total intensity */
  }
  /* normnalize and compute cumulated probability */
  for (j=0; j < (Sqw_Data->nb_w) ; j++) {
    Sqw_Data->DOS[j].I      /= val;
    Sqw_Data->DOS[j].cumul_I = Sqw_Data->DOS[j].I;
    /* cumulated probabilities for the random choice */
    if (j) Sqw_Data->DOS[j].cumul_I += Sqw_Data->DOS[j-1].cumul_I;
  }
  if (Sqw->verbose_output > 1)
    printf("Isotropic_Sqw: %s: Generated normalized DOS[%i]\n",
      Sqw->compname, Sqw_Data->nb_w);

  /* set up Q probability table per w bin */
  /* uses:  Sqw_restrict */
  Sqw_Data->Q_proba_per_w =
    (struct Sqw_Q_struct**)calloc(Sqw_Data->nb_w, sizeof(struct Sqw_Q_struct*));

  if (!Sqw_Data->Q_proba_per_w) {
    printf("Isotropic_Sqw: %s: Cannot allocate P(Q|w) array (%li bytes).\n"
           "ERROR          Exiting.\n",
      Sqw->compname, Sqw_Data->nb_w*sizeof(struct Sqw_Q_struct*));
    return(NULL);
  }
  for (j=0; j < (Sqw_Data->nb_w) ; j++) {

    double intensity=0;

    Sqw_Data->Q_proba_per_w[j]=
        (struct Sqw_Q_struct*)calloc(Sqw_Data->nb_q, sizeof(struct Sqw_Q_struct));

    if (!Sqw_Data->Q_proba_per_w[j]) {
      printf("Isotropic_Sqw: %s: Cannot allocate P(Q|w)[%i] (%li bytes).\n"
             "ERROR          Exiting.\n",
        Sqw->compname, j, Sqw_Data->nb_q*sizeof(struct Sqw_Q_struct));
      return(NULL);
    }
    /* set P(Q|W) and compute total intensity */
    for (i=0; i < (Sqw_Data->nb_q) ; i++) {
      q  = Table_Index(Table_Array[0], 0, i);
      Sqw_Data->Q_proba_per_w[j][i].Q     = q;
      Sqw_Data->Q_proba_per_w[j][i].proba = Table_Index(Sqw_restrict, i, j)*alpha;
      intensity                     += Sqw_Data->Q_proba_per_w[j][i].proba;
    }
    for (i=0; i < (Sqw_Data->nb_q) ; i++) {
      Sqw_Data->Q_proba_per_w[j][i].cumul_proba =
                        Sqw_Data->Q_proba_per_w[j][i].proba/intensity;
      if (i>0)     /* cumulate */
        Sqw_Data->Q_proba_per_w[j][i].cumul_proba +=
          Sqw_Data->Q_proba_per_w[j][i-1].cumul_proba;
    }
  }
  if (Sqw->verbose_output > 1)
    printf("Isotropic_Sqw: %s: Generated P(Q|w)\n",
      Sqw->compname);

  /* generate quick lookup tables for DOS and Q_proba_per_w */
  Sqw_Data->DOS_lookup = NULL;
  DOS_lookup = (long*)calloc(Sqw_Data->lookup_length, sizeof(long));
  if (!DOS_lookup) {
    printf("Isotropic_Sqw: %s: Cannot allocate DOS_lookup (%li bytes).\n"
           "Warning        Will be slower.\n",
      Sqw->compname, Sqw_Data->lookup_length*sizeof(long));
  }
  if (DOS_lookup) {
    for (i=0; i < Sqw_Data->lookup_length; i++) {
      w = (double)i/(double)Sqw_Data->lookup_length;
      DOS_lookup[i] = Sqw_search_DOS(*Sqw_Data, w);
    }
    Sqw_Data->DOS_lookup = DOS_lookup;
  }
  Sqw_Data->QW_lookup = NULL;
  QW_lookup=(long**)calloc(Sqw_Data->nb_w, sizeof(long*));
  if (!QW_lookup) {
    printf("Isotropic_Sqw: %s: Cannot allocate QW_lookup (%li bytes).\n"
           "Warning        Will be slower.\n",
      Sqw->compname, Sqw_Data->nb_w*sizeof(long*));
  }
  if (QW_lookup) {
    for (j=0; j < (Sqw_Data->nb_w) ; j++) {
      QW_lookup[j] = (long*)calloc(Sqw_Data->lookup_length, sizeof(long));
      if (!QW_lookup[j]) {
        printf("Isotropic_Sqw: %s: Cannot allocate QW_lookup[%i] (%li bytes).\n"
               "Warning        Will be slower.\n",
        Sqw->compname, j, Sqw_Data->lookup_length*sizeof(long));
        free(QW_lookup); QW_lookup = NULL;
      } else {
        for (i=0; i < Sqw_Data->lookup_length; i++) {
          w = (double)i/(double)Sqw_Data->lookup_length;
          QW_lookup[j][i] = Sqw_search_Q_proba_per_w(*Sqw_Data, w, j);
        }
      }
    }
    Sqw_Data->QW_lookup = QW_lookup;
  }
  if ((Sqw_Data->QW_lookup || Sqw_Data->DOS_lookup) && Sqw->verbose_output > 1)
    printf("Isotropic_Sqw: %s: Generated lookup tables with %i entries\n",
      Sqw->compname, Sqw_Data->lookup_length);

  Table_Free_Array(Table_Array);
  Table_Free(&Sqw_restrict);

  return(Sqw_Data);
} /* end Sqw_readfile */

/*****************************************************************************
* Sqw_init_read: Read coherent/incoherent Sqw data files
*   Returns Sqw total intensity, or 0 (error)
* Used in : INITIALIZE (1)
*****************************************************************************/
double Sqw_init(struct Sqw_sample_struct *Sqw, char *file_coh, char *file_inc)
{
  double ret=0;

  struct Sqw_Data_struct *d_inc, *d_coh;
  Sqw->type = 'i';
  d_inc = Sqw_readfile(Sqw, file_inc, &(Sqw->Data_inc));
  if (d_inc && !d_inc->intensity && Sqw->s_inc) {
    if (Sqw->verbose_output > 0)
      printf("Isotropic_Sqw: %s: Using Isotropic elastic incoherent scattering\n"
             "               (Vanadium like)\n", Sqw->compname);
    ret=1;
  }
  Sqw->type = 'c';
  d_coh = Sqw_readfile(Sqw, file_coh, &(Sqw->Data_coh));
  if (!d_inc || !d_coh) return(0);
  if (d_coh && !d_coh->intensity && Sqw->s_coh)
    printf("Isotropic_Sqw: %s: Coherent scattering Sqw intensity is null.\n"
           "Warning        Disabling coherent scattering.\n", Sqw->compname);
  if (d_inc && d_coh && d_inc->intensity && d_coh->intensity) {
    char msg[80];
    strcpy(msg, "");
    /* check dimensions/limits for Q, Energy in coh and inc Tables */
    if (d_inc->nb_q  != d_coh->nb_q)
      strcpy(msg, "Q axis size");
    if (d_inc->nb_w  != d_coh->nb_w)
      strcpy(msg, "Energy axis size");
    if ((d_inc->q_min != d_coh->q_min)
     || (d_inc->q_max != d_coh->q_max))
      strcpy(msg, "Q axis limits");
    if ((d_inc->w_min != d_coh->w_min)
     || (d_inc->w_max != d_coh->w_max))
      strcpy(msg, "Energy axis limits");
    if (strlen(msg)) {
      printf("Isotropic_Sqw: %s: Sqw data from files %s and %s do not match\n"
             "ERROR          wrong %s\n",
             Sqw->compname, file_coh, file_inc, msg);
      return(0);
    }
  }
  if (!ret) ret=d_inc->intensity+d_coh->intensity;
  return(ret);
}

/*****************************************************************************
*  Sqw_search_root: Search for the roots of A*x**2 + B*x + C
*   The roots are returned as x1 and x2
* Used in : TRACE (1)
*****************************************************************************/
int Sqw_search_root(double A, double B, double C,
             double *x1, double *x2)

{

  double delta, sqrtdelta, inv_2A;
  int    OK = 0;

  if (fabs(A) < 1E-10) {
    if (B) {  *x1 = -C/B; *x2 = *x1; OK=3; }
  } else {
    delta = B*B - 4*A*C;
    if (delta < 0) {
      /* Imaginary roots */
      OK=0;
    } else {
      /* Real roots */
      sqrtdelta = sqrt(delta);
      inv_2A = 1/(2*A);
      *x1 = inv_2A*(-B - sqrtdelta);
      *x2 = inv_2A*(-B + sqrtdelta);
      OK = 1;
    }
    return(OK);
  }
}

/*****************************************************************************
*  Sqw_search_DOS: For a given random number 'randnum', search for the bin
*   containing  the corresponding Sqw->DOS
* Used in : TRACE (1)
*****************************************************************************/
int Sqw_search_DOS(struct Sqw_Data_struct Sqw, double randnum)
{
  int i=1;

  if (Sqw.DOS_lookup) {
    i = Sqw.DOS_lookup[(long)floor(randnum*Sqw.lookup_length)];
  }

  while (i < Sqw.nb_w && (&(Sqw.DOS[i]) != NULL) && (randnum > Sqw.DOS[i].cumul_I))
      i++;
  if (i >= Sqw.nb_w) i = Sqw.nb_w;

  if (&(Sqw.DOS[i]) == NULL)
  {
      fprintf(stderr, "Isotropic_Sqw: No corresponding value in the DOS. randnum too big.\n");
      fprintf(stderr, "  i=%i ; randnum=%f ; Sqw.DOS[i-1].cumul_I=%f (Sqw_search_DOS)\n",
            i, randnum, Sqw.DOS[i-1].cumul_I);
      return i-1;
  }
  else
      return (i < Sqw.nb_w ? i : Sqw.nb_w-1);
}

/*****************************************************************************
*  Sqw_search_Q_proba_per_w: For a given random number randnum, search for
*   the bin containing the corresponding Sqw.DOS in the Q probablility grid
* Used in : TRACE (1)
*****************************************************************************/
int Sqw_search_Q_proba_per_w(struct Sqw_Data_struct Sqw,
                             double randnum, int index)
{
  int i=0;

  if (Sqw.QW_lookup && Sqw.QW_lookup[index]) {
    i = Sqw.QW_lookup[index][(long)floor(randnum*Sqw.lookup_length)];
  }

  while (i < Sqw.nb_q && (&(Sqw.Q_proba_per_w[index][i]) != NULL)
    && (randnum > Sqw.Q_proba_per_w[index][i].cumul_proba)) {
      i++;
  }

  if (&(Sqw.Q_proba_per_w[index][i]) == NULL)
    return -1;
  else
    return i;
}

#endif

%}


/*****************************************************************************/
/*****************************************************************************/

DECLARE
%{

  struct Sqw_sample_struct VarSqw;
  int    columns[8] = powder_format;

/* end DECLARE */
%}

/*****************************************************************************/
/*****************************************************************************/

INITIALIZE
%{
  int i;
  /* check for parameters */

  VarSqw.verbose_output= verbose;
  if (radius_o>0) {
    if (thickness>0 && thickness < radius_o) radius_i=radius_o-thickness;
    if (yheight>0) VarSqw.shape = 0; else VarSqw.shape = 2;
    if (VarSqw.verbose_output)
      printf("Isotropic_Sqw: %s: is a %s%s : R_i=%f R_o=%f h=%f \n",
            NAME_CURRENT_COMP, (radius_i > 0 ? "hollow ":""),
            (VarSqw.shape ? "Sphere" : "Cylinder"),
            radius_i,radius_o,yheight);
    if (radius_i >= radius_o) {
      printf("Isotropic_Sqw: %s:thickness/radius_i bigger than sample (radius_o) !\n"
             "ERROR          Disabling\n",
           NAME_CURRENT_COMP);
      exit(0);
    }
  } else if (xwidth>0 && yheight>0 && zthick>0) {
    VarSqw.shape = 1;
    if (2*thickness>=xwidth || 2*thickness>=yheight || 2*thickness>zthick) {
      printf("Isotropic_Sqw: %s:2*thickness bigger than sample (xwidth yheight zthick)!\n"
             "ERROR          Disabling\n",
           NAME_CURRENT_COMP);
      exit(0);
    }

    if (VarSqw.verbose_output)
      printf("Isotropic_Sqw: %s: is a %sBox : dx=%f dy=%f dz=%f \n",
            NAME_CURRENT_COMP, (thickness>0 ? "hollow" : ""), xwidth,yheight,zthick);
  } else  {
    printf("Isotropic_Sqw: %s: ERROR sample has no dimension\n"
           "               set some of radius_o xwidth yheight zthick \n",
           NAME_CURRENT_COMP);
    exit(0);
  }
  if (concentric &&
  !((VarSqw.shape==1 && thickness > 0) || (VarSqw.shape!=1 && radius_i > 0))) {
    printf("Isotropic_Sqw: %s:Can not use concentric mode\n"
           "ERROR          on non hollow shape.\n",
           NAME_CURRENT_COMP);
    exit(0);
  }

  strncpy(VarSqw.compname, NAME_CURRENT_COMP, 256);
  VarSqw.T2E       =(1/11.605);   /* Kelvin to meV */
  VarSqw.sqSE2K = (V2K*SE2V)*(V2K*SE2V);
  VarSqw.sqw_threshold = (threshold > 0 ? threshold : 0);
  VarSqw.s_abs     = sigma_abs;
  VarSqw.s_coh     = sigma_coh;
  VarSqw.s_inc     = sigma_inc; /* s_scatt member initialized in Sqw_init */
  VarSqw.maxloop   = 100;
  VarSqw.neutron_removed = 0;
  VarSqw.neutron_enter   = 0;
  VarSqw.neutron_pmult   = 0;
  VarSqw.neutron_exit    = 0;
  VarSqw.w_max     = (wmax > 0 ? wmax : FLT_MAX);
  VarSqw.w_min     = (wmin > 0 ? wmin : 0);
  VarSqw.q_max     = (qmax > 0 ? qmax : FLT_MAX);
  VarSqw.q_min     = (qmin > 0 ? qmin : 0);
  VarSqw.rho       = V_rho;
  VarSqw.sqw_norm  = norm_sqw;
  VarSqw.sqw_auto  = auto_qw;

  /* PowderN compatibility members */
  VarSqw.Dd        = powder_Dd;
  VarSqw.DWfactor  = powder_DW;
  VarSqw.Temperature= T;
  for (i=0; i< 9; i++) VarSqw.column_order[i] = columns[i];
  VarSqw.column_order[8] = (VarSqw.column_order[0] >= 0 ? 0 : 2);

  if (!Sqw_init(&VarSqw, Sqw_coh, Sqw_inc)) {
    printf("Isotropic_Sqw: %s: ERROR importing data files (Sqw_init)\n",NAME_CURRENT_COMP);
    exit(0);
  }

  /* optional ways to define rho */
  if (!VarSqw.rho && powder_Vc > 0)
    VarSqw.rho = 1/powder_Vc;

  if (VarSqw.rho <= 0) exit(printf("Isotropic_Sqw: %s: Null density (V_rho)\n",NAME_CURRENT_COMP));
  /* 100: convert from barns to fm^2 */
  VarSqw.my_a_v  =(VarSqw.rho*100*VarSqw.s_abs*2200);
  VarSqw.my_s    =(VarSqw.rho*100*(VarSqw.s_coh+VarSqw.s_inc));

  if ((VarSqw.s_coh || VarSqw.s_inc) && !VarSqw.Temperature
   && (VarSqw.Data_coh.intensity || VarSqw.Data_inc.intensity))
    printf("Isotropic_Sqw: %s: Sample temperature is zero (T).\n"
           "Warning        Disabling detailed balance.\n",NAME_CURRENT_COMP);

  if (!VarSqw.s_coh && !VarSqw.s_inc) {
    printf("Isotropic_Sqw: %s: Scattering cross section is zero\n"
           "ERROR          (sigma_coh, sigma_inc)\n",NAME_CURRENT_COMP);
    exit(0);
  }

  if (d_phi) d_phi = fabs(d_phi)*DEG2RAD;

  if (d_phi > PI) d_phi = 0; /* V_scatt on 4*PI */

  if (d_phi && order != 1) {
    printf("Isotropic_Sqw: %s: Focusing can only apply for single\n"
           "               scattering. Setting to 4*PI mode.\n",
           NAME_CURRENT_COMP);
    d_phi = 0;
  }

  VarSqw.w_max     = (wmax > 0 ? wmax : 0);
  VarSqw.w_min     = (wmin > 0 ? wmin : FLT_MAX);
  VarSqw.q_max     = (qmax > 0 ? qmax : 0);
  VarSqw.q_min     = (qmin > 0 ? qmin : FLT_MAX);

/* end INITIALIZE */
%}

/*****************************************************************************/
/*****************************************************************************/
TRACE
%{

int    intersect;     /* flag to continue/stop */
double t0,  t1,  t2,  t3; /* times for intersections */
double dt0, dt1, dt2, dt; /* time intervals */
double k=0, kf, kf1, kf2;
double v=0, vf;
double d_path;        /* total path length for straight trajectory */
double my_a;        /* absorption cross-section scaled to velocity (2200) */
double ws, wi;        /* probability for scattering/absorption and for */
                      /* interaction along d_path */
double tmp;           /* temporary var */
double omega=0;       /* energy transfert */
double q=0;           /* wavevector transfert */
long   index_w;       /* energy index for table look-up DOS */
long   index_q;       /* Q index for table look-up P(Q|w) */
double c1, R,theta, costheta, sintheta; /* for the choice of kf direction */
double u1x,u1y,u1z;
double u2x,u2y,u2z;
double u0x,u0y,u0z;
int    index_counter;
int    flag=0;
int    flag_concentric=0;
double my_t=0;
double p_mult=1;
double alpha, alpha0;
struct Sqw_Data_struct Data_sqw;


/* Store Initial neutron state */

VarSqw.ki_x = V2K*vx;
VarSqw.ki_y = V2K*vy;
VarSqw.ki_z = V2K*vz;
VarSqw.ti   = t;
VarSqw.vi   = 0;
VarSqw.ki   = 0;
VarSqw.type = '\0';

do { /* Main interaction loop. Ends with intersect=0 */

  /* Intersection neutron trajectory / sample (sample surface) */
  if (VarSqw.shape==1)
    intersect=box_intersect     (&t0,&t3, x,y,z,vx,vy,vz, xwidth,yheight,zthick);
  else if (VarSqw.shape==0)
    intersect=cylinder_intersect(&t0,&t3, x,y,z,vx,vy,vz, radius_o,yheight);
  else
    intersect=sphere_intersect(&t0,&t3, x,y,z,vx,vy,vz, radius_o);

  /* Computing the intermediate times */
  if (intersect) {
    tmp = 0;
    if (radius_i>0) {
      if ((VarSqw.shape==0 && cylinder_intersect(&t1,&t2, x,y,z,vx,vy,vz,
                                       radius_i,yheight))
      || (VarSqw.shape==2 && sphere_intersect(&t1,&t2, x,y,z,vx,vy,vz,
                                       radius_i)))
      tmp = 1; /* we also take into account the hollow part */
    } else if (thickness>0 && VarSqw.shape==1 && box_intersect(&t1,&t2, x,y,z,vx,vy,vz,
             xwidth-2*thickness, yheight-2*thickness, zthick-2*thickness))
      tmp = 1;
    if (!tmp) t1 = t2 = t3; /* no empty space inside */
  } else break; /* neutron does not hit sample: transmitted  */

  tmp = 0;
  if (intersect) { /* the neutron hits the sample */

    if (!v) {
      v  = vx*vx+vy*vy+vz*vz;
      v = sqrt(v);
    }
    k  = V2K*v;

    if (!VarSqw.vi) VarSqw.vi = v;
    if (!VarSqw.ki) VarSqw.ki = k;

    if (t0 > 0) {  /* we are before the sample */
      PROP_DT(t0); /* propagates neutron to the entry of the sample */
    } else if (t1 > 0 && t1 > t0) { /* we are inside first part of the sample */
      /* no propagation, stay inside */
    } else if (t2 > 0 && t2 > t1) { /* we are in the hole */
      PROP_DT(t2); /* propagate to inner surface of 2nd part of sample */
    } else if (t3 > 0 && t3 > t2) { /* we are in the 2nd part of sample */
      /* no propagation, stay inside */
    }

    dt0=t1-(t0 > 0 ? t0 : 0); /* Time in first part of hollow/cylinder/box */
    dt1=t2-(t1 > 0 ? t1 : 0); /* Time in hole */
    dt2=t3-(t2 > 0 ? t2 : 0); /* Time in 2nd part of hollow cylinder */

    if (dt0 < 0) dt0 = 0;
    if (dt1 < 0) dt1 = 0;
    if (dt2 < 0) dt2 = 0;

    /* initialize concentric mode */
    if (concentric && !flag_concentric && t0 >= 0
     && VarSqw.shape==0 && radius_i > 0) {
      flag_concentric=1;
    }

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

    if (!dt0 && !dt2) {
      intersect = 0; /* the sample was passed entirely */
      break;
    }

    VarSqw.neutron_enter++;

    if (!v) {
      printf("Isotropic_Sqw: %s: ERROR: Null velocity\n",NAME_CURRENT_COMP);
      VarSqw.neutron_removed++;
      ABSORB; /* should never occur */
    }

    /* check for scattering event */
    d_path = v*( dt0 +dt2 );      /* total path lenght in sample */
    my_a   = VarSqw.my_a_v / v; /* absorption 'mu' */

    my_t = my_a + VarSqw.my_s;  /* total scattering Xsect (tmp var) */
    if (!my_t) {
      printf("Isotropic_Sqw: %s: ERROR: Null Scattering cross section\n",NAME_CURRENT_COMP);
      VarSqw.neutron_removed++;
      ABSORB; /* should never occur */
    }

    /* Proba of scattering vs absorption */
    ws = VarSqw.my_s/my_t;  /* (inc+coh)/(inc+coh+abs) */
    /* Proba of interacing along length d_path */
    wi = 1 - exp(-my_t*d_path);

    flag = 0; /* flag used for propagation to exit point before ending */

    /* are we next to the exit ? probably no scattering */
    if (VarSqw.my_s*d_path <= 4e-7) {
      flag = 1;
      p_mult *= 1 - wi; /* No interaction before the exit */
    }

    /* MC choice: Interaction or transmission ? */
    if (!flag && wi > 0 && (wi == 1 || (tmp=rand01()) < wi)) { /* Interaction neutron/sample */
      p_mult *= ws; /* Update weight ; account for absorption and scatter */
    } else {
      flag = 1; /* Transmission : no interaction neutron/sample */
      if (!VarSqw.type) VarSqw.type = 't';
    }

    if (flag) { /* propagate to exit of sample and finish */
      intersect = 0;
      break;
    }
  } /* end if intersect the neutron hits the sample */
  else break;

  if (intersect) { /* Multiple scattering event */

    /* Decaying exponential distribution of the path length before scattering */
    /* Select a point at which to scatter the neutron, taking
         secondary extinction into account. */
    if (my_t*d_path < 1e-6)
    /* For very weak scattering, use simple uniform sampling of scattering
       point to avoid rounding errors. */
      dt = rand0max(d_path); /* length */
    else
      dt = -log(1 - rand0max((1 - exp(-my_t*d_path)))) / my_t; /* length */
    dt /= v; /* Time from present position to scattering point */

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

    /* Neutron propagation to the scattering point */
    PROP_DT(dt);

    /* choice between coherent/incoherent scattering */
    if (!tmp) tmp = rand01();
    else { /* re-use preceeding random shot */
      tmp /= wi;
    }
    tmp *= (VarSqw.s_coh+VarSqw.s_inc);

    flag=1;
    if (VarSqw.s_inc && tmp < VarSqw.s_inc) {
      /* CASE 1: incoherent case */
      if (!VarSqw.Data_inc.intensity) {
        /* CASE 1a: no incoherent Sqw from file, use isotropic V-like */
        if (d_phi && order == 1) {
          randvec_target_rect_angular(&u1x, &u1y, &u1z, &tmp,
              vx, vy, vz, 2*PI, d_phi, ROT_A_CURRENT_COMP);
          p_mult *= tmp/4*PI; /* weighted by focused range to total range */
        } else
          randvec_target_circle(&u1x, &u1y, &u1z, NULL, vx, vy, vz, 0);

        vx = u1x; vy = u1y; vz = u1z;
        vf = v; kf = k;
        flag = 0;
        if (!VarSqw.type) VarSqw.type = 'v';
        SCATTER;
      } else {
        /* CASE 1b: incoherent Sqw from file */
        if (VarSqw.Data_inc.intensity)
          Data_sqw = VarSqw.Data_inc;
        else flag=0;
        if (!VarSqw.type) VarSqw.type = 'i';
      }
    } else if (VarSqw.s_coh && tmp > VarSqw.s_inc) {
      if (VarSqw.Data_coh.intensity) {
        /* CASE2: coherent case */
        Data_sqw = VarSqw.Data_coh;
      } else flag=0;
      if (!VarSqw.type) VarSqw.type = 'c';
    }

    if (flag) { /* not true when V-like scattering */

      /* give us a limited number of tries for scattering */
      for (index_counter=VarSqw.maxloop; index_counter > 0 ; index_counter--) {

        /* MC choice: energy transfer in the DOS */
        index_w = Sqw_search_DOS(Data_sqw, rand01());
        if (&(Data_sqw.DOS[index_w]) != NULL)
          omega = Data_sqw.DOS[index_w].omega;
        else omega=0;

        /* MC choice: detailed balance: choice between w and -w */
        if (omega && VarSqw.Temperature
         && rand01() > 1/(1+exp(-omega/(VarSqw.Temperature*VarSqw.T2E))) ) {
          omega = -omega; /* anti Stokes */
        }                 /* else Stokes */

        /* MC choice: momentum transfer Q in P(Q|w) */
        index_q = Sqw_search_Q_proba_per_w(Data_sqw, rand01(), index_w);
        if (index_q == -1) {
          if (VarSqw.verbose_output >= 3)
            printf("Isotropic_Sqw: %s: Warning: No suitable P(q|w) transfert\n"
               "               for w=%g index_w=%i\n",
            NAME_CURRENT_COMP, omega, index_w);
          VarSqw.neutron_removed++;
          continue; /* no Q value for this w choice */
        }
        if (&(Data_sqw.Q_proba_per_w[index_w]) != NULL)
          q = Data_sqw.Q_proba_per_w[index_w][index_q].Q;
        else {
          if (VarSqw.verbose_output >= 3)
            printf("Isotropic_Sqw: %s: Warning: No suitable q transfert\n"
               "               for w=%g\n",
            NAME_CURRENT_COMP, omega);
          VarSqw.neutron_removed++;
          ABSORB;
        }

        /* Search for length of final wave vector kf */
        /* kf is such that : hbar*w = hbar*hbar/2/m*(k*k - kf*kf) */
        /* acceptable values for kf are kf1 and kf2 */
        if (!Sqw_search_root(1, 0, -k*k + VarSqw.sqSE2K*omega, &kf1, &kf2)) {
          if (VarSqw.verbose_output >= 3)
            printf("Isotropic_Sqw: %s: Warning: imagimary root for\n"
                 "               w=%g q=%g (triangle canm not close)\n",
            NAME_CURRENT_COMP, omega, q);
          VarSqw.neutron_removed++;
          continue; /* all roots are imaginary */
        }

        /* MC choice: between kf1 and kf2 */
        kf = fabs(rand01() < 0.5 ? kf1 : kf2);
        vf = K2V*kf;

        /* Search of the direction of kf such that : q = ki - kf */
        /* ks = c1*ki/norm(ki)^2 + u0 */
        /* u0 is in the plane orthogonal to ki */

        c1= (k*k+kf*kf-q*q)/(2*kf); /* unit is [k] */

        /* u1 and u2 : basis of the plane orthogonal to ki */
        u1x=0;           u1y=-vz;    u1z= vy; /* in vertical plane */
        u2x=vy*vy+vz*vz; u2y=-vx*vy; u2z=-vx*vz;

        /* circle radius = intersection with the plane defined by u1 and u2 */
        R = 1-c1*c1/(k*k); /* unit is [1] */
        if (R >= 0) R = sqrt(R);
        else {
          continue; /* No solution for R */
        }
        NORM(u1x, u1y, u1z);
        NORM(u2x, u2y, u2z);

        p_mult *= kf/k/4/PI;

        break; /* exit for loop on success */
      } /* end for index_counter */

      if (!index_counter) { /* for loop ended: failure for scattering */
        intersect=0; /* Could not scatter: finish multiple scattering loop */
        if (VarSqw.verbose_output >= 2)
          printf("Isotropic_Sqw: %s: Warning: No scattering [q,w] conditions\n"
               "               last try (%i): type=%c w=%g q=%g R=%g\n",
          NAME_CURRENT_COMP, VarSqw.maxloop, (VarSqw.type ? VarSqw.type : '-'), omega, q, R);
        VarSqw.neutron_removed++;
        if (order && SCATTERED != order) ABSORB;
        break;       /* finish multiple scattering loop */
      }

      /* Choose point on Debye-Scherrer cone */
      if (order == 1 && d_phi)
      { /* relate height of detector to the height on DS cone */
        tmp = sin(d_phi*DEG2RAD/2)/sin(2*theta);
/* If full Debye-Scherrer cone is within d_phi, don't focus */
        if (tmp < -1 || tmp > 1) d_phi = 0;
/* Otherwise, determine alpha to rotate from scattering plane
     into d_phi focusing area*/
        else alpha = 2*asin(tmp);
      } /* else d_phi=0 from INIT */
      if (d_phi) { /* Focusing */
        alpha = fabs(alpha);
        /* Trick to get scattering for pos/neg theta's */
        alpha0= 2*rand01()*alpha;
        if (alpha0 > alpha) {
          alpha0=PI+(alpha0-1.5*alpha);
        } else {
          alpha0=alpha0-0.5*alpha;
        }
        p_mult *= alpha/PI;
      }
      else
        alpha0 = PI*randpm1();

      costheta = cos(alpha0); sintheta = sin(alpha0);

      /* u0 = component of ks which is in the plane orthogonal to ki */
      u0x=R*(costheta*u1x + sintheta*u2x);
      u0y=R*(costheta*u1y + sintheta*u2y);
      u0z=R*(costheta*u1z + sintheta*u2z);

      tmp = v*k;

      u0x=c1*vx/tmp + u0x; /* normalized output velocity */
      u0y=c1*vy/tmp + u0y;
      u0z=c1*vz/tmp + u0z;

      /* Final velocity */
      vx = vf*u0x;
      vy = vf*u0y;
      vz = vf*u0z;

      v = vf; k = kf;
      /* nornalize to \int S or norm_sqw in selected [q,w] range */
      p_mult *= (norm_sqw > 0 ? norm_sqw : Data_sqw.intensity);
      if (p_mult > 1) VarSqw.neutron_pmult++;

      p *= p_mult;

      SCATTER;

      if      (fabs(omega) > VarSqw.w_max) VarSqw.w_max = fabs(omega);
      else if (fabs(omega) < VarSqw.w_min) VarSqw.w_min = fabs(omega);
      if      (fabs(q) > VarSqw.q_max) VarSqw.q_max = fabs(q);
      else if (fabs(q) < VarSqw.q_min) VarSqw.q_min = fabs(q);

      if (VarSqw.neutron_exit > VarSqw.maxloop
       && VarSqw.neutron_removed && VarSqw.sqw_auto) {
        if (VarSqw.verbose_output)
          printf("Isotropic_Sqw: %s:  %li events removed out of %li.\n"
                 "INFO:          Optimizing q=[%g:%g] w=[%g:%g]\n",
          NAME_CURRENT_COMP, VarSqw.neutron_removed, VarSqw.neutron_exit,
          VarSqw.q_min,   VarSqw.q_max,   VarSqw.w_min,   VarSqw.w_max);
        VarSqw.Data_inc.q_min = VarSqw.Data_coh.q_min = VarSqw.q_min;
        VarSqw.Data_inc.q_max = VarSqw.Data_coh.q_max = VarSqw.q_max;
        VarSqw.Data_inc.w_min = VarSqw.Data_coh.w_min = VarSqw.w_min;
        VarSqw.Data_inc.w_max = VarSqw.Data_coh.w_max = VarSqw.w_max;
        if (!Sqw_init(&VarSqw, Sqw_coh, Sqw_inc)) {
          printf("Isotropic_Sqw: %s: ERROR re-importing data files (Sqw_init from TRACE)\n",NAME_CURRENT_COMP);
          exit(0);
        }
        VarSqw.sqw_auto = 0;
      }

    } /* end if (flag) */

    /* test for a given multiple order */
    if (order && SCATTERED >= order) {
      intersect=0; /* reached required number of SCATTERing */
      break;       /* finish multiple scattering loop */
    }

  } /* end if (intersect) Multiple scattering event  */

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

/* Store Final neutron state */
VarSqw.kf_x = V2K*vx;
VarSqw.kf_y = V2K*vy;
VarSqw.kf_z = V2K*vz;
VarSqw.tf   = t;
VarSqw.vf = v;
VarSqw.kf = k;

if (SCATTERED) {

  VarSqw.dq = sqrt((VarSqw.kf_x-VarSqw.ki_x)*(VarSqw.kf_x-VarSqw.ki_x)
                  +(VarSqw.kf_y-VarSqw.ki_y)*(VarSqw.kf_y-VarSqw.ki_y)
                  +(VarSqw.kf_z-VarSqw.ki_z)*(VarSqw.kf_z-VarSqw.ki_z));
  VarSqw.dw = VS2E*(VarSqw.vf*VarSqw.vf - VarSqw.vi*VarSqw.vi);

  VarSqw.neutron_exit++;
} else VarSqw.dq=VarSqw.dw=0;

/* end TRACE */
%}

FINALLY
%{
  int  k;

  for (k=0; k < 2; k++) {
    struct Sqw_Data_struct Data_sqw;
    char tmp[1024];
    char type_ci[64];
    char ext[20];
    char filename[256];

    strcpy(type_ci, (k == 0 ? "coherent" : "incoherent"));
    strcpy(ext, (k == 0 ? "coh" : "inc"));
    Data_sqw =  (k == 0 ? VarSqw.Data_coh : VarSqw.Data_inc);

    if (save_sqw && Data_sqw.intensity) {
      double DOS[Data_sqw.nb_w];
      double SQW[Data_sqw.nb_q][Data_sqw.nb_w];
      long   i,j;

      if (VarSqw.verbose_output >= 1)
            printf("Isotropic_Sqw: %s: Writting %s S(q), DOS and S(q,w) files as monitors.\n", NAME_CURRENT_COMP, type_ci);

      for (j=0; j < (Data_sqw.nb_w) ; j++) {
        DOS[j] = Data_sqw.DOS[j].I;
        for (i=0; i < (Data_sqw.nb_q) ; i++)
          SQW[i][j] = Data_sqw.Q_proba_per_w[j][i].proba;
      }
      sprintf(tmp, "S(q) from %s (%s)", Data_sqw.filename, type_ci);
      sprintf(filename, "%s_%s.sq", NAME_CURRENT_COMP, ext);
      if (Data_sqw.nb_q > 1)
        DETECTOR_OUT_1D(tmp,"Wavevector [Angs-1]","S(q)","q",
        Data_sqw.q_min, Data_sqw.q_max,
        Data_sqw.nb_q,NULL,Data_sqw.Sq,NULL,
        filename);
      sprintf(tmp, "Normalized DOS from %s (%s)", Data_sqw.filename, type_ci);
      sprintf(filename, "%s_%s.dos", NAME_CURRENT_COMP, ext);
      if (Data_sqw.nb_w > 1)
        DETECTOR_OUT_1D(tmp,"Energy [meV]","g(w)","w",
        Data_sqw.w_min, Data_sqw.w_max,
        Data_sqw.nb_w,NULL,DOS,NULL,
        filename);
      sprintf(tmp, "S(q,w) from %s (%s), Norm=%g", Data_sqw.filename, type_ci,
        (norm_sqw > 0 ? norm_sqw : Data_sqw.intensity));
      sprintf(filename, "%s_%s.sqw", NAME_CURRENT_COMP, ext);
      if (Data_sqw.nb_q > 1 && Data_sqw.nb_w > 1)
        DETECTOR_OUT_2D(tmp,"Wavevector [Angs-1]","Energy [meV]",
        Data_sqw.q_min, Data_sqw.q_max, Data_sqw.w_min, Data_sqw.w_max,
        Data_sqw.nb_q,Data_sqw.nb_w,NULL,&SQW[0][0],NULL,
        filename);
    }
    /* test if Sqw data is larger than used [q,w] range */
    if (VarSqw.neutron_removed && Data_sqw.nb_q*Data_sqw.nb_w > 1)
    if ((VarSqw.w_max < Data_sqw.w_max || VarSqw.w_min > Data_sqw.w_min)
     || (VarSqw.q_max < Data_sqw.q_max || VarSqw.q_min > Data_sqw.q_min))
          printf("Isotropic_Sqw: %s: Data from %s file %s in the range\n"
                 "                  q=[%g:%g] w=[%g:%g]\n"
                 "                  is LARGER than accessible beam range.\n"
                 "               This may result in neutron events removed.\n"
                 "WARNING        Set qmin=%g, qmax=%g, wmin=%g, wmax=%g\n"
                 "               to Focus simulation optimally.\n",
            NAME_CURRENT_COMP, type_ci, Data_sqw.filename,
            Data_sqw.q_min, Data_sqw.q_max, Data_sqw.w_min, Data_sqw.w_max,
            VarSqw.q_min,   VarSqw.q_max,   VarSqw.w_min,   VarSqw.w_max);

    if (Data_sqw.DOS)           free(Data_sqw.DOS);
    if (Data_sqw.Q_proba_per_w) free(Data_sqw.Q_proba_per_w);
    if (Data_sqw.DOS_lookup)    free(Data_sqw.DOS_lookup);
    if (Data_sqw.QW_lookup)     free(Data_sqw.QW_lookup);
    if (Data_sqw.Sq)            free(Data_sqw.Sq);
  }

  if (VarSqw.neutron_removed)
    printf("Isotropic_Sqw: %s: %li neutron events (out of %li) that should have\n"
           "               scattered were removed because scattering conditions\n"
           "WARNING        could not be satisfied after %i tries.\n",
          NAME_CURRENT_COMP, VarSqw.neutron_removed, VarSqw.neutron_exit, VarSqw.maxloop);
  if (VarSqw.neutron_pmult)
    printf("Isotropic_Sqw: %s: %li neutron events (out of %li) reached\n"
           "               unrealistic weight. The S(q,w) norm is probably too high\n"
           "WARNING        could not be satisfied after %i tries.\n",
          NAME_CURRENT_COMP, VarSqw.neutron_pmult, VarSqw.neutron_exit, VarSqw.maxloop);


/* end FINALLY */
%}
/*****************************************************************************/
/*****************************************************************************/

MCDISPLAY
%{
  magnify("xyz");
  if(VarSqw.shape==1)
  {
    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);

    if (thickness>0) {
      xmin = -0.5*xwidth+thickness;
      xmax = -xmin;
      ymin = -0.5*yheight+thickness;
      ymax = -ymin;
      zmin = -0.5*zthick+thickness;
      zmax = -zmin;
      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(VarSqw.shape==0)
  {
    circle("xz", 0,  yheight/2.0, 0, radius_i);
    circle("xz", 0,  yheight/2.0, 0, radius_o);
    circle("xz", 0, -yheight/2.0, 0, radius_i);
    circle("xz", 0, -yheight/2.0, 0, radius_o);
    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);
    line(-radius_o, -yheight/2.0, 0, -radius_o, +yheight/2.0, 0);
    line(+radius_o, -yheight/2.0, 0, +radius_o, +yheight/2.0, 0);
    line(0, -yheight/2.0, -radius_o, 0, +yheight/2.0, -radius_o);
    line(0, -yheight/2.0, +radius_o, 0, +yheight/2.0, +radius_o);
  } else {
    circle("xy",0,0,0,radius_i);
    circle("xz",0,0,0,radius_i);
    circle("yz",0,0,0,radius_i);
    circle("xy",0,0,0,radius_o);
    circle("xz",0,0,0,radius_o);
    circle("yz",0,0,0,radius_o);
  }
/* end MCDISPLAY */
%}

/*****************************************************************************/
/*****************************************************************************/

END