/********************************************************************
*
* Instrument: ILL_H16_IN5
*
* %Identification
* Written by: E. Farhi, J. Ollivier, Celia Castan Guerrero
* Date: Jan 2004-2009
* Origin: ILL
*
* %INSTRUMENT_SITE: ILL
*
*   The IN5B instrument: chopper system + sample + PSD and tof detector
*
* %Description
*
*  The IN5@ILL TOF spectrometer from chopper system to final detector, with sample.
*  The detector model includes Fe housing and tube cross talk absorbing masks
*  with angle restriction (neutrons that scatter in Fe in front of a tube and
*  enter a different tube are absorbed). This model does not include the H16 guide.
*
* %Example: lambda=4.5 Detector: Det_PSD_I=1.3e6
*
* %Parameters
* lambda: [AA]       mean incident wavelength
* dlambda: [AA]      wavelength half width.
* speed: [rpm]       chopper speed (60*frequency)
* ratio: [1]         velocity of chopper3 = velocity of chopper1 * ratio
* housing: [string]  material used as detector housing. Use 0 or NULL for none.
* radius: [m]        radius of sample (outer).
* thickness: [m]     thickness of sample. 0=filled
* height: [m]        height of sample. 0=sphere
* coh: [string]      sample coherent data file or NULL. Use powder LAZ/LAU or SQW file.
* inc: [string]      sample incoherent Sqw data file or NULL
* order: [1]         order for scattering in sample. O=all, 1=single
* %Link
* The <a href="http://www.ill.eu/in5">IN5@ILL</a> cold time of flight instrument
*
* %E
************************************************************************/
DEFINE INSTRUMENT ILL_IN5(lambda=4.5,dlambda=0.05,speed=8500, ratio=0.5,
  string housing="Fe.laz", string coh="Y3Fe5O12_YIG.laz", string inc="NULL",
  thickness=0, height=0.025, radius=0.005, order=0)

DECLARE
%{

//---Guide (1st part of guide) --------------------------------------
  double L_Guide1, L_Guide21;
  double L_Guide22, L_Guide23;             // gerade dimensions
  double L_Guide3, L_Guide41;
  double L_Guide42, L_Guide43;
  double L_Guide44, L_Guide45;
  double L_Collimator, L_CollSample;
  double L_gap, disk_gap,mono_gap;
  //---Reactor & Krumm Guide coating (2nd part of guide) ------------
  double alt_Guide_Qc,alt_Guide_Ro,alt_Guide_alpha,alt_Guide_W;
  //---Neue Guide coating--------------------------------------------
  double Guide_Qc,Guide_Ro,Guide_alpha,Guide_W;
  //---Choppers-------------------------------------------------------
  double Ch_mean_R[7];               // <R> = rotation axis - guide axis
  double Ch_width[7],Ch_height[7];   // slots widht and height
  double Ch_alpha[7];                // angular aperture of choppers
  double Ch_phase[7];                // matching (phase) time for choppers
  double Ch_phase_angle[7];          // matching (angle) for choppers
  double Ch_Ltot[7];                 // distance - ref_time chopper
  double Ch_Vp[7];                   // angular velocity of choppers
  double disk_N;                     // slot number on the disks

//---Sample-by-itself---------------------------------------------------
  double L_sample;
//---Detector PSD-------------------------------------------------------
  double ang_ini, ang_fin, det_angle;
%}

INITIALIZE
%{
  int i;
  double v0;                         // neutron mean velocity

//==========================================================================
//                 Guide
//==========================================================================
  L_gap     = 0.2130;     // gap VTE+OT-H16
  L_Guide1  = 4.3900;     // for gerade Guide1
  L_Guide21 = 0.6950;     // for gerade Guide21
  L_Guide22 = 0.1300;     // for gerade Guide22
  L_Guide23 = 0.69500;    // for gerade Guide23
  disk_gap  = 0.02;       // full gap at choppers
  L_Guide3  = 5.5125;     // for gerade Guide3
  L_Guide41 = 0.7425;     // for gerade Guide41
  L_Guide42 = 0.0350;     // for gerade Guide42
  L_Guide43 = 0.7500;     // for gerade Guide43
  L_Guide44 = 0.0350;     // for gerade Guide44
  L_Guide45 = 0.7900;     // for gerade Guide45
  mono_gap  = 0.0300;     // gap for the 1st monitor
  L_Collimator= 0.1300;           // for gerade Collimator
  L_CollSample= 0.2400-0.025;     // the sample chamber size & keep

  printf("Instrument: ILL_IN5 (IN5 disk chopper time-of-flight spectrometer).\n  Wavelength lambda=%g [Angs]\n",
    lambda);

  // Alt Guide coating
  alt_Guide_Qc    = 0.021745; // for m=1 alpha and W aren't used.
  alt_Guide_Ro    = 0.995;
  alt_Guide_alpha = 6.07;
  alt_Guide_W     = 0.0023;
  // New Guide and super-mirors
  Guide_Qc    = 0.02275;
  Guide_Ro    = 0.996;
  Guide_alpha = 5.75;
  Guide_W     = 0.00125;

//==========================================================================
//                 Choppers
//==========================================================================
  Ch_mean_R[0]  =  0.285;            //
  Ch_mean_R[1]  =  0.285;            //
  Ch_mean_R[2]  =  0.285;            //
  Ch_mean_R[3]  =  0.299;            // <R> = rotation axis - guide axis
  Ch_mean_R[4]  =  0.299;            //
  Ch_mean_R[5]  =  0.304;            //
  Ch_mean_R[6]  =  0.304;            //

  Ch_height[0]  =  0.2;              //
  Ch_height[1]  =  0.17   ;          //
  Ch_height[2]  =  0.16813;          //
  Ch_height[3]  =  0.081;            // Height of the disk which "see" the guide
  Ch_height[4]  =  0.08031;          //
  Ch_height[5]  =  0.07069;          //
  Ch_height[6]  =  0.0700;           //

  Ch_alpha[0]   =  9.0;              //
  Ch_alpha[1]   =  9.0;              //
  Ch_alpha[2]   =  9.0;              //
  Ch_alpha[3]   =  9.5;              // angular apperture of choppers [degrees]
  Ch_alpha[4]   =  9.5;              //
  Ch_alpha[5]   =  3.25;             //
  Ch_alpha[6]   =  3.25;             //

  disk_N  = 2;

  for (i=1;i<=6;i++){
     Ch_Vp[i]    = 0.0;
     Ch_Ltot[i]  = 0.0;
     //printf("Ch%d: Rmin = %f, alpha = %f\n",i,Ch_mean_R[i],Ch_alpha[i]);
     //printf("Alpha_guide at Ch. %d = %f deg\n",i,Gu_alpha[i]*180/PI);
  }

  if (speed==0){
    printf("FATAL ERROR: Chopper speed = 0 !");
    exit(-1);
  }

  // set the choppers speed [rad/s]
  Ch_Vp[0]   =  speed*2*PI/60;
  Ch_Vp[1]   =  speed*2*PI/60;
  Ch_Vp[2]   =  speed*2*PI/60;
  Ch_Vp[3]   =  speed*2*PI/60*ratio;
  Ch_Vp[4]   =  speed*2*PI/60;
  Ch_Vp[5]   =  speed*2*PI/60;
  Ch_Vp[6]   =  speed*2*PI/60;

  v0  = 3956.035/lambda;

  //----------------------------------------------------------------
  // Compute the phases of each choppers
  //-------------------------------------
  // Zero time at chopper 0
  // 1st compute the distance from the zero time position for each chopper
  // 2nd compute the phase as distance/velocity, it means, the time delay
  //----------------------------------------------------------------


  Ch_Ltot[0] = 0;
  Ch_Ltot[1] = L_gap+L_Guide1+0.0003+L_Guide21+disk_gap/2.0;
  Ch_Ltot[2] = Ch_Ltot[1]+disk_gap+L_Guide22;
  Ch_Ltot[3] = Ch_Ltot[2]+disk_gap+L_Guide23+L_Guide3+L_Guide41+2*0.0003;
  Ch_Ltot[4] = Ch_Ltot[3]+disk_gap+L_Guide42;
  Ch_Ltot[5] = Ch_Ltot[4]+disk_gap+L_Guide43;
  Ch_Ltot[6] = Ch_Ltot[5]+disk_gap+L_Guide44;


  for (i=0;i<=6;i++)
  {
    Ch_phase[i]  =  Ch_Ltot[i]/v0;
    printf("Chopper %d: L=%lf [m] V=%lf [rad/s] Phase=%f [mu-sec]"
           "           op-times(%lf+n*%lf)s +/- %lf\n",
            i,Ch_Ltot[i],Ch_Vp[i], Ch_phase[i]*1.0e+6,
            Ch_phase[i],PI/Ch_Vp[i],2*PI/Ch_Vp[i]*Ch_alpha[i]/360.0);
  }

//========================================
//   Actual sample and detector
//========================================
  thickness    = 0.0125;
  radius    = 0.015;
  height     = 0.06;

  ang_ini = -11.9175;            //angular range of de detector in degrees
  ang_fin = 134.8172;            //
  det_angle = fabs(ang_fin-ang_ini)/2.0 + ang_ini;

%}



//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_//


TRACE

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

/*------------------*/
/*  SOURCE          */
/*------------------*/

REMOVABLE COMPONENT VCS = Source_gen(
  yheight  = 0.22,
  xwidth   = 0.14,
  focus_xw = 0.038,
  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 PREVIOUS

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

/*--------------------------------------------*/
/*    In the old version, here was the        */
/*    reset time for the choppers             */
/*    (t=0 for the choppers phase calculus)   */
/*                                            */
/*       Now,  t=0 is at Chopper0             */
/*--------------------------------------------*/
/////CHOPPER TIME-RESET/////////////////////////////////////////////////////
COMPONENT Chopper0  = DiskChopper(
  theta_0 = 20.222, radius = Ch_mean_R[0], yheight = 0.2,
  nu  = Ch_Vp[0]/2/PI, nslit = disk_N, delay = Ch_phase[0], isfirst=1)
AT (0,0,0.2) RELATIVE PREVIOUS
///////////////////////////////////////////////////////////////////////////

/*-------------------------*/
/*      GERADE Guide      */
/*-------------------------*/

COMPONENT Guide1  = Guide_channeled(
  w1 = 0.03000, h1 = 0.20000, w2 = 0.03000, h2 = 0.17415, l = L_Guide1,
  R0 = Guide_Ro, Qcx = Guide_Qc, Qcy = Guide_Qc,
  alphax = Guide_alpha, alphay = Guide_alpha,
  mx = 1, my=2 , W = Guide_W)
AT (0,0,L_gap) RELATIVE Chopper0

COMPONENT Guide21  = Guide_channeled(
  w1 = 0.03000, h1 = 0.17415, w2 = 0.03000, h2 = 0.17000, l = L_Guide21,
  R0 = Guide_Ro, Qcx = Guide_Qc, Qcy = Guide_Qc,
  alphax = Guide_alpha, alphay = Guide_alpha,
  mx = 1, my = 2 , W = Guide_W)
AT (0,0,L_Guide1+0.0003) RELATIVE Guide1


/*-------------------------------*/
/*     CHOPPER   - I -           */
/*-------------------------------*/

COMPONENT Chopper1  = DiskChopper(
  theta_0 = Ch_alpha[1], radius = Ch_mean_R[1], yheight = Ch_height[1],
  nu  = Ch_Vp[1]/2/PI, nslit = disk_N, delay = Ch_phase[1])
AT (0,0, L_Guide21+disk_gap/2) RELATIVE Guide21


//////GUIDE TO CHOPPER2//////////////////////////////////////////////
COMPONENT Guide22  = Guide_channeled(
  w1 = 0.03000, h1 = 0.17000, w2 = 0.03000, h2 = 0.16813, l = L_Guide22,
  R0 = Guide_Ro, Qcx = Guide_Qc, Qcy = Guide_Qc,
  alphax = Guide_alpha, alphay = Guide_alpha,
  mx = 2, my = 3 , W = Guide_W)
AT (0,0,L_Guide21+disk_gap) RELATIVE Guide21



/*-------------------------------*/
/*     CHOPPER   - II -          */
/*-------------------------------*/

COMPONENT Chopper2  = DiskChopper(
  theta_0 = Ch_alpha[2], radius = Ch_mean_R[2], yheight = Ch_height[2],
  nu   = Ch_Vp[2]/2/PI, nslit = disk_N, delay = Ch_phase[2])
AT (0,0, L_Guide22+disk_gap/2) RELATIVE Guide22


/*------------------*/
/*  MONOS - EINS-   */
/*                  */
/*  mid-way between */
/* Ch2 & L23        */
/*------------------*/
/*
COMPONENT M1 = Monitor_nD(xwidth=0.03, yheight=0.17,
  options="auto time")
AT (0,0, disk_gap/4+0.002) RELATIVE Chopper2
*/
/*-------------------------------*/
/*     Second guide part         */
/*                               */
/*-------------------------------*/

COMPONENT Guide23  = Guide_channeled(
  w1 = 0.03000, h1 = 0.16813, w2 = 0.02856, h2 = 0.15931, l = L_Guide23,
  R0 = Guide_Ro, Qcx = Guide_Qc, Qcy = Guide_Qc,
  alphax = Guide_alpha, alphay = Guide_alpha,
  mx = 2, my = 3 , W = Guide_W)
AT (0,0,L_Guide22+disk_gap) RELATIVE Guide22

COMPONENT Guide3  = Guide_channeled(
  w1 = 0.02856, h1 = 0.15931, w2 = 0.01733, h2 = 0.09041, l = L_Guide3,
  R0 = Guide_Ro, Qcx = Guide_Qc, Qcy = Guide_Qc,
  alphax = Guide_alpha, alphay = Guide_alpha,
  mx = 2, my = 3 , W = Guide_W)
AT (0,0,L_Guide23+0.0003) RELATIVE Guide23

COMPONENT Guide41  = Guide_channeled(
  w1 = 0.01733, h1 = 0.09041, w2 = 0.01579, h2 = 0.08100, l = L_Guide41,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha,
  Qcy = Guide_Qc, alphay = Guide_alpha, mx = 2, my = 3 , W = Guide_W)
AT (0,0,L_Guide3+0.0003) RELATIVE Guide3


/*-------------------------------*/
/*     CHOPPER   - III -         */
/*-------------------------------*/

COMPONENT Chopper3  = DiskChopper(
  theta_0 = Ch_alpha[3], radius = Ch_mean_R[3], yheight = Ch_height[3],
  nu   = Ch_Vp[3]/2/PI, nslit = disk_N, delay = Ch_phase[3])
AT (0,0, L_Guide41+disk_gap/2) RELATIVE Guide41

COMPONENT Guide42 = Guide_channeled(
  w1 = 0.01577, h1 = 0.08088, w2 = 0.01568, h2 = 0.08031, l = L_Guide42,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha,
  Qcy = Guide_Qc, alphay = Guide_alpha, mx = 2, my = 3, W = Guide_W)
AT (0,0,L_Guide41+disk_gap) RELATIVE Guide41

/*-------------------------------*/
/*     CHOPPER   - IV -          */
/*-------------------------------*/

COMPONENT Chopper4  = DiskChopper(
  theta_0 = Ch_alpha[4], radius = Ch_mean_R[4], yheight = Ch_height[4],
  nu   = Ch_Vp[4]/2/PI, nslit = disk_N, delay = Ch_phase[4])
AT (0,0, L_Guide42+disk_gap/2) RELATIVE Guide42

COMPONENT Guide43 = Guide_channeled(
  w1 = 0.01566, h1 = 0.08019, w2 = 0.01411, h2 = 0.07069, l = L_Guide43,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha, mx = 2,
  W  = Guide_W,  Qcy = Guide_Qc, alphay = Guide_alpha, my = 3)
AT (0,0,L_Guide42+disk_gap) RELATIVE Guide42

/*-------------------------------*/
/*     CHOPPER   - V -           */
/*-------------------------------*/

COMPONENT Chopper5  = DiskChopper(
  theta_0 = Ch_alpha[5], radius = Ch_mean_R[5], yheight = Ch_height[5],
  nu = Ch_Vp[5]/2/PI, nslit = disk_N, delay = Ch_phase[5])
AT (0,0,L_Guide43+disk_gap/2) RELATIVE Guide43

COMPONENT Guide44 = Guide_channeled(
  w1 = 0.01413, h1 = 0.07081, w2 = 0.01400, h2 = 0.0700, l = L_Guide44,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha, mx = 2,
  W  = Guide_W,  Qcy = Guide_Qc, alphay = Guide_alpha, my = 3)
AT (0,0,L_Guide43+disk_gap) RELATIVE Guide43


/*-------------------------------*/
/*     CHOPPER   - VI -           */
/*-------------------------------*/

COMPONENT Chopper6  = DiskChopper(
  theta_0 = Ch_alpha[6], radius = Ch_mean_R[6], yheight = Ch_height[6],
  nu = Ch_Vp[6]/2/PI, nslit = disk_N, delay = Ch_phase[6])
AT (0,0,L_Guide44+disk_gap/2) RELATIVE Guide44

COMPONENT Guide45 = Guide_channeled(
  w1 = 0.01400, h1 = 0.06983, w2 = 0.01400, h2 = 0.05663, l = L_Guide45,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha, mx = 2,
  W  = Guide_W,  Qcy = Guide_Qc, alphay = Guide_alpha, my = 3)
AT (0,0,L_Guide44+disk_gap) RELATIVE Guide44


/*//-----Here the actual 1st monitor*/

/*-----------------------*/
/*  KOLLIMATOR Guide    */
/*-----------------------*/
COMPONENT Collimator = Guide_channeled(
  w1 = 0.01400, h1 = 0.05617, w2 = 0.01400, h2 = 0.05400, l = L_Collimator,
  R0 = Guide_Ro, Qcx = Guide_Qc, alphax = Guide_alpha, mx = 2,
  W  = Guide_W,  Qcy = Guide_Qc, alphay = Guide_alpha, my = 3)
AT (0,0,L_Guide45+mono_gap) RELATIVE Guide45

/*----------------------------------------*/
/*              MONOS - ZWEI-             */
/*                                        */
/*             at sample position         */
/*                                        */
/*----------------------------------------*/

COMPONENT Det_sample_t = Monitor_nD(xwidth=0.014, yheight=0.054,
  options="auto t bins=20", restore_neutron=1)
AT (0,0,L_Collimator+0.0002) RELATIVE Collimator


/*-----------------------*/
/*       SAMPLE          */
/*-----------------------*/

COMPONENT arm2 = Arm()
AT (0,0,L_Collimator+L_CollSample+0.025)  RELATIVE Collimator
ROTATED (0,det_angle,0) RELATIVE Collimator

COMPONENT SAMPLE = Isotropic_Sqw(
  radius = radius, thickness=thickness, yheight = height,
  Sqw_coh=coh, Sqw_inc=inc, p_interact=0.9,
  order = order, d_phi = 180/PI*atan(1.5/4)*2, verbose=1)
AT (0,0,0) RELATIVE arm2
ROTATED (0,0,0) RELATIVE arm2
EXTEND
%{
   if(!SCATTERED) ABSORB;
%}

/*------------------------------*/
/*  DETECTOR AFTER SAMPLE       */
/*------------------------------*/

COMPONENT center_det = Arm()
AT (0,0,0) RELATIVE arm2
ROTATED (0,0,0) RELATIVE arm2

//--------------- DETECTOR IDEAL ----------------------------------------

COMPONENT Det_ideal_ay = Monitor_nD(xwidth=(4.0-0.0005-0.00002)*2, yheight=3,
  options="banana, theta limits=[-73.36735 73.36765] bins=100, y bins=100")
AT (0,0,0) RELATIVE center_det

//------------ Fe HOUSING------------------------------------------------

COMPONENT hous = PowderN(
  reflections=housing, radius = 4.0-0.00001, thickness = 0.0005,
  yheight = 3.0,p_transmit=0.8)
  //WHEN (housing && strcmp(housing,"0") && strcmp(housing,"NULL"))
AT (0,0,0) RELATIVE center_det
ROTATED (0,0,0) RELATIVE center_det

//------------ PSD Detector ---------------------------------------------

COMPONENT Det_PSD = PSD_Detector(
    yheight = 3.0, radius = 4.0, zdepth = 0.02600, awidth=(ang_fin-ang_ini)*PI/180*4.0,
    nx = 384, ny = 128, //type = "events",
    PressureConv = 4.75, PressureStop = 1.25, threshold=100,
    borderx=-1, bordery=-1, LensOn = 1, filename = "in5det.dat",
    FN_Conv="Gas_tables/He3inHe.table", FN_Stop="Gas_tables/He3inCF4.table")
AT (0,0,0) RELATIVE center_det
ROTATED (0,0,0) RELATIVE center_det

COMPONENT in5_t = Monitor_nD(
  options="banana, t limits=[0.0206 0.0216] bins=41, parallel, previous")
AT      (0,0,0) RELATIVE center_det
ROTATED (0,0,0) RELATIVE center_det

END