/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Template monochromator Diffractometer
*
* %Identification
* Written by: C. M. I Enrique, K. Lieutenant, E. Farhi and L. Cussen
* Date: 13 Apr 2006
* Origin: ILL
* %INSTRUMENT_SITE: ILL
*
* Simple monochromator Diffractometer for powders
*
* %Description
* The diffractometer D2B is characterised by the very high take-off angle (135 deg)
* for the monochromator, which has a relatively large mosaic spread of 20' to
* compensate for the corresponding intensity (dl/l) loss. It is 300 mm high,
* focusing vertically onto about 50 mm; this large incident vertical divergence
* is matched by 200 mm high detectors and collimators. A complete diffraction
* pattern is obtained after about 100 steps of 0.025 deg in 2theta, since the 64
* detectors are spaced at 2.5 deg intervals. Such scans take typically 30 minutes;
* they are repeated to improve statistics.
*
* D2B was designed for work on samples and high resolution of very large
* d-spacings using wavelengths of between 2.4 Å and 6 Å. Wavelengths can easily be
* changed under computer control, since they are all obtained by a simple rotation
* within the Ge[hhl] plane. A large graphite filter can be switched in to provide
* a very clean beam at 2.4 Angs, and a cold Be-filter can be used for longer
* wavelengths.
*
* This model implements as well the Caglioti UVW equations, that give estimates
* of the instrument resolution.
*
* Monochromator lattice parameter
* Ge       111 DM=3.266 AA
* Ge       311 DM=1.714 AA
* Ge       511 DM=1.089 AA
* Ge       533 DM=0.863 AA
*
* %Example: lambda=1 Detector: D2B_BananaPSD_I=1.4489E+04
*
* %Parameters
* lambda: [Angs]  Wavelength at monochromator, computed from DM and THETA_M if left as 0.
* coating: [1]    Super-mirror in-beam tube guide coating
* DM: [Angs]      d-spacing of monochromator, computed from lambda and THETA_M if left as 0.
* THETA_M: [deg]  Monochromator take-off angle, computed from lambda and DM if left as 0.
* TILT: [deg]     Monochromator additional tilt, for rocking curves
* RV: [m]         Monochromator vertical curvature, 0 for flat, -1 for automatic setting
* L1: [m]         Source-Monochromator distance
* L2: [m]         Monochromator-Sample distance
* L3: [m]         Sample-Detector distance
* Powder: [str]   File name for powder description
* verbose: [1]    Print DIF configuration. 0 to be quiet
* ALPHA1: [min]   Horizontal collimator divergence for L1 arm (before monochromator)
* ALPHA2: [min]   Horizontal collimator divergence for L2 arm (monochromator-sample)
* ALPHA3: [min]   Horizontal collimator divergence for L3 arm (sample-detector)
* ETA: [min]      Monochromator horizontal mosaic (gaussian)
* SM: [1]         Scattering sense of beam from Monochromator. 1:left, -1:right
* Dheight: [m]    Banana detector height
*
* %Link
* G. Caglioti, A. Paoletti, F.P. Ricci, Nucl. Instr. and. Meth. 3 (1958) 223
* %Link
* L.D. Cussen, Nucl. Instr. and. Meth. A 554 (2005) 406
* %Link
* M. Morhac, NIM A 600 (2009) 478
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT ILL_D2B(lambda=1.594, DM=0, string Powder="Na2Ca3Al2F14.laz",
	RV=-1, L1=16.05, L2=2.645, L3=1.3, verbose=1, ALPHA1=18,ALPHA2=11,ALPHA3=5,
	ETA=12, THETA_M=67.5, TILT=0, SM=1, Dheight=0.3,coating=0)

/* The DECLARE section allows us to declare variables or  small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
	double U,V,W;   /* Caglioti UVW parameters */
	double L0;
	char lamopts[128];
%}

/* The INITIALIZE section is executed when the simulation starts     */
/* (C code). You may use them as component parameter values.         */
INITIALIZE
%{
	double L;
	double KI, Vi, EI;
	double ALPHA3_Opt;

	if (!THETA_M && lambda && DM)
		THETA_M =asin(lambda/(2*DM))*RAD2DEG;
	else if (THETA_M && !lambda && DM)
		lambda = fabs(sin(THETA_M*DEG2RAD))*2*DM;
	else if (THETA_M && lambda)
		DM = fabs(lambda/sin(DEG2RAD*THETA_M)/2.0);

	THETA_M *= SM; /* take-off direction left or right */

	/* test input parameters */
	if (!THETA_M || !DM || !lambda) exit(fprintf(stderr,
	      "ILL_D2B: ERROR: Monochromator take-off, d-spacing or wavelength is null (THETA_M=%g, DM=%g, lambda=%g). Abort.\n",
	      THETA_M, DM, lambda));

	/* perform additional computations */
	if (!L1)    L = L2;
	else        L = 1/(1/L1+1/L2);
	if (RV < 0) RV=2*L*sin(DEG2RAD*THETA_M);
	KI = 2*PI/lambda;
	Vi = K2V*fabs(KI);
	EI = VS2E*Vi*Vi;

	if (ALPHA1<0)  ALPHA1=600;
	if (ALPHA2<0)  ALPHA2=600;
	if (ALPHA3<0)  ALPHA3=600;
	if (ETA<=0)    ETA   =30;

	/* compute Cagioloti U,V,W and L0 (intensity) */
	U = 4*(ALPHA1*ALPHA1 * ALPHA2*ALPHA2 + ALPHA1*ALPHA1 * ETA*ETA + ALPHA2*ALPHA2 * ETA*ETA)
	     /(tan(THETA_M*DEG2RAD)*tan(THETA_M*DEG2RAD) * (ALPHA1*ALPHA1 + ALPHA2*ALPHA2 + 4*ETA*ETA));

	V = 4*ALPHA2*ALPHA2*(ALPHA1*ALPHA1 +2*ETA*ETA)
	     /(tan(THETA_M*DEG2RAD) * (ALPHA1*ALPHA1 + ALPHA2*ALPHA2 + 4*ETA*ETA));

	W = ALPHA3*ALPHA3+(ALPHA1*ALPHA1 * ALPHA2*ALPHA2 + 4*ALPHA2*ALPHA2 * ETA*ETA)
	                 /(ALPHA1*ALPHA1 + ALPHA2*ALPHA2 + 4*ETA*ETA);
	L0 = (ALPHA1*ETA*ALPHA2*ALPHA3/tan(fabs(THETA_M)*DEG2RAD))
	   / sqrt(ALPHA1*ALPHA1 + ALPHA2*ALPHA2 + 4*ETA*ETA);

	ALPHA3_Opt = 1/sqrt(1/ALPHA1/ALPHA1+1/ALPHA2/ALPHA2+1/ETA/ETA);

	if (verbose) {
	  printf("ILL_D2B: Detailed DIF configuration\n");
	  printf("* Incoming beam: lambda=%.4g [Angs] EI=%.4g [meV]  KI=%.4g [Angs-1] Vi=%g [m/s]\n",
	    lambda, EI, KI, Vi);
	  printf("* Monochromator:  DM=%.4g [Angs] RV=%.4g [m] %s, take-off 2*THETA_M=%.4g [deg]\n",
	    DM, RV, (!RV ? "flat" : "curved"), THETA_M*2);
	  printf("* Resolution:     Caglioti              U=%g V=%g W=%g [deg^2]\n", U/60/60, V/60/60, W/60/60);
	  printf("                  Optimal          ALPHA3=%g [min]\n",  ALPHA3_Opt);
	  printf("                  Optimal resolution A1/2=%g [min]\n",  ALPHA3_Opt/sqrt(2));
	  printf("                  Normalized Intensity L0=%g [min^3]\n", L0);
	}

	RV = fabs(RV)*SM;
	sprintf(lamopts, "lambda limits=[%g, %g]",lambda*0.99, lambda*1.01);
%}

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

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

/* source with constant flux. Directly send neutron events to monochromator area */
COMPONENT Source = Source_gen(dist = 11, focus_xw = fabs(0.05*fabs(sin(DEG2RAD*THETA_M))), focus_yh = 0.3, radius=0.11,
	lambda0 = lambda, dlambda = lambda*0.01,
	T1=683.7,I1=0.5874e+13,T2=257.7,I2=2.5099e+13,T3=16.7 ,I3=1.0343e+12)
  AT (0, 0, 0) RELATIVE Origin

/* octogonal slit, and further rectangular slits */
COMPONENT Win0 = Slit(radius = 0.11)                   AT (0, 0, 2.462) RELATIVE Origin
COMPONENT Oct1 = Slit(xwidth = 0.15, yheight = 0.15)   AT (0, 0, 2.512) RELATIVE Origin
COMPONENT Oct2 = Slit(xwidth = 0.177, yheight = 0.177) AT (0, 0, 2.512) RELATIVE Origin ROTATED (0,0,45) RELATIVE Origin
COMPONENT Win1 = Slit(radius = 0.0565)  AT (0, 0, 3.5321) RELATIVE Origin
COMPONENT Win2 = Slit(radius = 0.0555)  AT (0, 0, 4.0421) RELATIVE Origin
COMPONENT Win3 = Slit(radius = 0.0575)  AT (0, 0, 4.087)  RELATIVE Origin
COMPONENT Win4 = Slit(radius = 0.0575)  AT (0, 0, 4.999)  RELATIVE Origin
COMPONENT Win5 = Slit(radius = 0.0625)  AT (0, 0, 5.006)  RELATIVE Origin
COMPONENT Win6 = Slit(radius = 0.0625)  AT (0, 0, 5.151)  RELATIVE Origin

COMPONENT SMguide = Guide(w1=0.07, h1=0.3, m=coating, l=5)
  WHEN (coating>0)
  AT (0,0,11) RELATIVE Origin

/* ALPHA1 collimator */
COMPONENT Alpha_One = Collimator_linear(
  xwidth = 0.07, yheight = 0.30, length = 0.70, divergence = ALPHA1, transmission = 1.0)
  WHEN (ALPHA1>0)
  AT (0, 0, L1-1.045) RELATIVE Origin

COMPONENT Slit_Mono = Slit(xwidth = 0.05, yheight = 0.3)
  AT (0, 0, L1-0.1) RELATIVE Origin

COMPONENT D2B_Mono_XY = Monitor_nD(
     options="x y",
     bins = 40, xwidth = 0.1, yheight = 0.3, restore_neutron=1)
   AT (0, 0, 0) RELATIVE PREVIOUS
COMPONENT D2B_ALPHA1_Div = Monitor_nD(
     options="dx limits=[-0.5 0.5] dy limits=[-1 1]",
     bins = 40, xwidth = 0.1, yheight = 0.3, restore_neutron=1)
   AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT D2B_Mono_Lambda = Monitor_nD(
     options=lamopts,
     bins = 40, xwidth = 0.1, yheight = 0.3, restore_neutron=1)
   AT (0, 0, 0) RELATIVE PREVIOUS

/* TIP: monochromator cradle */
COMPONENT mono_cradle = Arm()
  AT (0, 0, L1) RELATIVE Origin

/* TIP: curved monochromator with NH>1 NV>1 et RH>0 RV>0 */
SPLIT COMPONENT Monok = Monochromator_curved(
    width = 0.1, height = 0.3, NH = 1, NV = 28, RV=RV,
    mosaich = ETA, mosaicv = 12, DM = DM)
  AT (0, 0, 0) RELATIVE mono_cradle
  ROTATED (0, THETA_M+TILT, 0) RELATIVE mono_cradle

COMPONENT sphere = PSD_monitor_4PI(filename="sphere", radius=0.5, restore_neutron=1)
  AT (0, 0, 0) RELATIVE mono_cradle

  
/* TIP: positioning diffraction direction for monok (order 1) */
COMPONENT mono_out = Arm()
  AT (0, 0, 0) RELATIVE mono_cradle
  ROTATED (0, 2*THETA_M, 0) RELATIVE mono_cradle

COMPONENT Shield_Hole = Guide_channeled(
    w1 = 0.05, h1 = 0.274, w2 = 0.044, h2 = 0.201, l = 0.752, R0 = 0, nslit = 1)
  AT (0, 0, 0.304) RELATIVE mono_out

/* ALPHA2 collimator */
COMPONENT Alpha_Two = Collimator_linear(
    xwidth = 0.05, yheight = 0.2, length = 0.2, divergence = ALPHA2, transmission = 1.0)
  WHEN (ALPHA2>0)
  AT (0, 0, L2-0.25) RELATIVE mono_out

COMPONENT Slit_Sample = Slit(xwidth = 0.01, yheight = 0.05)
  AT (0, 0, L2-0.02) RELATIVE mono_out

COMPONENT D2B_ALPHA2_Div = Monitor_nD(
     options="dx limits=[-0.5 0.5] dy limits=[-1 1]",
     bins = 40, xwidth = 0.02, yheight = 0.05, restore_neutron=1)
   AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT D2B_Sample_Lambda = Monitor_nD(
     options=lamopts,
     bins = 40, xwidth = 0.02, yheight = 0.05, restore_neutron=1)
   AT (0, 0, 0) RELATIVE PREVIOUS

/* sample position */
SPLIT COMPONENT Arm_Sample = Arm()
  AT (0, 0, L2) RELATIVE mono_out

COMPONENT Container_in = PowderN(reflections="V.lau", radius = 0.0025+0.00051, thickness = 0.0005, yheight = 0.04,
                         concentric = 1, p_interact=0.1, d_phi=RAD2DEG*atan2(Dheight+0.05,L3))
  AT (0, 0, 0) RELATIVE Arm_Sample

COMPONENT Sample = PowderN(
    reflections = Powder, radius = 0.0025, d_phi=RAD2DEG*atan2(Dheight+0.05,L3),
    yheight = 0.04)
  AT (0, 0, 0) RELATIVE Arm_Sample

COMPONENT Container_out = COPY(Container_in)(concentric=0)
  AT (0, 0, 0) RELATIVE Arm_Sample

/* ALPHA3 radial collimator */
COMPONENT collimador_radial = Collimator_radial(
    yheight=Dheight, length=.30,
    divergence=ALPHA3,verbose=1,approx=1,
    theta_min=-165, theta_max=-5, radius=L3-0.3-0.01)
  WHEN (ALPHA3>0)
	AT (0, 0, 0) RELATIVE Sample

/* perfect detector: 1D(theta) */
COMPONENT D2B_BananaTheta = Monitor_nD(
    options = "banana, angle limits=[5 165], bins=3200",
    radius = L3, yheight = Dheight, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample

COMPONENT D2B_BananaThetaCenter = Monitor_nD(
    options = "banana, angle limits=[5 165], bins=3200",
    radius = L3, yheight = 0.1, restore_neutron=1)
  AT (0, 0, 0) RELATIVE Sample

/* perfect detector: 2D(theta,y) to see diffraction rings */
COMPONENT D2B_BananaPSD = Monitor_nD(
    options = "banana, theta limits=[-165 -5] bins=320, y bins=25",
    radius = L3*1.005, yheight = Dheight)
  AT (0, 0, 0) RELATIVE Sample

/* The END token marks the instrument definition end */
END