/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: FLEX (after full upgrade)
*
* %Identification
* Written by: M. Skoulatos and K. Habicht
* (up)Date: 19 October 2010
* Origin: Helmholtz-Zentrum Berlin
*
* Instrument short description
* Primary and secondary spectrometer for the FLEX upgrade
*
* %Parameters
* kI: [inv AA] incident wavevector
*
* %End
*******************************************************************************/

/* Change name of instrument and input parameters with default values */
DEFINE INSTRUMENT FLEX2112(kI=1.55, wVS=0.03, tilt=0, SA=-1, A4=50, L3=1.00, L4=1.00, rho_AH=1)

/* The DECLARE section allows us to declare variables or small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
#include <math.h>

// unit conversions
double m_n=1.674928e-27;     // in kg
double hbar=1.054572e-34;    // in J sec
double meV_to_SI=1.6022e-22; // conversion
double SI_to_meV=6.2414e21; 

// guide coating properties
double MGUIDE=3.2;
double W_para = 0.0025;
double R0_para = 0.99;
double Qc_para = 0.0217;
double alpha_para=3.90;

// curved guide properties
double LGUIDE_1=18.65;
double LGUIDE_2=17.8;
double beta_guide_1, s_guide_1;
double beta_guide_2, s_guide_2;

// monochromator properties
double RMH, RMV, rho_MH, rho_MV;
double lambdaI, EI;

// analyser properties
// double RAH, RAV, rho_AH, rho_AV;
double RAH;
double lambdaF, EF;

// angles
double A1, A2, A3, A5, A6;

// distances
double L1=1.75;
double L2=1.75;
// double L3=1.10;
// double L4=1.00;


//velocity selector frequency
double SelFreq;

// monitor/source properties
double dlambdaI,dkI,dEI,l_min,l_max,e_min,e_max;

%}

/* The INITIALIZE section is executed when the simulation starts     */
/* (C code). You may use them as component parameter values.         */
INITIALIZE
%{
#include <math.h>

int    SM=-1;      // monochromator scattering sense is always negative
int    SS=+1;
// int    SA=-1;
double DM = 3.355; // PG 002 monochromator d-spacing in inv AA   
double DA = 3.355; // PG 002 analyser d-spacing in inv AA

/* -------calculate spectrometer angles and monochromator/analyser curvatures-------*/    
// ----- monochromator
  A1 = SM*RAD2DEG*asin(PI/(DM*kI));   
  A2 = 2*A1;

//  RMH = L2/sin(A1*DEG2RAD);
  RMH = (L1+L2)/(2*sin(A1*DEG2RAD));
//  rho_MH=1/RMH;

//  RMV = 2*L2*sin(A1*DEG2RAD);
  RMV = (L1+L2)*sin(A1*DEG2RAD);
//  rho_MV=1/RMV;

  printf("monochromator angles are:\n A1=%7.3f, A2=%7.3f [deg]\n",A1,A2);
  printf("radius of curvature RMH=%7.3f [m]\nradius of curvature RMV=%7.3f [m] \n",RMH,RMV);

//  printf("curvature rhoMH=%7.3f [m^-1]\ncurvature rhoMV=%7.3f [m] \n",rho_MH,rho_MV);

// ----- sample
  A3 = 0;
//  A4 = 70;

printf("sample angles are:\n A3=%7.3f, A4=%7.3f [deg]\n",A3,A4);

// ----- analyser
  A5 = SA*RAD2DEG*asin(PI/(DA*kI));   
  A6 = 2*A5;
  RAH = 1/rho_AH;
//  need to put this in matlab:  RAH = (L3+L4)/(2*sin(A5*DEG2RAD));
//  RAV = (L3+L4)*sin(A5*DEG2RAD);
  printf("analyser angles are:\n A5=%7.3f, A6=%7.3f [deg]\n",A5,A6);
  printf("radius of curvature RAH=%7.3f [m] \n",RAH);
//  printf("radius of curvature RAV=%7.3f [m] \n",RAV);

/* ------------------------------calculate EI, lambdaI and so on ------------------*/    
  lambdaI=2*PI/kI;  // in AA
  EI=SI_to_meV*hbar*hbar*kI*kI*1e20/2/m_n;
  
  SelFreq = 3956.06/(lambdaI*(360/(19.7+(2.1*tilt)))*0.25);
//  empirical formulas for source emmitance to save time on simulations!
//  1st attempt: linear
//  dlambdaI=0.01*(7-lambdaI)*lambdaI;
//  2nd attempt, 1/x, better one:
//  dlambdaI=((0.082/(lambdaI+0.4))-0.002)*lambdaI;
//  3rd attempt, 6th order polynomial, even better (with enough parameters we can model an elephant):
 dlambdaI = (1.75004E-05*pow(lambdaI,6) - 5.13903E-04*pow(lambdaI,5) + 6.04404E-03*pow(lambdaI,4) - 3.66989E-02*pow(lambdaI,3) + 1.22399E-01*pow(lambdaI,2) - 2.19626E-01*lambdaI + 2.0E-01)*lambdaI;

  l_min=lambdaI-dlambdaI;
  l_max=lambdaI+dlambdaI;
  e_min=81.81/l_max/l_max;
  e_max=81.81/l_min/l_min;

  printf("velocity selector frequency is %7.3f [Hz]\n",SelFreq);  
  printf("lambda=%7.3f [A]\nkI=%7.3f [A^-1] \nEI=%7.3f [meV]\n",lambdaI,kI,EI);
  printf("Maxwellian source forced to emit only between l_min=%7.3f and l_max=%7.3f [A] in order to save simulation time\n",l_min,l_max);
  printf("correspondingly, the energy monitors at sample position are detecting between %7.3f and %7.3f [meV]\n",e_min,e_max);


// guide properties 
  beta_guide_1=LGUIDE_1/2800*180/PI;
  beta_guide_2=LGUIDE_2/2800*180/PI;
  s_guide_1 = 2800-sqrt(2800*2800-LGUIDE_1*LGUIDE_1);
  s_guide_2 = 2800-sqrt(2800*2800-LGUIDE_2*LGUIDE_2);
%}

/* Here comes the TRACE section, where the actual      */
/* instrument is defined as a sequence of components.  */
TRACE

/* The Arm() class component defines reference points and orientations  */
/* in 3D space. Every component instance must have a unique name. Here, */
/* Origin is used. This Arm() component is set to define the origin of  */
/* our global coordinate system (AT (0,0,0) ABSOLUTE). It may be used   */
/* for further RELATIVE reference, Other useful keywords are : ROTATED  */
/* EXTEND GROUP PREVIOUS. Also think about adding a neutron source !    */
/* Progress_bar is an Arm displaying simulation progress.               */

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

// note: keep source target window with fixed width 0.20 cm seperate from entrance window

COMPONENT Gen_Source = Source_gen(
    radius = 0.0775, dist = 1.53, xw = 0.2, yh = 0.125,
    Lmin = l_min, Lmax = l_max, I1 = 1e10,
    T1 = 45)
  AT (0, 0, 0) RELATIVE Origin
  ROTATED (0, 0, 0) RELATIVE Origin

// START complete neutron guide system in front of curved guide

// assumption: new extraction (in-pile) part is not rotated and not shifted w.r.t. NL2old axis
COMPONENT NL1A_NL1B_NL2_NL3_InPile_Entrance_Window = Arm()
  AT (0, 0, 1.530) RELATIVE Origin ROTATED (0, 0, 0) RELATIVE Origin

// this is the extraction part - with 4 mm symmetrically increased width and height
COMPONENT NL1A_NL1B_NL2_NL3_InPile = Guide(
    w1 = 0.170696, h1 = 0.129, w2 = 0.340522, h2 = 0.129, l = 1.870, 
    R0=R0_para, Qc=Qc_para, alpha=alpha_para, m=MGUIDE, W=W_para)
  AT (0, 0, 0) RELATIVE NL1A_NL1B_NL2_NL3_InPile_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1A_NL1B_NL2_NL3_InPile_Entrance_Window

// the NL1B_Straight1_Entrance_Window is shifted ans rotated according to T.Krists proposal 25.08.2008 - no gaps for simplicity
COMPONENT NL1B_Straight1_Entrance_Window = Arm()
  AT (-0.1277, 0, 3.400) RELATIVE Origin ROTATED (0, -2.15, 0) RELATIVE Origin  

// in the following straight guide sections the gaps are neglected and so only a single straight guide is considered!
COMPONENT NL1B_Straight1 = Guide(
    w1 = 0.06, h1 = 0.125, w2 = 0.06, h2 = 0.125, l = 1.549,
    R0=R0_para, Qc=Qc_para, alpha=alpha_para, m=MGUIDE, W=W_para)
  AT (0, 0, 0) RELATIVE NL1B_Straight1_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Straight1_Entrance_Window

COMPONENT NL1B_Curved_Entrance_Window = Arm()
  AT (0, 0, 1.549) RELATIVE NL1B_Straight1_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Straight1_Entrance_Window  

// END complete neutron guide system in front of curved guide
// START curved guide

COMPONENT NL1B_Curved_Guide_1 = Guide_curved(               
    w = 0.06, h = 0.125, l = LGUIDE_1, R0 = R0_para, Qc = Qc_para, 
    alpha = alpha_para, m = MGUIDE, W = W_para, r = 2800)
  AT (0, 0, 0) RELATIVE NL1B_Curved_Entrance_Window ROTATED (0, 0, 180) RELATIVE NL1B_Curved_Entrance_Window

//following arm translates and rotates to allow for the 18.65m guide. It also rotates by an extra 90 deg around z, to allow for a 12 o'clock velocity selector window.

COMPONENT NL1B_Velocity_Selector_Gap_Entrance_Window = Arm()
  AT (-s_guide_1, 0, LGUIDE_1) RELATIVE NL1B_Curved_Entrance_Window ROTATED (0, -beta_guide_1, 0) RELATIVE NL1B_Curved_Entrance_Window  

/************************************** 12h selector **************************************/

COMPONENT SELEC = Selector(
xmin=-0.0625,
xmax=0.0625,
ymin=-0.03,
ymax=0.03,
len=0.25,
num=72,
width=0.0004,
radius=0.123,
alfa=19.7,
feq=SelFreq
)

  AT (0,0,0.2) RELATIVE NL1B_Velocity_Selector_Gap_Entrance_Window 
  ROTATED (0,tilt,90)    RELATIVE NL1B_Velocity_Selector_Gap_Entrance_Window

COMPONENT NL1B_Curved_2_Entrance_Window = Arm()
  AT (0, 0, 0.45) RELATIVE NL1B_Velocity_Selector_Gap_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Velocity_Selector_Gap_Entrance_Window  

COMPONENT NL1B_Curved_Guide_2 = Guide_curved(               
    w = 0.06, h = 0.125, l = LGUIDE_2, R0 = R0_para, Qc = Qc_para, 
    alpha = alpha_para, m = MGUIDE, W = W_para, r = 2800)
  AT (0, 0, 0) RELATIVE NL1B_Curved_2_Entrance_Window ROTATED (0, 0, 180) RELATIVE NL1B_Curved_2_Entrance_Window

COMPONENT  NL1B_Straight2_Entrance_Window = Arm()
  AT (-s_guide_2, 0, LGUIDE_2) RELATIVE NL1B_Curved_2_Entrance_Window ROTATED (0,-beta_guide_2,0) RELATIVE NL1B_Curved_2_Entrance_Window

COMPONENT NL1B_Straight2 = Guide(
    w1 = 0.06, h1 = 0.125, w2 = 0.06, h2 = 0.125, l = 3.500, 
    R0=R0_para, Qc=Qc_para, alpha=alpha_para, m=MGUIDE, W=W_para)
  AT (0, 0, 0) RELATIVE NL1B_Straight2_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Straight2_Entrance_Window

// START elliptical guide section 

COMPONENT NL1B_Elliptical_Entrance_Window = Arm()
  AT (0, 0, 3.500) RELATIVE NL1B_Straight2_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Straight2_Entrance_Window  

COMPONENT elliptical_piece = Guide_tapering(
    option = "elliptical", segno = 800, w1 = 0.06, h1 = 0.125,
    l = 2.5, linw = 2.9, loutw = 0.4, R0 = 0.99, Qcx = 0.0217,
    Qcy = 0.0217, alphax = 4.00, alphay = alpha_para, W = 0.002,
    mx = 5.2, my = MGUIDE)
  AT (0, 0, 0) RELATIVE NL1B_Elliptical_Entrance_Window
  ROTATED (0, 0, 0) RELATIVE NL1B_Elliptical_Entrance_Window

// END elliptical guide section 

// XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
COMPONENT Virtual_source = Slit(
    width = wVS, height = 0.2)
  AT (0, 0, 2.900) RELATIVE NL1B_Elliptical_Entrance_Window
// XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */

// START GUIDE before monochromator 

COMPONENT NL1B_Vertical_Guide_Entrance_Window = Arm()
  AT (0, 0, 2.925) RELATIVE NL1B_Elliptical_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Elliptical_Entrance_Window

COMPONENT NL1B_Vertical_Guide = Guide_channeled(
    w1 = 0.03, h1 = 0.125, w2 = 0.15, h2 = 0.125, l = 1.475,
    R0 = R0_para, Qcx = Qc_para, Qcy = Qc_para,
    alphax = alpha_para, alphay = alpha_para, W = W_para, k = 1,
    mx = 0, my = MGUIDE)
  AT (0, 0, 0) RELATIVE NL1B_Vertical_Guide_Entrance_Window
  ROTATED (0, 0, 0) RELATIVE NL1B_Vertical_Guide_Entrance_Window

COMPONENT NL1B_Guide_Exit = Arm()
  AT (0, 0, 1.475) RELATIVE NL1B_Vertical_Guide_Entrance_Window ROTATED (0, 0, 0) RELATIVE NL1B_Vertical_Guide_Entrance_Window  

// COMPONENT psd_endguide = PSD_monitor(
//     nx = 200, ny = 200, filename = "psd-endguide.dat",
//     restore_neutron = 1, xwidth = 0.17, yheight = 0.135)
//   AT (0, 0, 0) RELATIVE NL1B_Guide_Exit
// 
COMPONENT energy_endguide = E_monitor(
    nchan = 200, filename = "energy-endguide.dat", restore_neutron = 1,
    xwidth = 0.17, yheight = 0.135, Emin = e_min, Emax = e_max)
  AT (0, 0, 0) RELATIVE NL1B_Guide_Exit

// COMPONENT hdiv_endguide = Hdiv_monitor(
//     nh = 200, filename = "hdiv-endguide.dat",
//     restore_neutron = 1, xwidth = 0.17, yheight = 0.135,
//     h_maxdiv = 3)
//   AT (0, 0, 0) RELATIVE NL1B_Guide_Exit
// 
// COMPONENT div_endguide = Divergence_monitor(
//     nh = 200, nv = 200, filename = "div-endguide.dat",
//     restore_neutron = 1, xwidth = 0.17, yheight = 0.135,
//     h_maxdiv = 3, v_maxdiv = 3)
//   AT (0, 0, 0) RELATIVE NL1B_Guide_Exit

// START monochromator part

COMPONENT Mono_center = Arm()
  AT (0, 0, 0.25) RELATIVE NL1B_Guide_Exit ROTATED (0,A1,0) RELATIVE NL1B_Guide_Exit

// set monochromator transmission to zero
COMPONENT Monochromator = Monochromator_curved(
    zwidth = 0.02,
    yheight = 0.02, 
    gap = 0.0005, 
    NH = 15,
    NV = 7, 
    mosaich = 40, 
    mosaicv = 40, 
    r0 = 1.0,
    t0 = 0.00001,
    RV = RMV, 
    RH = RMH
    )
  AT (0, 0, 0) RELATIVE Mono_center ROTATED (0,0,0) RELATIVE Mono_center

COMPONENT Mono_sample_arm = Arm()
  AT (0,0,0) RELATIVE Mono_center ROTATED (0,A2,0) RELATIVE NL1B_Guide_Exit

// COMPONENT psd_mono = PSD_monitor(
//     nx = 200, ny = 200, filename = "psd-mono.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.15)
//   AT (0, 0, 0.15) RELATIVE Mono_sample_arm
// 
COMPONENT energy_mono = E_monitor(
    nchan = 200, filename = "energy-mono.dat", restore_neutron = 1,
    xwidth = 0.2, yheight = 0.125, Emin = e_min, Emax = e_max)
  AT (0, 0, 0.15) RELATIVE Mono_sample_arm
// 
// COMPONENT hdiv_mono = Hdiv_monitor(
//     nh = 200, filename = "hdiv-mono.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.125,
//     h_maxdiv = 4)
//   AT (0, 0, 0.15) RELATIVE Mono_sample_arm
// 
// COMPONENT div_mono = Divergence_monitor(
//     nh = 200, nv = 200, filename = "div-mono.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.125,
//     h_maxdiv = 4, v_maxdiv = 4)
//   AT (0, 0, 0.15) RELATIVE Mono_sample_arm
// 
COMPONENT energy_pre_sample = E_monitor(
    nchan = 200, filename = "energy-pre-sample.dat", restore_neutron = 1,
    xwidth = 0.08, yheight = 0.08, Emin = e_min, Emax = e_max)
  AT (0, 0, L2-0.1) RELATIVE Mono_sample_arm


// END monochromator part


// Secondary spectrometer --->
//----------------
// start sample part

COMPONENT Sample_center = Arm()
  AT (0, 0, L2) RELATIVE Mono_sample_arm ROTATED (0,A3,0) RELATIVE Mono_sample_arm

COMPONENT vanabarba = V_sample(
    radius_i = 0.000, radius_o = 0.005, h = 0.08,
    focus_xw = 0.21, focus_yh = 0.135,
    target_index = +3)
// target should be "COMPONENT Analyser_center = Arm()"
  AT (0, 0, 0) RELATIVE Sample_center ROTATED (0,0,0) RELATIVE Sample_center
  
COMPONENT Sample_analyser_arm = Arm()
  AT (0,0,0) RELATIVE Sample_center ROTATED (0,A4,0) RELATIVE Mono_sample_arm

// COMPONENT psd_sam = PSD_monitor(
//     nx = 200, ny = 200, filename = "psd-sam.dat",
//     restore_neutron = 1, xwidth = 0.03, yheight = 0.1)
//   AT (0, 0, 0.1) RELATIVE Sample_analyser_arm
// 
COMPONENT energy_sam = E_monitor(
    nchan = 200, filename = "energy-sam.dat", restore_neutron = 1,
    xwidth = 0.03, yheight = 0.1, Emin = e_min, Emax = e_max)
  AT (0, 0, 0.11) RELATIVE Sample_analyser_arm
// 
// COMPONENT hdiv_sam = Hdiv_monitor(
//     nh = 200, filename = "hdiv-sam.dat",
//     restore_neutron = 1, xwidth = 0.03, yheight = 0.1,
//     h_maxdiv = 15)
//   AT (0, 0, 0.12) RELATIVE Sample_analyser_arm
// 
// COMPONENT div_sam = Divergence_monitor(
//     nh = 200, nv = 200, filename = "div-sam.dat",
//     restore_neutron = 1, xwidth = 0.03, yheight = 0.1,
//     h_maxdiv = 5, v_maxdiv = 5)
//   AT (0, 0, 0.13) RELATIVE Sample_analyser_arm

// end sample part
//----------------
// start analyser part

COMPONENT Analyser_center = Arm()
  AT (0, 0, L3) RELATIVE Sample_analyser_arm ROTATED (0,A5,0) RELATIVE Sample_analyser_arm

// set analyser transmission to zero
COMPONENT Analyser = Monochromator_curved(
    zwidth = 0.01287,
    yheight = 0.125, 
    gap = 0.0005, 
    NH = 15,
    NV = 1, 
    mosaich = 30, 
    mosaicv = 30, 
    r0 = 1.0,
    t0 = 0.00001,
    RV = 0, 
    RH = -RAH
    )
  AT (0, 0, 0) RELATIVE Analyser_center ROTATED (0,0,0) RELATIVE Analyser_center

COMPONENT Analyser_detector_arm = Arm()
  AT (0,0,0) RELATIVE Analyser_center ROTATED (0,A6,0) RELATIVE Sample_analyser_arm

// COMPONENT psd_ana = PSD_monitor(
//     nx = 200, ny = 200, filename = "psd-ana.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.2)
//   AT (0, 0, 0.1) RELATIVE Analyser_detector_arm
// 
COMPONENT energy_ana = E_monitor(
    nchan = 200, filename = "energy-ana.dat", restore_neutron = 1,
    xwidth = 0.2, yheight = 0.2, Emin = e_min, Emax = e_max)
  AT (0, 0, 0.1) RELATIVE Analyser_detector_arm
// 
// COMPONENT hdiv_ana = Hdiv_monitor(
//     nh = 200, filename = "hdiv-ana.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.2,
//     h_maxdiv = 5)
//   AT (0, 0, 0.1) RELATIVE Analyser_detector_arm
// 
// COMPONENT div_ana = Divergence_monitor(
//     nh = 200, nv = 200, filename = "div-ana.dat",
//     restore_neutron = 1, xwidth = 0.2, yheight = 0.2,
//     h_maxdiv = 5, v_maxdiv = 5)
//   AT (0, 0, 0.1) RELATIVE Analyser_detector_arm

// end analyser part

//final monitors at detector position

// COMPONENT psd_det = PSD_monitor(
//     nx = 200, ny = 200, filename = "psd-det.dat",
//     restore_neutron = 1, xwidth = 0.049, yheight = 0.150)
//   AT (0, 0, L4) RELATIVE Analyser_detector_arm

COMPONENT energy_det = E_monitor(
    nchan = 200, filename = "energy-det.dat", restore_neutron = 1,
    xwidth = 0.049, yheight = 0.150, Emin = e_min, Emax = e_max)
  AT (0, 0, L4) RELATIVE Analyser_detector_arm

// COMPONENT hdiv_det = Hdiv_monitor(
//     nh = 200, filename = "hdiv-det.dat",
//     restore_neutron = 1, xwidth = 0.049, yheight = 0.150,
//     h_maxdiv = 5)
//   AT (0, 0, L4) RELATIVE Analyser_detector_arm
// 
// COMPONENT div_det = Divergence_monitor(
//     nh = 200, nv = 200, filename = "div-det.dat",
//     restore_neutron = 1, xwidth = 0.049, yheight = 0.150,
//     h_maxdiv = 5, v_maxdiv = 5)
//   AT (0, 0, L4) RELATIVE Analyser_detector_arm

/* This section is executed when the simulation ends (C code). Other    */
/* optional sections are : SAVE                                         */
FINALLY
%{
%}
/* The END token marks the instrument definition end */
END