/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Instrument: ILL_IN4
*
* %Identification
* Written by: <a href="mailto:farhi@ill.fr">Emmanuel Farhi</a>
* Date: March 5th 2015.
* Origin: <a href="http://www.ill.fr">ILL (France)</a>
* %INSTRUMENT_SITE: ILL
*
* The IN4 thermal Time-of-Flight spectrometer at the ILL (H12 tube).
*
* %Description
* IN4C is a high-flux time-of-flight spectrometer used for the study of excitations
* in condensed matter. It works in the thermal neutron energy range (10-100 meV).
*
* Primary spectrometer
*
* The main components of the beam conditioning part are the two background
* choppers, the double curvature mono-chromator with four faces and the Fermi
* chopper. The background choppers are rapidly pulsating beam shutters which act
* as a low-pass filter. Thus they eliminate from the beam most of the fast
* neutrons and gamma rays that would give background noise in the spectra. The
* modular shielding encloses the background choppers in separate compartments in
* order to cut off these undesired neutrons as early as possible. A suitable
* energy is selected from the thermal neutron spectrum with the crystal
* monochromator. The monochromator, an assembly of 55 crystal pieces,
* concentrates the divergent incident beam onto a small area at the sample
* position. The full use of the available solid angle gives a high incident
* flux. The vertical curvature is fixed, and the horizontal
* variable curvature of the monochromator is essential in controlling
* the time and space focussing conditions for optimal performance (see H. Mutka,
* Nucl. Instr. and Meth. A 338 (1994) 144). The Fermi chopper rotates at speeds
* of up to 40000 rpm. It transmits short neutron pulses (10 ... 50 µs) to the
* sample. The time-of-flight of neutrons between the chopper and the sample (1
* ... 5 ms) can be measured by using precise electronic circuitry.
* A sapphire (Al2O3) filter can be inserted in the beam to remove the fast neutrons
* background.
*
* Monochromators:
* PG       002 DM=3.355 AA (Highly Oriented Pyrolythic Graphite)
* PG       004 DM=1.677 AA (used for lambda=1.1)
* PG       006 DM=1.118 AA
* Cu       220 DM=1.278 AA
* Cu       111 DM=2.095 AA
* Take-off:       39-65 deg
* flux at sample: 5e5 n/s/cm2 (at 1.1 Angs)
*
* Secondary spectrometer
*
* The sample environment is designed to accommodate standard
* cryostats and furnaces. A radial collimator around the sample position is used
* to cut the scattering from the sample environment. The secondary flight-path
* is in vacuum to avoid parasitic scattering of the transmitted neutrons. The
* detector bank covers scattering angles of up to 120°. In addition to the 3He
* detector tubes (length 300 mm, width 30 mm, elliptical section, pressure 6
* bar) a 3He filled multidetector (eight sectors with 12 radial cells each;
* outer diameter Phi 60 cm) will allow us to observe forward scattering. The
* time-of-flight spectra measured at various angles are further treated in order
* to obtain the scattering function S(Q,w) using e.g. LAMP.
*
* In this model, the sample is a cylindrical liquid/powder/glass scatterer
* surrounded by a container and an Al cryostat.
*
* %Example: lambda=1.2 DM=1.677 Detector: sample_flux_I=4.43306e+06
*
* %Parameters
* lambda: [Angs]         wavelength
* dlambda: [Angs]        wavelength HALF spread at source
* DM: [Angs]             monochromator d-spacing
* ETAM: [arcmin]         monochromator mosaic FWHM
* RMH: [m]               Monochromator horizontal curvature. Use -1 for auto.
* ratio: [1]             Disk Chopper ratio (nu=nu_FC/ratio)
* dE: [meV]              Inelastic energy for time focusing
* Sapphire_present: [1]  Flag, when 1 the Al2O3 filter is active
* sample_coh: [str]      Sample coherent dynamic structure factor (.sqw) or NULL
* sample_inc: [str]      Sample incoherent dynamic structure factor (.sqw) or NULL
* order: [1]             The number of multiple orders at the monochromator
*
* %Link
* H. Mutka, Nucl. Instr. and Meth. A 338 (1994) 144
* %Link
* http://www.ill.eu/fr/instruments-support/instruments-groups/instruments/in4c
* %End
*******************************************************************************/
DEFINE INSTRUMENT ILL_IN4(lambda=2.2, dlambda=0.1, DM=3.355, ETAM=35, RMH=-1, ratio=4, dE=0, Sapphire_present=1, string sample_coh="Dirac2D.sqw", string sample_inc="NULL", int order=1)

DECLARE %{
  double A1, RMH, RMV=-1.9;
  double LFS       = 0.675; /* FC to Sample distance */
  double LSD       = 2;     /* Sample to Detector distance */
  double LVS       = 0.583; /* BC1 to VS distance */
  double sample_width          = 0.05;
  double sample_thickness      = 1e-3;
  double sample_height         = 0.05;
  double environment_radius    = 0.05;
  double environment_thickness = 1e-3;
  #pragma acc declare create( environment_thickness )
  double container_thickness   = 5e-4;
  #pragma acc declare create( container_thickness )
  char   environment[]         = "Al.laz";
  char   container[]           = "V.laz";

  /* the following variables are computed in the IN4 configuration */
  double d0, d2; /* IN4c distances used in NoMad/IN4 doc (Ross/Rols) */
  double phase1F   = 0;
  double nu        = 0;
  double LMS       = 0; /* [m]    Mono-sample distance (aka d1) */
  double LRM       = 0; /* [m]    Distance from source to monochromator. */
  double bctr      = 0; /* [m]    background chopper BC2 translation from BC1 */
  double phase12   = 0; /* [deg]  Chopper phase BC2 wrt BC1 */
  char   mon_sqw[1024];
  char   mon_qe[1024];

  double Lmin, Lmax;

  double t1F=0;
  double t12=0;

  double Ei=0;

  /* Write_Sqw: writes an Sqw file for Isotropic_Sqw as a Dirac 2D array
     INPUT:
       dq:        [1]   momentum binning between adjacent Dirac peaks. 0 to get bins/10.
       dw         [1]   energy   binning between adjacent Dirac peaks. 0 to get bins/10.
       max_q:     [Angs] maximum momentum transfer to generate [Angs]
       max_w:     [meV] maximum energy transfer to generate [meV]
       bins:      [1]   S(q,w) size is bins x bins.
       filename:  [str] output filename
     OUTPUT:
       a file containing S(q,w)
       the return value is the number of Durac peaks generated

     Simple usage: Write_Sqw_Dirac(0, 0, 20, 500, 1000, "Dirac2D.sqw") */

  long Write_Sqw_Dirac(long dq, long dw,
    double max_q, double max_w, long bins, char *filename)
  {

    double index_q, index_w;
    double min_q=0,    min_w=0;
    long   bin_q=0,    bin_w=0;
    long   count=0;

    FILE *fid = NULL;

    fid = fopen(filename, "w+");
    if (!fid) return(0);

    /* check binning */
    if (bins <= 0) bins=1000; /* 1000x1000 makes a 8 Mb array */
    if(dq <= 0)    dq=bins/10;
    if(dw <= 0)    dw=bins/10;
    bin_q=bins; bin_w=bins;

    /* write header */
    fprintf(fid,
      "# Sqw data file for Isotropic_Sqw\n"
      "# model S(q,w) as a set of Dirac peaks, to obtain the 2D resolution function\n"
      "# on a (q,w) grid. (q,w) grid=%g Angs, %g meV.\n"
      "# filename: %s\n"
      "#\n"
      "# Physical parameters:\n"
      "# sigma_coh 1    coherent scattering cross section in [barn]\n"
      "# sigma_inc 0    incoherent scattering cross section in [barn]\n"
      "# sigma_abs 0    absorption scattering cross section in [barn]\n"
      "# density   1    in [g/cm^3.5]\n"
      "# weight    1    in [g/mol]\n"
      "# nb_atoms  1    in [atoms/unit cell]\n"
      "#\n", (max_q-min_q)/bin_q, (max_w-min_w)/bin_w, filename);

    /* write q axis */
    fprintf(fid,
      "# WAVEVECTOR vector of m=%i values %g:%g in [Angstroem-1]: q\n",
      bin_q, min_q, max_q);
    for (index_q=0; index_q < bin_q; index_q++) {
      double q = min_q+index_q*(max_q-min_q)/bin_q;
      fprintf(fid, "%g ", q);
    }
    fprintf(fid, "\n");

    /* write w axis */
    fprintf(fid,
      "# ENERGY vector of n=%i values %g:%g in [meV]: w\n",
      bin_w, min_w, max_w);
    for (index_w=0; index_w < bin_w; index_w++) {
      double w = min_w+index_w*(max_w-min_w)/bin_w;
      fprintf(fid, "%g ", w);
    }
    fprintf(fid, "\n");

    /* write Sqw matrix as zero except when on grid */
    fprintf(fid,
      "# matrix of S(q,w) values (m=%i rows x n=%i columns), one row per q value: sqw\n",
      bin_q, bin_w);
    for (index_q=0; index_q < bin_q; index_q++) {
      for (index_w=0; index_w < bin_w; index_w++) {
        double sqw = 0;
        if (fmod(index_q, dq) == 0 && fmod(index_w, dw) == 0) {
          sqw=1; count++; }
        fprintf(fid, "%g ", sqw);
      }
      fprintf(fid, "\n");
    }
    fprintf(fid, "\n");
    fprintf(fid, "# end of Sqw file %s\n", filename);

    fclose(fid);
    return(count); /* size of S(q,w) generated */
  } /* Write_Sqw_Dirac */
%}

USERVARS %{
  char   flag_sample_choice;
  char flag_source_order;
  char order_extend;
  char flag_sample;
  char flag_env;
  double vix;
  double viy;
  double viz;
%}

INITIALIZE %{
  double Vi, Ki;
  double thetaB;
  double dSx=5.556;   /* BC1 to sample along beam tube axis */
  double dSy=1.3;     /* lateral position of sample */
  double L1F, L1M, LMF;

  if (DM <= 0) {
    if      (lambda < 1.0) DM=1.118;
    else if (lambda < 1.8) DM=1.677;
    else                   DM=3.355;
  }

  thetaB = -asin(lambda/DM/2);
  A1     = thetaB*RAD2DEG;

  Ki = 2*PI/lambda;
  Vi = K2V*fabs(Ki);
  Ei = VS2E*Vi*Vi;

  /* IN4c configuration */


  /* compute distances for IN4 */
  d2 = fabs(dSy/tan(2*thetaB));
  LMS= fabs(dSy/sin(2*thetaB)); /* Monok to Sample = d1 */
  d0 = dSx - d2 - LVS;    /* VS to Monok */

  /* set distances for IN4c */
  LRM = 6.1 + LVS + d0;
  L1M= LRM - 6.1;      /* BC1 to Monok */
  LMF= LMS-LFS;        /* Monok to Fermi */
  L1F= L1M + LMF;      /* BC1 to Fermi */

  /* FC rotation frequency (Master) */
  nu  = K2V/( fabs(DM*cos(thetaB)) * (LFS+LSD*pow( 1-dE/Ei, -1.5)) *(1-LMS/d0) );

  phase12 = 22.5;
  /* compute position of BC2 and phases */
  t12 = phase12/360/(nu/ratio); /* time delay [s] */

  bctr  = t12*Vi;
  if (bctr > 2.965) { /* chopper BC2 at its maximum position, from BC1 */
    bctr = 2.965;
  }

  /* compute back the phases (in case bctr has changed) */
  phase12 = -bctr*360*(nu/ratio)/Vi;
  /* distance BC1-FC */
  t1F     = ((LMS-LFS)+(LRM-6.1))/Vi;
  phase1F = -t1F*360*nu;

  if (RMH < 0) {
    double L = 1/(1/d0+1/LMS); /* Monok optical focusing distance */
    RMH= 2*L/sin(DEG2RAD*A1);  /* RH = 2*L/sin(DEG2RAD*A1); */
  }

  if (dlambda <= 0) dlambda = lambda*.95;

  Lmin = lambda-dlambda;
  Lmax = lambda+dlambda;
  if (Lmin < 0)   Lmin = 0.1;
  if (Lmax > 3.5) Lmax = 3.5;

  MPI_MASTER(
  /* print some information when starting */
  printf("%s: Thermal ToF spectrometer\n", NAME_INSTRUMENT);
  printf("  Divergence at the lead shutter: dX    =%g [deg]\n", atan2(0.2,5.2)*RAD2DEG);
  printf("  Take-off at monochromator:      A1    =%g [deg] ; DM=%g [Angs]\n", A1, DM);
  printf("  Incident energy:                Ei    =%g [meV] ; Ki=%g [Angs-1]\n", Ei, Ki);
  printf("  Incident velocity:              Vi    =%g [m/s]\n", Vi);
  printf("  Source-Mono distance:           LRM   =%g [m]\n", LRM);
  printf("  Virtual Source-Mono distance:   d0    =%g [m]\n", d0);
  printf("  Mono-Sample distance:           LMS   =%g [m] (d1)\n", LMS);
  printf("  Curvature at monochromator:     RMH   =%g [m] ; RMV=%g [m]\n", RMH, RMV);
  printf("  Fermi Chopper Frequency:        nu    =%g [Hz] ; rpm=%g [rpm]\n", nu, nu*60);
  printf("  BC2 phase wrt BC1:              PhiBC2=%g [deg] L12=%g [m] delay t12=%g [s] (BCTR)\n",
    phase12, bctr, t12);
  printf("  FC phase wrt BC1:               PhiFC =%g [deg] L1F=%g [m] delay t1F=%g [s]\n",
    phase1F, L1F, t1F);

  /* print a visual representation of distances */
  printf("Distances: in [m]\n");
  printf("[H12] %g [BC1] %g [VS] %g [BC2] %g [PG] %g [FC] %g [Spl] %g [Det]\n",
    6.1, LVS, bctr - LVS, LRM - (6.1 + bctr), LMS - LFS, LFS, LSD);

  /* generate Sqw Dirac array */
  printf("\nGenerate %li Dirac peaks on (q,w) grid.\n",
    Write_Sqw_Dirac(0, 0, 3*Ki, 5*Ei, 1000, "Dirac2D.sqw"));
  );

  if (nu < 0 || nu > 700) exit(printf("ERROR: Invalid: Fermi master frequency nu=%g [Hz]. Change DM ?\n", nu));
  if (fabs(A1) < 5)      exit(printf("ERROR: Invalid: mono take-off angle A1=%g [deg]. Change DM ?\n", A1));
  if (bctr < 0)           exit(printf("Invalid: BC1-BC2 distance bctr=%g [m]. Change DM ?\n", bctr));
  if (fmod(phase12, 45.0) < 22 )
    printf("%s: WARNING: The choppers are NOT in closed position phase12=%g [deg].\n",
      NAME_INSTRUMENT, phase12);

  sprintf(mon_sqw, "user1 limits=[0 %g], user2 limits=[%g %g]", 3*Ki, -Ei, 4*Ei);
  sprintf(mon_qe,  "banana, angle limits=[-150 150], energy limits=[0 %g]", 4*Ei);

  #pragma acc update device( environment_thickness )

  #ifdef USE_MPI
  /* In case of MPI mode, add an MPI_Barrier to ensure no node runs on an uncompleted/non-existentt Dirac2D.sqw */
  MPI_Barrier(MPI_COMM_WORLD);
  #endif
%}


TRACE

COMPONENT Origin = Progress_bar()
  AT (0, 0, 0) ABSOLUTE
  EXTEND %{
    flag_source_order = floor(rand01()*order_extend*.99);
    /* tests for consistency */
    int ord=INSTRUMENT_GETPAR(order);
    if ((ord <= 0 || ord > 4)) {
      ord = 4;
    }
    order_extend = ord;
  %}

COMPONENT Thermal = Source_gen(
  radius   = 0.10/2,
  focus_xw = 0.1,
  focus_yh = 0.1,
  dist  =5.2,
  lambda0=lambda, dlambda=dlambda,
  T1=683.7,I1=0.5874e+13,T2=257.7,I2=2.5099e+13,T3=16.7 ,I3=1.0343e+12,
  verbose  = 1)
  WHEN (flag_source_order == 0)
  AT (0, 0, 0) RELATIVE Origin

COMPONENT Thermal2 = COPY(Thermal)
  (lambda0=lambda/2)
  WHEN (flag_source_order == 1)
  AT (0, 0, 0) RELATIVE Origin

COMPONENT Thermal3 = COPY(Thermal)
  (lambda0=lambda/3)
  WHEN (flag_source_order == 2)
  AT (0, 0, 0) RELATIVE Origin

COMPONENT Thermal4 = COPY(Thermal)
  (lambda0=lambda/4)
  WHEN (flag_source_order >= 3)
  AT (0, 0, 0) RELATIVE Origin

/* could we put a guide in pile ? max l=2 */

/* bouchon barillet Phi=100mm, at 5.2 m */
COMPONENT Obt1 = Monitor_nD(
  xwidth=0.1, yheight=0.1, options="disk, slit, x y", bins=100)
  AT (0, 0, 5.2) RELATIVE Thermal
  EXTEND %{
    p *= order_extend;
  %}

COMPONENT Obt1_lambda = Monitor_nD(xwidth=0.1, yheight=0.1, options="lambda limits=[.1 3.5]", bins=100)
  AT (0, 0, 0.01) RELATIVE Obt1

/* sapphire filter to remove fast neutrons
   c along beam axis, e=90 60x110 mm
 */

COMPONENT SapphireFilter = Filter_gen(xwidth=0.12, yheight=0.12,
  filename="Al2O3_sapphire.trm")
  WHEN (Sapphire_present)
  AT (0,0,0.2) RELATIVE Obt1

COMPONENT Win1 = Monitor_nD(
  xwidth=0.12, yheight=0.12, options="disk, slit, x y", bins=100)
  AT (0,0, 0.4) RELATIVE Obt1

COMPONENT Win1_lambda = Monitor_nD(xwidth=0.06, yheight=0.1, options="lambda limits=[.1 3.5]", bins=100)
  AT (0, 0, 0.01) RELATIVE Win1

/* BC1 should be as early as possible. Opening of slits: 22.5 deg i.e. 6 cm */
COMPONENT BC1  = DiskChopper(radius=0.603/2, nslit=8, isfirst=1,
  theta_0=22.5, nu=nu/ratio, yheight=0.10, phase=0)
  AT (0,0, 6.1) RELATIVE Thermal

/* the "Virtual Source" (which is just a slit) */
/* the monochromator makes an image of it onto the sample */
COMPONENT VS   = Slit(xwidth=0.08, yheight=0.20)
AT (0,0, 0.583) RELATIVE BC1

COMPONENT BC2_slit = Slit(xwidth=0.07, yheight=.12)
  AT (0,0, bctr-0.02) RELATIVE BC1

COMPONENT BC2_t = Monitor_nD(
  xwidth=0.07, yheight=0.12, options="t limits=[0.0001 0.0015]", bins=100)
  AT (0,0, bctr-0.01) RELATIVE BC1

COMPONENT BC2  = DiskChopper(radius=.643/2, nslit=8, theta_0=22.5,
  nu=nu/ratio, delay=t12, yheight=0.14)
AT (0,0, bctr) RELATIVE BC1

COMPONENT BC2_t_post = COPY(BC2_t)
  AT (0,0, 0.01) RELATIVE BC2

COMPONENT Cradle = Monitor_nD(
  options="x y", bins=50, xwidth=.25, yheight=.25, restore_neutron=1)
  AT (0,0,LRM) RELATIVE Thermal

COMPONENT Cradle_lambda = Monitor_nD(
  options="lambda limits=[.1 3.5]", bins=100, xwidth=.25, yheight=.25, restore_neutron=1)
  AT (0,0,LRM+0.01) RELATIVE Thermal

COMPONENT Cradle_t = Monitor_nD(
  options="t limits=[0.0005 0.0019]", bins=100, xwidth=.25, yheight=.25,
  restore_neutron=1)
  AT (0,0,LRM) RELATIVE Thermal

COMPONENT Mono_xy = Monitor_nD(
  options="x y", bins=50, xwidth=.22, yheight=.2, restore_neutron=1)
  AT (0,0,0) RELATIVE Cradle
  ROTATED (0, A1, 0) RELATIVE Cradle

SPLIT COMPONENT Mono = Monochromator_curved(
  width=0.22, height=0.2, NH=11, NV=5, RV=RMV, RH=RMH,
  DM=DM, mosaic=ETAM, reflect="HOPG.rfl", transmit="HOPG.trm")
  AT (0,0,0) RELATIVE Cradle
  ROTATED (0, A1, 0) RELATIVE Cradle
  EXTEND %{
    if (!SCATTERED) ABSORB; /* remove transmitted beam */
  %}

COMPONENT Mono_out = Arm()
  AT (0,0,0) RELATIVE Cradle
  ROTATED (0, 2*A1, 0) RELATIVE Cradle

COMPONENT Mono_t = COPY(Cradle_t)
  AT (0,0,0.1*(LMS-LFS)) RELATIVE Mono_out

COMPONENT FC_Pos = Monitor_nD(
  yheight=0.064, xwidth=0.03, options="t limits=[0.001 0.0023]", bins=100,
  restore_neutron=1)
AT (0,0,LMS-LFS) RELATIVE Mono_out

COMPONENT FC_Slit = Slit(yheight=0.064, xwidth=0.03)
AT (0,0,-0.03) RELATIVE FC_Pos

COMPONENT Fermi = FermiChopper(delay=t1F, radius=0.025, nu=-nu,
   yheight=0.064, xwidth=0.023, nslit=50, length=0.025)
AT (0,0,0) RELATIVE FC_Pos
EXTEND %{
    if (!SCATTERED) ABSORB;
    vix=vx,viy=vy,viz=vz;
  %}

/* sample position */
COMPONENT Sample_pos = Arm()
  AT (0,0,LMS) RELATIVE Mono_out

COMPONENT Sample_rot = Arm()
  AT (0,0,0) RELATIVE Sample_pos
  ROTATED (0,45,0) RELATIVE Sample_pos

SPLIT 100 COMPONENT sample_flux = Monitor_nD(
    xwidth = 0.06, yheight = 0.06, options = "x y",
    restore_neutron=1, bins=60)
  AT (0, 0, 0) RELATIVE Sample_pos

COMPONENT reset = Arm()
AT (0,0,0) RELATIVE PREVIOUS
EXTEND %{
  flag_sample = 0;
  flag_env = 0;
  flag_sample_choice = (rand01() > 0.5 ? 1 : 2);
  p *= 2; /* compensate for MC choice on 2 samples */
%}

COMPONENT sample_tof = Monitor_nD(
  xwidth = 0.1, yheight = 0.1, options = "x, time limits=[0.0018 0.0019]",
  bins=100, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample_pos

COMPONENT sample_lambda = Monitor_nD(
  xwidth = 0.1, yheight = 0.1, options = "lambda limits=[.1 3.5]",
  bins=300, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample_pos

COMPONENT sample_w = Monitor_nD(
  xwidth = 0.1, yheight = 0.1, options = "energy", min=0, max=2*Ei,
  bins=100, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample_pos

COMPONENT sample_qxy = Monitor_nD(
  xwidth = 0.1, yheight = 0.1, options = "kx limits=[-0.19 0.19] ky limits=[-0.16 0.16]",
  bins=100, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample_pos

/* sample environment and cell */

/* external shield */
COMPONENT Environment_in=Isotropic_Sqw(
  radius = environment_radius, yheight = 0.1, thickness=environment_thickness,
  Sqw_coh=environment, concentric=1, verbose=0, order=1, d_phi=2*RAD2DEG*atan2(1, LSD)
) WHEN (environment_thickness > 0)
AT (0, 0, 0) RELATIVE Sample_rot
EXTEND %{
  flag_env += SCATTERED;
%}

/* sample container */
COMPONENT Container_in=Isotropic_Sqw(
  xwidth  = sample_width+1e-4+container_thickness,
  zdepth  = sample_thickness+1e-4+container_thickness,
  yheight = sample_height, thickness=container_thickness,
  Sqw_coh=container, concentric=1, verbose=0, order=1, d_phi=2*RAD2DEG*atan2(1, LSD)
  ) WHEN(container_thickness > 0)
AT (0, 0, 0) RELATIVE Sample_rot
EXTEND
%{
  flag_env += SCATTERED;
%}

COMPONENT SampleS=Isotropic_Sqw(
  xwidth = sample_width, zdepth=sample_thickness, yheight = sample_height,
  Sqw_coh= sample_coh, Sqw_inc= sample_inc, p_interact=0.9,
  d_phi=2*RAD2DEG*atan2(1, LSD), order=1)
WHEN (flag_sample_choice == 1)
AT (0, 0, 0)   RELATIVE Sample_rot
EXTEND
%{
  flag_sample += SCATTERED;
%}

COMPONENT SampleV=Incoherent(xwidth = sample_width, zdepth=sample_thickness, yheight = sample_height,
  focus_ah = 2*RAD2DEG*atan2(1, LSD), focus_aw=150.0)
WHEN (flag_sample_choice == 2)
AT (0, 0, 0)   RELATIVE Sample_rot
EXTEND
%{
  flag_sample += SCATTERED;
%}

COMPONENT Container_out=COPY(Container_in)(concentric=0)
WHEN(container_thickness > 0)
AT (0, 0, 0) RELATIVE Sample_rot
EXTEND
%{
  flag_env += SCATTERED;
%}

/* external shield */
COMPONENT Environment_out=COPY(Environment_in)(concentric=0)
WHEN (environment_thickness > 0)
AT (0, 0, 0) RELATIVE Sample_rot
EXTEND %{
  flag_env += SCATTERED;
%}

COMPONENT Detector = Monitor_nD(radius=LSD, yheight=2, restore_neutron=1,
  options="theta limits=[-15 135] bins=300, y bins=100, banana")
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_sample = COPY(Detector)
  WHEN (flag_sample)
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_env = COPY(Detector)
  WHEN (flag_env)
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_Sqw = Sqw_monitor(
  nq=1000, nE=1000, vix="vix",viy="viy",viz="viz",
    qmin=0, qmax=10, Emin=-50, Emax=50, filename="Sqw_full")
  WHEN (flag_sample && flag_sample_choice == 1)
  /* "user1 limits=[0 10], user2 limits=[-50 50]" */
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_qe = Monitor_nD(
  radius=LSD, yheight=2, bins=500, restore_neutron=1,
  options=mon_qe)
  WHEN (flag_sample && flag_sample_choice == 1)
  /* "banana, angle limits=[-150 150], energy limits=[0 50]" */
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_SqwV = COPY(Detector_Sqw)(filename="Sqw_filtered")
  WHEN (flag_sample && flag_sample_choice == 2)
  AT (0,0,0) RELATIVE Sample_pos

COMPONENT Detector_qeV = COPY(Detector_qe)
  WHEN (flag_sample && flag_sample_choice == 2)
  AT (0,0,0) RELATIVE Sample_pos

/* display static shapes for viz puposes only =============================== */
COMPONENT Fuel_centre = Arm()
  AT (0.47, 0, -0.22) RELATIVE Thermal

COMPONENT H12_tube = Shape(radius=0.1/2, yheight=5)
  AT (0, 0, 2.5)     RELATIVE Thermal
  ROTATED (90, 0, 0) RELATIVE Thermal

COMPONENT Fuel = Shape(radius=.2, yheight=1.3)
  AT (0, 0, 0) RELATIVE Fuel_centre

COMPONENT D2O_vessel = Shape(radius=1.3, yheight=1)
  AT (0, 0, 0) RELATIVE Fuel_centre

COMPONENT ILL5_wall = Shape(xwidth=4, yheight=3, zdepth=.5)
  AT (0,0,27) RELATIVE Fuel_centre

COMPONENT H2O_vessel = Shape(radius=3, yheight=1)
  AT (0, 0, 0) RELATIVE Fuel_centre

COMPONENT Concrete_wall = Shape(radius=4.75, yheight=1)
  AT (0, 0, 0) RELATIVE Fuel_centre

COMPONENT Barillet = Shape(radius=0.5, yheight=0.5)
  AT (0,0, 4.7+0.5/2) RELATIVE Thermal
  ROTATED (90, 0, 0)  RELATIVE Thermal

COMPONENT Pillar = Shape(xwidth=0.5, yheight=2, zdepth=0.5)
  AT (1.5, 0, 17.5)   RELATIVE Thermal

END