/**************************************HANARO-HRPD****************************
*%Description

*The HANARO-HRPD installed at ST2 beam tube which has a 1.28e14 (n/cm2.s) at beam nose.
*The Real (measured) flux at monochromator is 2.7e8 (n/cm2.s).
*The Real (measured) flux at Sample table is around 1.3e6 (n/cm2.s) for 1.836 A wavelength.

*The HANARO-HRPD consists four collimators ( inpile collimator + 3 outer collimators ) with (30', 10, 30', 10') divergence angles.

*The monochromator has a large mosaic spread of 20', two-fold type crystals with two take off angles (Theta_M) to give four wavelengths as following:
Ge 511 + Theta_M = 45  ---------> 1.539 A
Ge 511 + Theta_M = 59  ---------> 1.866 A
Ge 331 + Theta_M = 45  ---------> 1.836 A
Ge 331 + Theta_M = 59  ---------> 2.224 A
*The monochromator located at L0=6.06 m from the beam nose, where the sample table located at L1=3.0 m from the monochromator.

NOTE: this code is only to investigate the magnitude of flux (n/cm2.s) at both; Monochromator and Sample position.

NOTE2: in every section below you will find a short desription of components and their justified parameters as they are on ground*/

/*******************************************************************************************************************/

DEFINE INSTRUMENT HANARO_HRPD_1(I=1.28e14, Wlmin = 0.1 , Wlmax = 10 , Temp=37,
Alpha0=30, Alpha1=10, Alpha2=30,
 Beta=20, Theta_M=45, NV=9, NH=2, 
				L0=6.06, L1= 3.0, COLL_MOUNTED=1,GeMode=0, Narrow=0, RVinput=0, Reflectivity_M  = 0.301)


DECLARE
%{

  double DM, L, RV, Reflectivity_M, Lambda;
  double lminuse,lmaxuse;

%}


INITIALIZE
%{

  if (GeMode==0) {
    printf("Mono running in Ge(331) mode with these parameters\n");
    DM=1.29799;    
  } else if (GeMode==1) {
    printf("Mono running in Ge(511) mode with these parameters\n");
    DM=1.08885;
  } else {
    fprintf(stderr,"Unsupported monochromator mode!! Exiting!\n");
    exit(-1);
  }

  Lambda=2*DM*sin(Theta_M*DEG2RAD);
  printf("DM=%g, Reflectivity_M=%g, Theta_M=%g ==> Lambda=%g\n",DM,Reflectivity_M,Theta_M,Lambda);

  if (Narrow==0) {
    lminuse=Wlmin;
    lmaxuse=Wlmax;
  } else if (Narrow==1) {
    lminuse=0.99*Lambda;
    lmaxuse=1.01*Lambda;
  } else {
    fprintf(stderr,"Unsupported value of Narrow input parameter - only 0 or 1 possible!! Exiting!\n");
    exit(-1);
  }
  printf("Wavelength band produced at the source: [Lmin:Lmax]=[%g : %g ] AA\n",lminuse,lmaxuse);
  
  
/* ----calculation of the Monochromator's radius of curveture----- */
 L= 1 / (1/L0+1/L1);
 if (RVinput== 0) {
   RV=2*L*sin(Theta_M*DEG2RAD); 
 } else {
   RV=RVinput;
 }
 
 printf("Used radius of curvature is RV=%g\n",RV);
%}


TRACE

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

/* +++++++++++++++++++  Source +++++++++++++++++++*/

/* The source is a thermal maxwellian source with (0.07 x 0.12) m2 size focuses on the reactor wall opening of ( 0.056 x 0.12 ) m2 size at ( 3.54m ) distance */
/* NOTE: I use the  ( lambda0 = Lambda, dlambda = 0.01*Lambda )  below as a method  to examine the reflected wavelength and flux at sample position only,
 where I make it comment  when I'm trying to examine the flux at Monochromator position*/

COMPONENT S = Source_Maxwell_3(
    yheight = 0.12, xwidth = 0.07, Lmax = lmaxuse, Lmin = lminuse,
    //lambda0 = Lambda, dlambda = 0.01*Lambda,
    dist = 3.54 , focus_xw = 0.056, focus_yh = 0.12, T1 = Temp +273 ,                   
    I1 =I/4/PI)
  AT (0, 0, 0) RELATIVE origin


/* ++++++++++++++++ In-pile Collimator+++++++++++++++*/

/*NOTE: the in-pile collimator has the same dimension of the focusing window but it is located at 3.8 m from the beam nose with following characteristics*/

COMPONENT C0 = Collimator_linear(
    xmin = -0.028, xmax = 0.028, ymin = -0.06, ymax = 0.06,                           
    length = 0.40, divergence =Alpha0, transmission = 0.56)
  AT (0, 0,3.8) RELATIVE origin                   



/* +++++++++++++++ First Collimator +++++++++++++++++*/

/*NOTE: The ( +-0.001 and +-0.014 ) to get the real beam size (as measured in HANARO) at Monochromator */

COMPONENT C1 = Collimator_linear(
    xmin = -0.02+0.001, xmax = 0.02-0.001, ymin = -0.1+0.014, ymax = 0.1-0.014,                                   
    length = 0.350, divergence =Alpha1, transmission = 0.90)
  AT (0, 0, 4.916) RELATIVE origin                   


/* +++++++++++++++ Monochromator++++++++++++++++ */

/* The monochromator is vertically focusing type of 120mm in width and 212mm in height located at 6.06 m from beam nose, 
It is  (9 rows x 2 columns) of Ge crystals, each has (60mm in width  x 20mm in height) */

/* It has a large mosaic spread of 20', two-fold type crystals with two take off angles (Theta_M) to give four wavelengths as following:
Ge 511 + Theta_M = 45  ---------> 1.539 A 
Ge 511 + Theta_M = 59  ---------> 1.866 A
Ge 331 + Theta_M = 45  ---------> 1.836 A
Ge 331 + Theta_M = 59  ---------> 2.224 A */

/* The maximum reflectivity (r0) has been specified to be (r0= 0.301 for "Ge 331", and  r0= 0.218 for "Ge 511" ) i.e (1.539 A and 1.834 A both have the same r0)*/

/*NOTE: eventhough  I believe that r0 should be different for each wavelength but I inserted it as I they defined it in the HANARO_instrumentations documents */

/* NOTE2: I used a psd_monitor with a window size of (1cm x 1cm) in order to get the intensity measured by McStas in flux unit (n/cm2.s) */

COMPONENT Before_Mono = PSD_monitor(
				    filename = "Before_Mono1.psd" ,  xwidth = 0.01, yheight = 0.01, restore_neutron=1)
  AT (0, 0, L0) RELATIVE origin

  COMPONENT Before_Mono_Capture = Monitor_nD(filename = "Before_Mono_Capture.psd" ,  xwidth = 0.01, yheight = 0.01, options="capture, x bins=10, y bins=10", restore_neutron=1)
  AT (0, 0, L0) RELATIVE origin


COMPONENT Mono_Arm_1 = Arm(
    )
  AT (0, 0, L0) RELATIVE origin    

SPLIT 50 COMPONENT Mono = Monochromator_curved(
    width = 0.12, height = 0.212,                         
    zwidth = 0.06, yheight= 0.02,
    NH =NH, NV =NV,
    DM = DM, mosaich =Beta , mosaicv = Beta , RV = RV, verbose=1,
    order= 1,
    r0 = Reflectivity_M)
  AT (0, 0, 0) RELATIVE Mono_Arm_1
  ROTATED (0, Theta_M, 0) RELATIVE Mono_Arm_1 
EXTEND %{
  /* Remove the "direct beam", only transport reflected neutrons further */
  if (!SCATTERED) ABSORB;
%}



  COMPONENT Mono_Arm_2 = Arm()
  AT (0, 0, 0) RELATIVE Mono_Arm_1 
  ROTATED (0,2*Theta_M, 0) RELATIVE Mono_Arm_1 

COMPONENT After_Mono = PSD_monitor(
    filename = "After_Mono1.psd", restore_neutron = 1, xwidth = 0.01, yheight = 0.01)
  AT (0, 0, 1) RELATIVE  Mono_Arm_2


/*+++++++++++++++Second Collimator++++++++++++++++++*/

 /* NOTE: The (+-0.003 and -+0.014 ) to get the real beam size (as measured in HANARO) at sample */

COMPONENT C2 = Collimator_linear(
    xmin = -0.0175, xmax = 0.0175,                           
    ymin = -0.075, ymax = 0.075,                               
    length = 0.350, divergence =Alpha2, transmission = 0.9)
WHEN (COLL_MOUNTED)  AT (0, 0, 0.865) RELATIVE  Mono_Arm_2



/*+++++++++++++++ Sample Table ++++++++++++++ */

/* The sample table located at L1=3.0 m from the monochromator.*/
/* NOTE: I used a psd_monitor with a window size of (1cm x 1cm) in order to get the intensity measured by McStas in flux unit (n/cm2.s) */

COMPONENT Sample_Arm = Arm(
    )
  AT (0, 0, L1) RELATIVE  Mono_Arm_2

COMPONENT Before_Sa_BIG = PSD_monitor(
    filename = "Before_Sa1_BIG.psd", restore_neutron = 1, xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0) RELATIVE  PREVIOUS    

COMPONENT Before_Sa = PSD_monitor(
    filename = "Before_Sa1.psd", restore_neutron = 1, xwidth = 0.01, yheight = 0.01)
  AT (0, 0, 0) RELATIVE  PREVIOUS     

 COMPONENT Before_Sa_capture = Monitor_nD(options="capture, x bins=10, y bins=10", 
    filename = "Before_Sa1_capture.psd", restore_neutron = 1, xwidth = 0.01, yheight = 0.01)
  AT (0, 0, 0) RELATIVE  PREVIOUS                                                        

COMPONENT Before_Sa_L = L_monitor(
    filename = "Before_Sa1.L", Lmin = 0, nL=100,
    Lmax = Lambda*1.02, restore_neutron = 1, xwidth = 0.01, yheight = 0.01)
  AT (0, 0, 0) RELATIVE PREVIOUS                                

COMPONENT Before_Sa_div = Divergence_monitor(
					     filename = "Before_Sa1.Div", maxdiv_v=2, maxdiv_h=2, nh=100, nv=100, restore_neutron = 1, xwidth = 0.01, yheight = 0.01)
  AT (0, 0, 0) RELATIVE PREVIOUS                                


END