/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: ILL_BRISP
*
* %Identification
* Written by: E. Farhi and N. Formissano [formisan@ill.fr]
* Date: June 4th, 2009
* Origin: ILL
* %INSTRUMENT_SITE: ILL
*
* Time of Flight Neutron Spectrometer for Small Angle Inelastic Scattering BRISP
*
* %Description
* BRISP is a new concept thermal neutron Brillouin scattering spectrometer which
* exploits the time-of-flight technique and is optimized to operate at small
* scattering angles with good energy resolution.
*
* Keywords in the design of the BRISP spectrometer were :
*
* Thermal neutron energies: allowing for investigations in systems characterized
*   by sound velocities up to 3000 m/s (three different incident energies between
*   20 and 80 meV are presently available).
* Easy small-angle access: enabling low-Q spectroscopy with thermal neutrons.
*   Elastic wavevector transfer values Qel as low as 0.03 Å -1 at 20 meV incident
*   energy can be reached. The position of the two-dimensional detector can be
*   adjusted to cover different small-angle ranges between 1° and 15°.
* Time-of-Flight technique: for an efficient data collection allowing also for
*   accurate neutron measurements as a function of external parameters such as
*   temperature, pressure and magnetic field.
* Carefull optimization of monochromator-collimators-Fermi chopper:  leading to
*   0.5 meV energy resolution and 0.02 Å-1 Q resolution in a typical
*   configuration (20 meV incident energy and 4 m sample-detector distance),
*   along with acceptable counting rates (flux at the sample 104 n s-1 cm-2).
*   For this purpose, innovatory solutions were specially developed for some of
*   the BRISP components.
*
* Main components:
*
* a Soller collimator defining the beam impinging on the monochromator, with a
*   collimation angle of 0.4°
* two focusing multi-crystal monochromators, PG and Cu(111), that allow for the
*   selection of three incident energies in the range from 20 to 80 meV.
*   Fixed/variable curvatures are adopted in/outside the Brisp vertical scattering plane.
* a disk chopper used for background reduction and selection of the desired
*   monochromator reflection through proper phasing with the Fermi chopper.
* three honeycomb converging collimators [1] to define the incident beam on the
*   sample with a collimation angle of 0.4°, and to optimize convergence at three
*   detector positions (2, 4, 6 m from the sample). A coarse resolution option
*   is also available, without honeycomb collimator.
* a Fermi chopper producing short neutron pulses which enable the time-of-flight
*   analysis.
* a high-vacuum sample chamber possibly equipped with 1.5-300 K MAXI Orange
*   cryostat (100 mm) and 300-1900 K furnace
* a ~2 m2-area position sensitive gas detector (3He) whose distance from the
*   sample can be varied between 2 and 6 m in order to access the required Q-range.
*   A huge vacuum tank hosts the detector. An elastobore – polyethylene shielding
*   surrounds the vacuum tank to reduce the environmental background.
* the long vacuum line ensures an under-vacuum neutron flight path from the
*   background chopper to the detector.
*
* Configurations:
* crystal d-spacing (Å)   lambda0 (Å)    E0(meV)
* PG(002) 3.355(nominal)  1.977(expt.)   20.9 (expt.)
* Cu(111) 2.087           1.28  (expt.)  49.9 (expt.)
* PG(004) 1.677(nominal)  0.989(expt.)   83.6 (expt.)
*
* In this model, the sample is a plate of thickness e=4 mm, surrounded by an
* Al or Nb container, inside an Al shield (phi=10 cm).
*
* %Example: coh="V.lau" Detector: Detector_I=31.55
*
* %Parameters
* DM: [Angs]        Monochromator d-spacing. Use 3.355 for PG002, 1.677 for PG004 and 2.087 for Cu111.
* coh: [str]        Sample coherent specification (use laz, lau or Sqw file, or NULL to disable). Sample is a 5x5 cm plate, e=4 mm.
* inc: [str]        Sample incoherent specification (use laz, lau or Sqw file, or NULL to scatter isotropically, using cross sections read from the coherent file)
* container: [str]  sample container material. Thickness is .2 mm. Use NULL, Al or Nb.
* LSD: [m]          Distance sample-detector
* Frequency: [Hz]   Disk-chopper frequency
* LCC:       [m]    Distance between the Disk-chopper and the Fermi-chopper
* DELAY:     [s]    Fermi phase wrt. Disk-chopper opening
*
* %Link
* The BRISP spectrometer at the ILL <a href="www.ill.fr/brisp">BRISP</a>
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT ILL_BRISP(DM=3.355, string coh="V.lau", string inc="NULL", string container="NULL", LSD=4.5, Frequency=0, LCC=0, DELAY=0)

/* The DECLARE section allows us to declare variables or  small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
  double LCE      = 1.0;       /* distance fermi chopper <--> sample   */
  double LME;                  /* distance monochromator <--> sample   */

  double lambda   =0;
  double A1       =17.5;
  double RATIO    =4;

  double Ei;
  #pragma acc declare create(Ei)
  double lsd;
  double RH, RV;
  
%}

USERVARS %{
  double flag_container;
  double flag_sample;
  double flag_env;
  double vix;
  double viy;
  double viz;
  double Ef;
  double Q;
  double dE;
%}

/* The INITIALIZE section is executed when the simulation starts     */
/* (C code). You may use them as component parameter values.         */
INITIALIZE
%{
  double Ki, vi;
  double RPM    =14960;	       /* Fermi chopper rotation speed, Disk=RPM/RATIO [rpm] */

  LME = 5.540+LCE;
  LCC = LME-LCE-1.5;
  lsd=LSD;

  lambda= sin(A1*DEG2RAD)*2*DM;

  Ki = 2*PI/lambda;
  vi = K2V*fabs(Ki);
  Ei = VS2E*vi*vi;
  #pragma acc update device(Ei)

  DELAY   = LCC/vi;           /* neutron propagation time between choppers [s] */

  printf("Instrument ILL_BRISP:\n\tTime of Flight Neutron Spectrometer for Small Angle Inelastic Scattering\n");
  printf("Monochr. Bragg angle     [deg] %g\n",  A1);
  printf("Wavelength               [AA]  %g\n",  lambda);
  printf("Neutron velocity         [m/s] %g\n",  vi);
  printf("Incident Energy          [meV] %g\n",  Ei);

  printf("Travel time: Disk./Fermi [us]  %g\n",  DELAY*1e6);
  printf("CHOP Disk  Speed         [rpm] %g\n",  RPM/RATIO);
  printf("CHOP Fermi Speed         [rpm] %g        (%g [Hz], %g [rad/s])\n",
                                                 RPM, RPM/60, RPM/60*2*PI);
  printf("CHOP Disk  Phase         [deg] %g\n",  DELAY*RPM/60*360);
  printf("CHOP Fermi Phase         [deg] %g\n",  0.0);

  printf("Dist. Monochromator-Fermi [m]  %g\n",  LME-LCE);
  printf("Dist. Disk-Fermi          [m]  %g\n",  LCC);
  printf("Dist. Sample-Detector     [m]  %g\n",  LSD);

  Frequency=RPM/60; /* [rpm/60=rps=Hz]  */

  double L=LME;

  RV = 2*L*sin(DEG2RAD*A1);
  RH = 2*L/sin(DEG2RAD*A1);
%}

TRACE

COMPONENT Origin = Progress_bar()
  AT (0,0,0) ABSOLUTE ROTATED (-2*A1,0,0) ABSOLUTE

/* thermal tube IH3 17.5 deg inclinaison */
COMPONENT ThermalSource = Source_gen(
    radius = 0.1, focus_xw = 0.06, focus_yh = 0.06, dist=6,
    lambda0 = lambda, dlambda = lambda*0.01,
    T1=683.7,I1=0.5874e+13,T2=257.7,I2=2.5099e+13,T3=16.7 ,I3=1.0343e+12)
  AT (0, 0, 0) RELATIVE Origin

/* Soller Collimator                alpha1 = 0.4°
   700 mm length - 5 mm slit spacing
   16 Gd2O3 coated Kapton sheets, 75 µm thick
 */
COMPONENT C1 = Collimator_linear(xwidth=0.06, yheight=0.06,
          length=0.7, divergence=30,transmission=1.0)
AT (0,0,6) RELATIVE ThermalSource

COMPONENT MonoCradle1 = Arm()
AT (0,0,0.7+0.4) RELATIVE PREVIOUS

COMPONENT MonoCradle2 = Arm()
AT (0,0,0) RELATIVE PREVIOUS ROTATED (0,0,90) RELATIVE PREVIOUS

/* put back beam in horizontal pane */
/* Double focusing – horizontal variable focusing.
   Incident beam area at the monochromator: 6 x 6 cm2.
   Monochromator surface: 20.8 x 8.6 cm2
   20 crystals in a 4 x 5 matrix: crystal size 4 x 2 cm2, mosaic spread 0.4° (PG) and 0.25° (Cu)
 */

SPLIT COMPONENT Monok = Monochromator_curved(width=0.208, height=0.086, gap=0.0005,
          NH=4, NV=5, mosaich=30.0, mosaicv=30.0, r0=0.8, DM=DM)
AT (0,0,0) RELATIVE MonoCradle2 ROTATED (0,-A1,0) RELATIVE MonoCradle2

COMPONENT MonokOut = Arm()
AT (0,0,0) RELATIVE MonoCradle2 ROTATED (0,-2*A1,90) RELATIVE MonoCradle2

COMPONENT LmonokOUT = Monitor_nD(
    options = "lambda limits=[1.7 2.1] bins=40", xwidth=0.06, yheight=0.06)
  AT (0, 0, 1.7) RELATIVE MonokOut

COMPONENT Monitor1 = Monitor_nD(
    options = "x y, slit", xwidth=0.06, yheight=0.06)
  AT (0, 0, 0.01) RELATIVE PREVIOUS

/* Disk chopper: 8 slits width 6 cm on radius 0.21=8*16 deg=128 deg. Transmission 128/360=35 % */
COMPONENT DC=DiskChopper(radius=.240, theta_0=15, nslit=8,
  nu=Frequency/RATIO,yheight=0.06,isfirst=1)
  AT (0, 0, 1.73) RELATIVE MonokOut ROTATED (0,0,90) RELATIVE MonokOut
/* flux outgoing from Disk is in [n/s] (integrated in time) */

COMPONENT DC_OUT_T = Monitor_nD(
    options = "time limits=[5.5e-4 8.5e-4] bins=100", xwidth=0.06, yheight=0.06)
AT (0,0,0.2) RELATIVE PREVIOUS ROTATED (0,0,0) RELATIVE MonokOut

/* Honeycomb Collimator       alpha2 = 0.4°
- 2000 mm length, Gd2O3-coated aluminum
- honeycomb arrangement of converging tubes of hexagonal section
- 3 different collimators for convergence at 2, 4(recommended), 6 m from the sample position
- collimators are maintained under vacuum and selected by means of a 4-sector revolver
- one sector is used for coarse collimation using 2 Cd diaphragms (in: 81x81mm2, out: 67x67mm2)
 */
COMPONENT HC = Guide_honeycomb(
    w1 = .085, w2 = .066, l = 2, R0 = 0, nslit = 6)
  AT (0, 0, 0.9) RELATIVE DC

COMPONENT HC_OUT_PSD = Monitor_nD(
    options = "x y, all bins=50", xwidth=0.07, yheight=0.07)
  AT (0, 0, 2+0.1) RELATIVE PREVIOUS

/* this monitor has an integrated flux in [n/s] */
/* shows only one pulse from Disk chopper, but standing for all slits in a period
 * In reality, per pulse we must divide by nu=250.0/4
 */
COMPONENT FC_IN_T = Arm()
  AT (0, 0, LME-LCE-0.2) RELATIVE MonokOut
EXTEND %{
  /* we must make so that integrated flux is in [n/s] */
  PROP_Z0;
  /* we fisrt put all events in the Fermi time transmission window */
  /* Fermi opening time is atan(w/length)/PI/frequency */
  double Tfermi=atan2(0.031/136,0.013)/INSTRUMENT_GETPAR(Frequency)/PI; /* opening time for Fermi [s] */
  t  = INSTRUMENT_GETPAR(DELAY)*(INSTRUMENT_GETPAR(LCC)-0.2)/INSTRUMENT_GETPAR(LCC)+(rand01()-0.5)*Tfermi;    /* center time on DELAY at Fermi position */

  /* geometric integrated transmission of Fermi=0.55 % */
  /* 2 sides*angular opening = 2*atan2(0.031/136,0.013)*RAD2DEG/360 */
  p *= 2*atan2(0.031/136,0.013)*RAD2DEG/360;

%}

/* Dist: Monochromator – Fermi chopper 5540 mm */
COMPONENT FC=FermiChopper(xwidth=.031, yheight=.101, nslit=136, nu=Frequency,
  phase=0, radius=.135/2,length=0.013,delay=DELAY)
  AT (0, 0, LME-LCE) RELATIVE MonokOut ROTATED (0,0,90) RELATIVE MonokOut
/* Dist: Fermi chopper – Sample 1000 mm */

COMPONENT FC_OUT_T = Monitor_nD(
    options = "time limits=[2.15e-3 2.175e-3] bins=20", xwidth=0.06, yheight=0.06)
  AT (0, 0, 0.2) RELATIVE PREVIOUS ROTATED (0,0,0) RELATIVE MonokOut


/* Sample position ********************************************************** */
/* external shield Al, container and sample in concentric geometry            */
COMPONENT SampleXY=Monitor_nD(
  xwidth=.03, yheight=.03,
  options="x y, per cm2", restore_neutron=1)
  AT (0, 0, LME) RELATIVE MonokOut

SPLIT COMPONENT SamplePos=Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
EXTEND %{
  flag_container=flag_sample=flag_env=0;
  vix=vx; viy=vy; viz=vz;
%}

COMPONENT Environment_in=Isotropic_Sqw(
  radius = 0.05, yheight = 0.1, thickness=0.002,
  Sqw_coh="Al.lau", concentric=1, d_phi=15, order=1,
  p_interact=0.05)
  AT (0, 0, 0) RELATIVE SamplePos
EXTEND %{
  flag_env=SCATTERED;
%}

COMPONENT Container_in=Isotropic_Sqw(xwidth=0.05+2e-4, yheight=0.05+2e-4, zdepth=4e-3+2e-4,
  thickness=2e-4, Sqw_coh=container, concentric=1, d_phi=15, order=1,
  p_interact=0.05)
  AT (0, 0, 0) RELATIVE SamplePos
EXTEND %{
  flag_container=SCATTERED;
%}

COMPONENT Sample=Isotropic_Sqw(xwidth=0.05, yheight=0.05, zdepth=4e-3,
  Sqw_coh=coh, Sqw_inc=inc, p_interact=0.99, order=1, d_phi=30)
  AT (0, 0, 0) RELATIVE SamplePos
EXTEND %{
  if (SCATTERED)
    flag_sample=SCATTERED*(VarSqw.type == 'c' ? 1 : -1);
%}

COMPONENT Container_out=COPY(Container_in)(concentric=0)
  AT (0, 0, 0) RELATIVE SamplePos
EXTEND %{
  if (SCATTERED) flag_container=1;
%}

COMPONENT Environment_out=COPY(Environment_in)(concentric=0)
  AT (0, 0, 0) RELATIVE SamplePos
EXTEND %{
  if (SCATTERED) flag_env=1;
%}

/* End Sample position ****************************************************** */

/* Detector model *********************************************************** */
/* Dist: Sample – Detector 4500 mm */
COMPONENT DetectorPos=Arm()
AT (0, 0, LSD) RELATIVE SamplePos
EXTEND %{
  PROP_Z0;
  /* masks: 4 absorbing panels on corners: width=.203 height=(1118-600)/2=.259 */
  if (fabs(x) > (160*12.7/1000)/2-0.203 && fabs(y)>(1.118/2-.259) ) ABSORB;
  if (x*x+y*y < 0.05*0.05) ABSORB; /* central beamstop: radius 0.05 */

  Ef = VS2E*(vx*vx+vy*vy+vz*vz);
%}

/* Global monitor */
COMPONENT Detector_pre = Arm()
AT(0,0,0) RELATIVE PREVIOUS
EXTEND
%{
  Q=V2K*sqrt( (vx-vix)*(vx-vix)+(vy-viy)*(vy-viy)+(vz-viz)*(vz-viz) );
  dE=Ef-Ei;
%}

COMPONENT Detector=Monitor_nD(  /* SANS XY detector */
  options="x y", xwidth=160*12.7/1000, yheight=1.118, bins=50, restore_neutron=1)
  WHEN (flag_sample || flag_container || flag_env)
  AT (0, 0, 0.01) RELATIVE PREVIOUS

COMPONENT Detector_Q=Monitor_nD( /* I(q) detector */
  options="user1 limits=[0 1.5]", user1="Q",username1="Q transfer [Angs-1]", xwidth=160*12.7/1000, yheight=1.118, bins=50, restore_neutron=1)
  WHEN (flag_sample || flag_container || flag_env )
  AT (0, 0, 0.01) RELATIVE PREVIOUS

COMPONENT Detector_E=Monitor_nD( /* Energy distribution */
  options="user1", xwidth=160*12.7/1000, yheight=1.118, restore_neutron=1, bins=100,
  min=-Ei, max=2*Ei,
  user1="dE", username1="Energy transfer [meV]")
  WHEN (flag_sample || flag_container || flag_env)
  AT (0, 0, 0.01) RELATIVE PREVIOUS

COMPONENT Detector_QE=Monitor_nD(
  options="user1 limits=[0 1.5], user2 limits=[-10 10]", xwidth=160*12.7/1000, yheight=1.118,
  restore_neutron=1, bins=50,
  user1="Q", username1="Q transfer [Angs-1]",
  user2="dE", username2="Energy transfer [meV]")
WHEN (flag_sample || flag_container || flag_env)
AT (0, 0, 0.01) RELATIVE PREVIOUS

COMPONENT EndOfInstrument=Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
EXTEND %{
  /* make sure this is the end of the simulation (do not propagate to following volumes)*/
  ABSORB;
%}

/* Additional cosmetic shapes =============================================== */
/* reactor vessel: radius 1.25 m  . Core 0.220*/
/* Beam tube IH3 Phi=.108 at dx=-.454 dz=.279 from reactor core , inclined 17.5 deg */
/* put large detector housing phi=2.8 m to be pretty l=2-6 m from sample */

COMPONENT DetectorTube = Shape(radius=1.4, yheight=8.0)
AT (0,0,0) RELATIVE Detector ROTATED (90,0,0) RELATIVE Detector

COMPONENT ToMonok = Arm() /* to avoid misleading lines */
AT (0,0,0) RELATIVE MonoCradle1

COMPONENT ToOrigin = Arm() /* to avoid misleading lines */
AT (0,0,0) ABSOLUTE

COMPONENT Core = Shape(radius=.202, yheight=1.0)
AT (.454, 0, -.279) ABSOLUTE

COMPONENT Vessel = Shape(radius=1.25, yheight=2.0)
AT (0,0,0) RELATIVE PREVIOUS

END