/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Reflectometer
*
* %Identification
* Written by: Anette Vickery, contact: anette.vickery@fys.ku.dk
* Date: November 2011
* Origin: KU
* %INSTRUMENT_SITE: Templates
*
* Horizontal reflectometer, multi-angle of incidence
*
* %Description
* A horizontal reflectometer. The instrument is built of a 2.5 m long mirror and an inclined elliptical guide. The grazing angle is defined by a set of
* slits. The sample size is 4cm x 4cm (horizontal sample).
* The role of the mirror is to provide an angle of incidence of up to 4 deg and at the same time avoid a direct line of sight to the source.
* The reflectivity is the ratio between the direct beam reflected beam intensities. Therefore two simulations are needed to get the reflectivity curve:
*
* Example: mcrun Reflectometer.instr -n1e10 directbeam=1,thetasample=0.4,Qmin=0,Qmax=0.15 -d directbeam
* Example: mcrun Reflectometer.instr -n1e10 directbeam=0,thetasample=0.4,Qmin=0,Qmax=0.15 -d reflectedbeam
*
* %Parameters
* directbeam: []        If 1, the sample is a mirror with reflectivity 1. If 0, the mirror reflects like a D2O surface with zero roughness
* Lam_min: [AA]        Minimum wavelength emitted by the cold moderator
* Lam_max: [AA]        Maximum wavelength emitted by the cold moderator
* deltatheta: [%]      The angular resolution
* theta_sample: [deg]  The grazing angle
* defineAngle: []       If 0, the angle defining slits are wide open. If 1 the angular resolution is deltatheta
* Qmin: [AA^-1]        The minimum value of the interesting Qrange
* Qmax: [AA^-1]        The maximum value of the interesting Qrange
* straight: []          If 1, the guide walls are straight
* width: [m]           The width of the guide in the case of straight guide walls
* directbeam: []       
* defineAngle: []      
* straight: []         
*******************************************************************************/
DEFINE INSTRUMENT Reflectometer( directbeam=1, Lam_min=2.0, Lam_max=10, deltatheta=2, theta_sample=1, defineAngle=1, Qmin=0, Qmax=0.5, straight=1, width=0.05 )
DECLARE
%{
double Y_sample=-0.000516, Y_mirror=-6e-3, Y_guide=-0.008385, MINDIV=-4.5,MAXDIV=-0.1;
double focus_inv=1.581542, focus_outv=1.957539, small_axis_h=0.25, T_mirror=1.136359, TT_guide=2.140996;
double focus_inh=1.582,focus_outh=1.958,small_axis_w=0.25;
double deltaLambda=1e-3;
double foot=40,fxw=0.13,fyh=0.1;
double sample_Z=0.0; // negative sample_Z: move monitor upstream
double zlos=0.59,cutoffguide=0,slitGravity=1,slitGravCor=7.2075 ;// the slit height are adjusted due to gravity. Mean displacements are as for 7.21AA neutrons
double vCor,deltaY_Grav1,deltaY_Grav2; // slitheight correction due to gravity.
double h[2], w[2];
double focus_s_w; double focus_e_w;
double focus_s_h; double focus_e_h;
double elength_h; double elength_w;
// parameters for guide calculation
double guide_start=1.27+1e-6;
double W_par=0.003, R0_par=0.99,Qcrit=0.0217,mvalue=6,alpha_par=6.07;

// parameters for slit settings
 double S1gap=2,S2gap=2;
double Y_Slit1,Y_Slit2;// relative to guideaxis

// parameters for the source
double Pulse_width=0.00286;
double frequency=14;
double guide_dist, Pulse_freq;

// Parametes for the mirror
double L_mirror=2.5;         //length
double h_mirror=0.12;        //width

double dist_moderator_detector=30, dist_moderator_mirrorbegin=2, dist_guide_S1=0.02, dist_S1_S2=1.5, dist_sample_detector=3;
double guide_length,dist_guide_S2,dist_guide_sample,dist_S2_sample;

double sampleposition=27-3.25;// distance mirror-arm--sample
// Parameters for defining source focusing:
double focus_yh; // height of focusing rectangle
char optionstring1 [128];
char optionstring2 [128];
char optionstring3 [128];
char optionstring4 [128];
char optionstring5 [128];
char optionstring51 [128];
char optionstring6 [128];
char optionstring6x [128];
double deltaTheta;
double LOSSlitheight;
%}

USERVARS %{
  double mirflag;
  double guide1flag;
  double guide2flag;
  double reflectflag;
%}

INITIALIZE
%{

if (!Lam_max){
    Lam_max=Lam_min+deltaLambda;
}

 if (straight){
   focus_inh=1000;
   focus_outh=1000;
   small_axis_w=width;
 }
printf("-------sampleposition= %g m \n", sampleposition);
guide_length=30-(2+2.5+2+3);
printf("--------------------------guide_length= %g m \n", guide_length);
/* -------------------- square guide cross section: -------------------- */
/*  naming convention: focus_inh is the horizontal focus
    and focus_inv is the vertical focus. Distance from guide to focalpoints.
    The horizontal corrosponds to the width and vertical the height. */
    printf("-------------------------- focus_outh = %g m \n",focus_outh);


    focus_s_h=-focus_inv;
    focus_e_h=guide_length+focus_outv;

    focus_s_w=-focus_inh;
    focus_e_w=guide_length+focus_outh;

    printf("-------------------------- focus_s_h = %g m \n",focus_s_h);
    printf("-------------------------- focus_e_h = %g m \n",focus_e_h);

    elength_h=focus_e_h-focus_s_h;
    elength_w=focus_e_w-focus_s_w;

    h[0]=small_axis_w*sqrt(1-(((0-focus_s_h)-elength_h/2)/(elength_h/2))*(((0-focus_s_h)-elength_h/2)/(elength_h/2)));

    h[1]=small_axis_w*sqrt(1-(((1*guide_length-focus_s_h)-elength_h/2)/(elength_h/2))*(((1*guide_length-focus_s_h)-elength_h/2)/(elength_h/2)));

    w[0]=small_axis_h*sqrt(1-(((0-focus_s_w)-elength_w/2)/(elength_w/2))*(((0-focus_s_w)-elength_w/2)/(elength_w/2)));

    w[1]=small_axis_h*sqrt(1-(((1*guide_length-focus_s_w)-elength_w/2)/(elength_w/2))*(((1*guide_length-focus_s_w)-elength_w/2)/(elength_w/2)));

    printf("-------------------------- h0, h1 = %g, %g m \n",h[0],h[1]);
    printf("-------------------------- w0, w1 = %g, %g m \n",w[0],w[1]);

dist_S2_sample=sampleposition-(guide_start+guide_length+dist_S1_S2+dist_guide_S1);
deltaTheta= deltatheta*DEG2RAD*theta_sample/100;

if (defineAngle){
   S2gap= foot*sin(DEG2RAD*theta_sample)*1e-3; // slitheight to illuminate footprint only
   S1gap = sqrt((deltaTheta*dist_S1_S2*1e3/0.68)*(deltaTheta*dist_S1_S2*1e3/0.68)-S2gap*S2gap)*1e-3;
   }
printf("--------------------------deltaTheta= %g rad \n",deltaTheta);

printf("--------------------------dist_guide_S1= %g mm \n",dist_guide_S1*1e3);
printf("--------------------------dist_S1_S2= %g mm \n",dist_S1_S2*1e3);
printf("--------------------------dist_S2_sample= %g mm \n",dist_S2_sample*1e3);

printf("--------------------------sampleposition= %g mm \n",sampleposition*1e3);
printf("--------------------------S1gap= %g mm \n",S1gap*1e3);
printf("--------------------------S2gap= %g mm \n",S2gap*1e3);


//position of slits not taking gravity into account:
Y_Slit1=-((dist_S1_S2+dist_S2_sample)*(sin(TT_guide*PI/180)-cos(TT_guide*PI/180)*tan(theta_sample*PI/180))-Y_sample);// relative to guide_axis
Y_Slit2=-((dist_S2_sample)*(sin(TT_guide*PI/180)-cos(TT_guide*PI/180)*tan(theta_sample*PI/180))-Y_sample);// relative to guide_axis
printf("--------------------------Y_Slit1= %g m \n", Y_Slit1);
printf("--------------------------Y_Slit2= %g m \n", Y_Slit2);
vCor=3.956e3/slitGravCor; // m/s velocity for slitGravCor wavelength neutrons

//ys1 = g*(zs1^2)./(2.*(v(i).*cos(theta(h))).^2) + tan(theta(h))*zs1;
deltaY_Grav1= -9.81*(dist_S1_S2+dist_S2_sample)*cos(TT_guide*PI/180)*(dist_S1_S2+dist_S2_sample)*cos(TT_guide*PI/180)/(2*vCor*vCor*cos(theta_sample*PI/180)*cos(theta_sample*PI/180) );
deltaY_Grav2= -9.81*(dist_S2_sample)*cos(TT_guide*PI/180)*(dist_S2_sample)*cos(TT_guide*PI/180)/(2*vCor*vCor*cos(theta_sample*PI/180)*cos(theta_sample*PI/180) );
if (slitGravity){
 Y_Slit1=Y_Slit1+deltaY_Grav1; // updating to take gravity into account
 Y_Slit2=Y_Slit2+deltaY_Grav2;// updating to take gravity into account
}

printf("--------------------------deltaY_Grav1 = %g m \n", deltaY_Grav1*1e6);
printf("--------------------------deltaY_Grav2 = %g m \n", deltaY_Grav2*1e6);
printf("--------------------------Y_Slit1= %g m \n", Y_Slit1);
printf("--------------------------Y_Slit2= %g m \n", Y_Slit2);

focus_yh = L_mirror*sin(T_mirror*PI/180); // height of focusing rectangle
printf("--------------------------focus_yh= %g m \n", focus_yh);


printf("--------------------------focus_inh= %g m \n", focus_inh);
printf("--------------------------focus_outh= %g m \n", focus_outh);
printf("--------------------------focus_inv= %g m \n", focus_inv);
printf("--------------------------focus_outv= %g m \n", focus_outv);
printf("--------------------------guide_length= %g m \n",guide_length);
printf("--------------------------guide_start= %g m \n",guide_start);
printf("--------------------------widthOfEllipsev= %g m \n",small_axis_h*0.5);
printf("--------------------------widthOfEllipseh= %g m \n",small_axis_w*0.5);
printf("--------------------------R0=%g m \n",R0_par);
printf("--------------------------Qc=%g m \n",Qcrit);
printf("--------------------------alpha=%g m \n",alpha_par);
printf("--------------------------m=%g m \n",mvalue);
printf("--------------------------W=%g m \n",W_par);

%}

TRACE

COMPONENT Origin = Progress_bar()
  AT (0,0,0) ABSOLUTE
EXTEND %{
guide1flag=0;
guide2flag=0;
reflectflag=0;
mirflag=0;
%}

COMPONENT cold_source = ESS_butterfly(
    Lmin = Lam_min, Lmax = Lam_max, dist = 2,
    focus_xw = fxw, focus_yh =fyh)
  AT (0, 0, 0) RELATIVE Origin


//Position of mirror
  COMPONENT mirror_pos1=Arm()
  AT (0,Y_mirror+(1.25-zlos*0.5)*tan(T_mirror*PI/180), 2+zlos*0.5) RELATIVE Origin
ROTATED (90+T_mirror,0,0) RELATIVE Origin

COMPONENT mirror1 = Mirror(
    xwidth = h_mirror, yheight = zlos, center = 1, m=6, alpha=3.5,
    transmit = 0)
  AT (0, 0, 0) RELATIVE mirror_pos1
EXTEND
   %{
   if (SCATTERED) {
      mirflag=1;
      reflectflag=1;
    }
    %}

COMPONENT mirror_pos0=Arm()
AT (0,Y_mirror, 3.25) RELATIVE Origin
ROTATED (90+T_mirror,0,0) RELATIVE Origin
COMPONENT mirror0 = Mirror(
			   xwidth = h_mirror, yheight = 2*(1.25-zlos), center = 1, m=6, alpha=3.5,
    transmit = 0)
  AT (0, 0, 0) RELATIVE mirror_pos0
EXTEND
   %{
   if (SCATTERED) {
      mirflag=1;
      reflectflag=1;
    }
    %}

  COMPONENT mirror_pos2=Arm()
  AT (0,Y_mirror-(1.25-zlos*0.5)*tan(T_mirror*PI/180), 4.5-zlos*0.5) RELATIVE Origin
ROTATED (90+T_mirror,0,0) RELATIVE Origin

COMPONENT mirror2 = Mirror(
    xwidth = h_mirror, yheight = zlos, center = 1, m=6, alpha=3.5,
    transmit = 0)
  AT (0, 0, 0) RELATIVE mirror_pos2
EXTEND
   %{
   if (SCATTERED) {
      mirflag=1;
      reflectflag=1;
    }
    %}

COMPONENT guide_TT=Arm()
AT (0,Y_guide, 3.25) RELATIVE Origin
ROTATED (TT_guide,0,0) RELATIVE Origin


COMPONENT guide1 = Elliptic_guide_gravity(
    l=guide_length*0.5-1e-6,
    linxw = focus_inh,
    loutxw = focus_outh+guide_length*0.5,
    linyh = focus_inv,
    loutyh= focus_outv+guide_length*0.5,
    minorAxisxw=small_axis_h*0.5,
    minorAxisyh=small_axis_w*0.5,
    R0=R0_par,
    Qc=Qcrit,
    alpha=alpha_par,
    m=mvalue,
    W=W_par)
AT (0, 0,guide_start) RELATIVE guide_TT
EXTEND
   %{
   if (SCATTERED) {
      guide1flag=1;
      reflectflag=1;
   }
    %}


COMPONENT guide2 = Elliptic_guide_gravity(
    l=guide_length*0.5-1e-6-cutoffguide,
    linxw = focus_inh+guide_length*0.5,
    loutxw= focus_outh+cutoffguide,
    linyh = focus_inv+guide_length*0.5,
    loutyh= focus_outv+cutoffguide,
    minorAxisxw=small_axis_h*0.5,
    minorAxisyh=small_axis_w*0.5,
    R0=R0_par,
    Qc=Qcrit,
    alpha=alpha_par,
    m=mvalue,
    W=W_par)
AT (0, 0,guide_start+guide_length*0.5) RELATIVE guide_TT
EXTEND
   %{
   if (SCATTERED) {
      reflectflag=1;
      guide2flag=1;
    }
    %}

/* Insert vertical slits:*/
COMPONENT Slit1Arm=Arm()
  AT (0, 0, guide_start+guide_length+dist_guide_S1) RELATIVE guide_TT
ROTATED (0,0,0) RELATIVE ABSOLUTE

COMPONENT Slit11Arm=Arm()
  AT (0, Y_Slit1, 0) RELATIVE Slit1Arm

COMPONENT Slit1 = Slit(
    xwidth = 0.2, yheight=S1gap)
    AT (0,0,1e-5) RELATIVE Slit11Arm

COMPONENT Slit2Arm=Arm()
  AT (0, 0, guide_start+guide_length+dist_guide_S1+dist_S1_S2) RELATIVE guide_TT
ROTATED (0,0,0) RELATIVE ABSOLUTE

COMPONENT Slit22Arm=Arm()
  AT (0, Y_Slit2, 0) RELATIVE Slit2Arm

COMPONENT Slit2 = Slit(
    xwidth = 1, yheight=S2gap)
    AT (0,0,1e-5) RELATIVE Slit22Arm

//----------------position sample :
COMPONENT Hor_sampleArm=Arm()
  AT (0, 0, sampleposition) RELATIVE guide_TT
ROTATED (-(TT_guide),0,0) RELATIVE guide_TT

COMPONENT Ver_sampleArm=Arm()
  AT  (0,Y_sample,0) RELATIVE Hor_sampleArm
ROTATED (90,0,0) RELATIVE Hor_sampleArm

COMPONENT mirrorDirect = Mirror(
    reflect = "Reflectometer_directbeam.txt", xwidth = 4e-2, yheight = 4e-2,center=1,
    transmit = 0) WHEN (directbeam)
  AT (0, sample_Z, 0) RELATIVE Ver_sampleArm
ROTATED (0,0,90) RELATIVE Ver_sampleArm
EXTEND
   %{
   if (!SCATTERED) {
	ABSORB;
    }
    %}

COMPONENT mirror = Mirror(
    reflect = "Reflectometer_reffile.txt", xwidth = 4e-2, yheight = 4e-2,center=1,
    transmit = 0) WHEN (!directbeam)
  AT (0, sample_Z, 0) RELATIVE Ver_sampleArm
ROTATED (0,0,90) RELATIVE Ver_sampleArm
EXTEND
   %{
   if (!SCATTERED) {
	ABSORB;
    }
    %}

COMPONENT TOF2QcylPSD = TOF2Q_cylPSD_monitor(
    nQ = 100, filename = "TOF2QcylPSD_1", ny = 1, radius = 3,
    yheight = 1, Qmin=Qmin, Qmax=Qmax, T_zero=0.5*2.86*1e-3, L_flight=30, restore_neutron = 1,theta=DEG2RAD*theta_sample)
  AT (0, 0, 0) RELATIVE  Hor_sampleArm
ROTATED (0,0,90) RELATIVE ABSOLUTE

FINALLY
%{
%}
END