/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Instrument: Light_H15_IN6
*
* %Identification
* Written by: <a href="mailto:farhi@ill.fr">Emmanuel Farhi</a>
* Date: 17th Jan 2005.
* Origin: <a href="http://www.ill.fr">ILL (France)</a>
* %INSTRUMENT_SITE: ILL
*
* The IN6 Time-of-Flight simulation at the ILL (instrument only).
*
* %Description
*
* IN6 is a time focussing time-of-flight spectrometer designed for quasielastic and
* inelastic scattering for incident wavelengths in the range of 4 to 6 Angs.
*
* An intense beam is extracted from the H15 guide by a vertically focussing
* monochromator array. It consists of three composite pyrolytic graphite
* monochromators using the full height (20 cm) of the guide and focussing the beam
* at the sample position. In order to minimise the interference with the
* subsequent instruments, the monochromator can deliver only four wavelengths:
* 4.1; 4.6; 5.1; and 5.9 Angs. The second order reflection from the graphite
* monochromator is removed by a beryllium-filter cooled at liquid nitrogen
* temperature.
* To achieve the time-focussing condition, the beam is pulsed by a Fermi chopper.
* It has a small slot length to ensure a good transmission. The normal distance
* between the Fermi chopper and the sample is 38 cm. To prevent frame-overlap when
* the chopper is rotating faster than 7500 rpm, a suppressor chopper is placed
* before the Fermi chopper and rotates in phase with the latter.
*
* The secondary spectrometer consists first of an evacuated sample area. The
* detector bank is entirely covered with detector boxes, thus avoiding the
* inconvenience of moving the counters.
*
* This instrument model contains a cold source a triple monochromator
* (using the GROUP), two Fermi Choppers (including one background chopper), a
* liquid sample handling coherent and incoherent processes (elastic and inelastic)
* with multiple scattering, customized monitors, and the SPLIT mechanism to
* improve the statistics. The H15 guide is not described in this model.
*
* %Example: lambda=4.14 Detector: M_theta_t_all_I=120000
*
* %Parameters
* lambda: [Angs]                          wavelength within 4.14|4.6|5.12|5.92
* dlambda: [Angs]                         wavelength HALF spread. default is 0.075
* SPEED: [rpm]                            Fermi chopper speed.  -1=auto, 0=stopped in open pos.
* RATIO: [1]                              Suppressor speed ratio. -1=no suppressor.
* PHASE: [deg]                            Fermi phase w/r/ to Suppressor. -360=auto
* M1: [coder values]                      monochromator motor 1 position. -1=auto
* M2: [coder values]                      monochromator motor 2 positinn. -1=auto
* M3: [coder values]                      monochromator motor 3 position. -1=auto
* MONITOR: [something    like time in s]  monitor preset
* CHA_WIDTH: [us]                         channel width. -1=auto
* TOF_CHA_RESOL: [1]                      number of channels.
* TOF_DELAY: [us]                         TOF delay. -1=auto
* ELPEAK: [1]                             elastic peak channel. -1=auto
* mFC: [1]                                super mirror FermiChopper coating m-value
* Sqw_coh: [str]                          sample coherent S(q,w) file name
* Sqw_inc: [str]                          sample incoherent S(q,w) file name
* radius: [m]                             outer radius of sample hollow cylinder
* thickness: [m]                          thickness of sample hollow cylinder
*
* %Link
* <a href="http://www.ill.fr/in6">The IN6@ILL Yellow Pages</a>
* %Link
* R.Scherm et al, "Another Time of Flight Spectrometer", ILL Report 76S235, 1976
* %Link
* R.Scherm, "A high-resolution spectrometer ...", report Jul-295-NP, Kernforschungsanlage Julich, 1965
* %Link
* Y.Blanc, "Le spectrometre a temps de vol IN6", ILL Report 83BL21G, 1983
* %Link
* K.H.Beckurts et al, Neutron physics, Kernforschungsanlage Karlsruhe, 1964 (p317)
* %Link
* R.Scherm and T.Springer, "Proposal of a multiple Chopper...", Kernforschungsanlage Julich, 19xx
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT ILL_IN6(lambda=4.14, dlambda=0.075, SPEED=-1, M1=-1, M2=-1, M3=-1, MONITOR=1, CHA_WIDTH=-1, TOF_DELAY=-1, TOF_CHA_RESOL=128, ELPEAK=-1, RATIO=1, mFC=0, PHASE=-360, string Sqw_coh="Rb_liq_coh.sqw", string Sqw_inc="Rb_liq_inc.sqw", radius=0.01, thickness=0.005)

DECLARE
%{

/* capture flux positions from moderator: 21.4    28.4    61.2 */

/* variables for IN6 */
  double DM;        /*   mono d-spacing (Angs)  */
  double mos;
  double RV;

/* Bragg angles of the 3 monochromators */
  double A1;
  double A2;
  double A3;
  double LME;   /* distance monochromator 2 <--> sample */
  double LMM; /* distance between 2 monochromators    */
  double LED; /* distance sample <--> detector        */
  double LCE; /* distance fermi chopper <--> sample   */
  double LCC;   /* [m] Chopper1-Chopper2 distance */
  double Frequency, vi;
  #pragma acc declare create(Frequency)
  double ref_phas, phas_ferm;

  double iTOF_DELAY;   /* time of arrival at sample position, from source */
  double iCHA_WIDTH;
  double iTOF_CHA_RESOL;
  double iELPEAK;
  double iRATIO;
  #pragma acc declare create(iRATIO) 
  double iSPEED;
  double A2cradle;
  double iPHASE, period;

/* flags per event type */
  char flag_ci, flag_co, flag_ct; /* cryo-in/out container */
  char flag_single, flag_multi;

/* monitoring sample env */
  double ki_x, ki_y, ki_z;
  double kf_x, kf_y, kf_z;
  double dq, dw, vf;
  char optL[256];
  char optT[256];
  char opt1[256]; /* options for Monitor_nD */
  char opt2[256];
  char opt3[256];
%}

USERVARS
%{
  double vix;
  double viy;
  double viz;
  double monok_index;
%}

INITIALIZE
%{
  /* capture flux positions from moderator: 21.4    28.4    61.2 */

/* variables for IN6 */
  DM = 3.355;        /*   mono d-spacing (Angs)  */
  mos= 40;
  RV = 3;

/* Bragg angles of the 3 monochromators */
  LME = 2.1;   /* distance monochromator 2 <--> sample */
  LMM = 0.030; /* distance between 2 monochromators    */
  LED = 2.483; /* distance sample <--> detector        */
  LCE = 0.395; /* distance fermi chopper <--> sample   */
  LCC = 0.2;   /* [m] Chopper1-Chopper2 distance */
  ref_phas=0;
  dq=0;
  dw=0;

  double chopper_const = 252.77;  /*constant in chopper SPEED formula*/
  double Ki, Ei, theta;
  double tmin, tmax;
  double dE  = 0.0;   /* energy transfer */

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

/* IN6: calculate theta angles for 3 monochromators*/
  theta = asin(lambda/DM/2);
  A2 = theta*2;

  A1 = atan2(LME*sin(A2),(LME*cos(A2)+LMM))*RAD2DEG;
  A3 = atan2(LME*sin(A2),(LME*cos(A2)-LMM))*RAD2DEG;
  A2 *=RAD2DEG;

  A2cradle = A2;

  RV = 2*LME*sin(theta);

  if (A1<0.0) A1=180+A1;
  if (A2<0.0) A2=180+A2;
  if (A3<0.0) A3=180+A3;

  if (M1 == 0) A1 = 0; else
    if (M1>=0) A1 = -0.0210199*M1+178.55; else  M1 = -(A1-178.55)/0.0210199;
  if (M2 == 0) A2 = 0; else
    if (M2>=0) A2 = -0.0210302*M2+182.558; else M2 = -(A2-182.558)/0.0210302;
  if (M3 == 0) A3 = 0; else
    if (M3>=0) A3 = -0.0206945*M3+187.566; else M3 = -(A3-187.566)/0.0206945;

/* IN6: compute Tof settings from Light.Custom.Light_Custom_IN6_Calc_TOF_Choppers */
  {
    double el_t_resol  = 0.125;   /* [us] Electronic Time Base */
    ref_phas    = 0;      /* [deg] Reference Phase */
    double phase_offset= 0;       /* [deg] Phase Offset (added to Fermi phase) */
    double el_delay    = 44.875;  /* [us]  Default Electronic Delay */

    double speed, chan_width, dead_time, time_of_flight, trav_time;
    double delta_phase, el_peak_O, delay;

    if (TOF_CHA_RESOL<=0) iTOF_CHA_RESOL=128; else iTOF_CHA_RESOL=TOF_CHA_RESOL;
    if (RATIO <= 0)       iRATIO = 1;         else iRATIO        =RATIO;
    #pragma acc update device(iRATIO)
    if (ELPEAK >= 0 && ELPEAK<=iTOF_CHA_RESOL) iELPEAK  = ELPEAK;
    else iELPEAK=ceil(iTOF_CHA_RESOL/2);

    speed       = 60*K2V/(DM*cos(theta)*(LCE+(LED*pow((1-dE/Ei),-1.5))));
    period      = 0.5e6 * 60 * iRATIO/speed;

    chan_width  = floor(period/el_t_resol/iTOF_CHA_RESOL)*el_t_resol;
    dead_time   = period-(iTOF_CHA_RESOL*chan_width);
    time_of_flight = (0.07+LCC+LCE+LED)/vi*1e6;
    trav_time   = LCC/vi*1e6;
    delta_phase = (trav_time/period)*180;
    phas_ferm   = ref_phas + delta_phase;
    if (fmod(iRATIO, 2) == 0) phas_ferm *= 2;
    phas_ferm  += phase_offset;
    el_peak_O   = floor((time_of_flight + el_delay)/chan_width);
    delay       =  (el_peak_O-iELPEAK) * chan_width;
    if (iELPEAK >= el_peak_O) delay +=  period;
    if (delay <= 1) delay = 2;

    if (PHASE>-180 && PHASE <360) iPHASE=PHASE;    else iPHASE=-phas_ferm;
    if (CHA_WIDTH <=0)    iCHA_WIDTH=chan_width;   else iCHA_WIDTH=CHA_WIDTH;
    if (TOF_DELAY <=0)    iTOF_DELAY=delay;        else iTOF_DELAY=TOF_DELAY;
    if (SPEED     <0)     iSPEED=speed;            else iSPEED=SPEED;
    Frequency = iSPEED/60;
    #pragma acc update device(Frequency)
    printf("Instrument Simulation %s (%s)\n", instrument_name, instrument_source);
    printf("  using computed monochromator take-off angles: %g %g %g [deg]\n", A1, A2, A3);
    printf("Wavelength               [AA]  %g\n",  lambda);
    printf("Neutron velocity         [m/s] %g\n",  vi);
    printf("Monochr. Bragg angle     [deg] %g\n",  A2/2);
    printf("Incident Energy          [meV] %g\n",  Ei);
    printf("Focusing Energy transfer[meV] %g\n",  dE);
    printf("Travel time: Supp./Fermi [us]  %g\n",  trav_time);
    printf("Travel time: Supp./Det.  [us]  %g\n",  time_of_flight);
    printf("TOF Delay                [us]  %g\n",  delay);
    printf("TOF Dead Time            [us]  %g\n",  dead_time);
    printf("TOF Period (1 cycle)     [us]  %g\n",  period);
    printf("TOF Channel width        [us]  %g\n",  chan_width);
    printf("CHOP Fermi Phase         [deg] %g\n",  iPHASE);
    printf("CHOP Suppressor Phase    [deg] %g\n",  ref_phas);
    printf("CHOP Fermi Speed         [rpm] %g\n",  iSPEED);
    printf("CHOP Suppressor Speed    [rpm] %g\n",  iSPEED/iRATIO);
    printf("Number of time channels        %g\n",  iTOF_CHA_RESOL);
    printf("Current Elastic Peak Ch.       %g\n",  iELPEAK);
    printf("Elast. peak ch. for 0-delay    %g\n",  el_peak_O);

    /* chopper to detector */
    tmin = time_of_flight*1e-6 - (iCHA_WIDTH*iELPEAK-iTOF_DELAY)*1e-6;
  }


 /* distance to cover to detector from chopper: LCE+LED
      Center time on 0 at Fermi center
      LCE     from chopper to sample pos
      LED     from sample to detector
      propagation time t_p =(LCE+LED)/vi;
      falls on ELPEAK channel.
      Tmin = t_p-iCHA_WIDTH*1e-6*iELPEAK
      Tmax = Tmin +N_CHan...
   */


  tmax = tmin+iTOF_CHA_RESOL*iCHA_WIDTH*1e-6;

  printf("Time window:                   min=%g max=%g delay=%g tof-width=%g [ms]\n", tmin*1000, tmax*1000, iTOF_DELAY*1e-3, (tmax-tmin)*1000);

  sprintf(opt2, "kxy limits=[0 5] bins=50, energy limits=[%g %g] bins=40, banana, parallel",
  (Ei-20 < 0 ? 0 : Ei-20), Ei+20);
  sprintf(opt1, "angle limits=[0 180] bins=180, energy limits=[%g %g] bins=40, banana, parallel",
  (Ei-20 < 0 ? 0 : Ei-20), Ei+20);

  sprintf(opt3,"user1 limits=[0.5,3.5] bins=9, lambda limits=[%g %g] bins=20, square, per cm2",0.97*lambda,1.03*lambda);

  sprintf(optL,"lambda limits=[%g %g] bins=50",lambda-dlambda,lambda+dlambda);
  sprintf(optT,"t slit limits=[%g %g]",tmin, tmax);
%}

/* -------------------------------- TRACE -------------------------------- */
TRACE

REMOVABLE COMPONENT Origin = Progress_bar()
  AT (0,0,0) ABSOLUTE

REMOVABLE COMPONENT VCS = Source_gen(
  yheight  = 0.2,
  xwidth   = 0.03,
  focus_xw = 0.03,
  focus_yh = 0.2,
  lambda0  = lambda,
  dlambda  = dlambda,
  T1=216.8,I1=1.24e+13,	/* VCS parameters */
  T2=33.9, I2=1.02e+13,
  T3=16.7 ,I3=3.0423e+12,
  verbose  = 1)
  AT (0, 0, 0) RELATIVE Origin

REMOVABLE COMPONENT SourceTarget = Arm()
AT (0,0,2.55) RELATIVE PREVIOUS

/* ----------------------- IN6 Monochromators GROUP ----------------------- */

SPLIT COMPONENT Cradle = Arm()
AT (0,0,0.15) RELATIVE PREVIOUS
EXTEND
%{
  monok_index=0;
%}
/* triple-monochromator description:
 *   7 blades, vertically focusing RV=3 m, fixed.
 *   Each blade is 54 mm width, 29 mm heigh. Vertical angle +/- 3 deg.
 *   mosaic 23 to 40 min. Motors 0.012 deg/step
 *   distance between each crystal ensemble 4 cm
 */

COMPONENT Mono1 = Monochromator_curved(
  RV     = RV, NV = 7,NH=1,
  zwidth = 0.054, yheight = 0.029,
  DM     = 3.355, gap = 0.001,
  mosaic = 40, r0=1, t0=1,
  reflect="HOPG.rfl", transmit="HOPG.trm")
AT (0,0, -LMM) RELATIVE Cradle ROTATED (0,A1/2,0) RELATIVE Cradle
EXTEND
%{
  if (SCATTERED) { monok_index=1; }
%}

COMPONENT Mono2 = Monochromator_curved(
  RV     = RV, NV = 7,NH=1,
  zwidth = 0.054, yheight = 0.029,
  DM     = 3.355, gap = 0.001,
  mosaic = 40, r0=1, t0=1,
  reflect="HOPG.rfl", transmit="HOPG.trm")
AT (0,0, 0) RELATIVE Cradle ROTATED (0,A2/2,0) RELATIVE Cradle
EXTEND
%{
  if (SCATTERED) { monok_index=2; }
%}

COMPONENT Mono3 = Monochromator_curved(
  RV     = RV, NV = 7,NH=1,
  zwidth = 0.054, yheight = 0.029,
  DM     = 3.355, gap = 0.001,
  mosaic = 40, r0=1, t0=1,
  reflect="HOPG.rfl", transmit="HOPG.trm")
AT (0,0, +LMM) RELATIVE Cradle ROTATED (0,A3/2,0) RELATIVE Cradle
EXTEND
%{
  if (SCATTERED) { monok_index=3; }
%}

/* sample position direction  */
COMPONENT mono_out = Arm()
AT (0,0,0) RELATIVE Cradle ROTATED (0,A2cradle,0) RELATIVE Cradle

/* --------------------------- IN6 Suppressor ------------------------ */

COMPONENT SuppPos = Arm()
AT (0,0,LME-LCE-LCC) RELATIVE mono_out

COMPONENT Mon_SuppInL = Monitor_nD(
  xwidth = 0.052, yheight = 0.098,
  options=optL)
AT (0,0,-0.07-0.002) RELATIVE SuppPos

COMPONENT Mon_SuppInT = Monitor_nD(
  xwidth = 0.052, yheight = 0.098,
  options=optT, bins=iTOF_CHA_RESOL)
AT (0,0,-0.07-0.001) RELATIVE SuppPos
EXTEND
%{
  double Vi=sqrt(vx*vx+vy*vy+vz*vz); /* K2V*2*PI/sLambda; */
  double ratio=PI*Frequency/iRATIO/atan(0.052/0.14);
  //The opening angle for one neutron times number of openings to total closed angle.
  if (iRATIO && Vi) {
    /* we fisrt put all events in the Fermi time transmission window */
    /* Fermi opening time is atan(w/length)/PI/frequency */
    double Tfermi=atan2(0.052,0.14)/Frequency/PI; /* opening time for Fermi [s] */
    t  = -(0.071/Vi)+(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 *= 4*atan2(0.052,0.14)*RAD2DEG/360;
  }
%}

/* Suppressor Chopper position. */
COMPONENT Suppressor = FermiChopper(radius=0.07, nu=Frequency/iRATIO,
   yheight=0.098, xwidth=0.052, nslit=1, R0=0, phase=0,
   length=0.012, eff=1, verbose=1)
  WHEN (iRATIO > 0)
AT (0,0,0) RELATIVE SuppPos

COMPONENT Mon_SuppOutT = Monitor_nD(
  xwidth = 0.052, yheight = 0.098,
  options="t slit limits=[-700e-6 700e-6]", bins=iTOF_CHA_RESOL)
AT (0,0,+0.07+0.001) RELATIVE SuppPos

/* --------------------------- IN6 Fermi ------------------------ */

COMPONENT SlitFC = Slit(
    xwidth = 0.0441, yheight = 0.0641)
  AT (0, 0, LME-LCE-0.041) RELATIVE mono_out

COMPONENT FermiPos = Arm()
AT (0,0,LME-LCE) RELATIVE mono_out

COMPONENT FermiM = FermiChopper(phase=iPHASE, radius=0.04, nu=Frequency,
   yheight=0.064, xwidth=0.044, nslit=200.0, R0=.99,
   Qc=(mFC < 1 && mFC ? mFC*0.02176 : 0.02176), alpha=2.33, m=mFC, length=0.012, eff=1.0, verbose=1)
AT (0,0,0) RELATIVE FermiPos

COMPONENT Mon_FermiOutdT = Monitor_nD(
  xwidth = 0.044, yheight = 0.064,
  options="t slit limits=[180e-6 365e-6]", bins=iTOF_CHA_RESOL)
AT (0,0,+0.06+0.001) RELATIVE FermiPos

/* --------------------------- IN6 Fermi END --------------------- */

COMPONENT MonokMonitor = Monitor_nD(
    xwidth = 0.2, yheight = 0.2, user1="monok_index", username1="Monoch. index",
    options=opt3)
AT (0,0,+0.06+0.002) RELATIVE FermiPos

/* sample position (at 2.1 m from monoks) */

COMPONENT Mon_SampleInT = Monitor_nD(
  xwidth = 0.05, yheight = 0.05,
  options=" t limits=[255e-6 425e-6] parallel, per cm2", bins=iTOF_CHA_RESOL)
AT (0,0,LME-.273) RELATIVE mono_out

COMPONENT Mon_SampleInXY = Monitor_nD(
  xwidth = 0.06, yheight = 0.1,
  options="x y parallel, per cm2 bins=50")
AT (0,0,0) RELATIVE PREVIOUS

/* BEGIN ********************************** Sample environment and sample */

COMPONENT sample_pos = Arm()
AT (0,0,LME) RELATIVE mono_out
EXTEND %{
  vix=vx;
  viy=vy;
  viz=vz;
%}
  
SPLIT COMPONENT Sample=Isotropic_Sqw(
  radius = radius, thickness=thickness, yheight = 0.055,
  Sqw_coh=Sqw_coh, Sqw_inc=Sqw_inc, p_interact=0.9
) AT (0, 0, 0) RELATIVE sample_pos
EXTEND
%{
  if (!SCATTERED) ABSORB;
%}

COMPONENT Out = PSD_monitor_4PI(filename="out")
  AT (0,0,0) RELATIVE sample_pos

COMPONENT M_theta_t_all = Monitor_nD(
 xwidth=2.5, yheight=1,
 options=opt1,
 bins=100)
AT (0,0,0) RELATIVE sample_pos

COMPONENT Detector_nM = Sqw_monitor(
 radius=1.3, yheight=1,
 nq=100, qmin=0, qmax=6,
 nE=100, Emin=-20, Emax=20,
 vix="vix", viy="viy", viz="viz",
 filename="Detector_nM")
 AT (0, 0, 0) RELATIVE PREVIOUS

END