/************************************************
*
* McStas model of the DMC powder diffractometer at PSI, CH.
*
* %Identification
* Written by: Peter Willendrup (Risoe), Uwe Filges (PSI), Lukas Keller (PSI)
* 
* Date: May 7th, 2008
* Origin:PSI
* %INSTRUMENT_SITE: PSI
*
* McStas model of the DMC powder diffractometer at PSI, CH.
*
* %Description
* McStas model of the DMC powder diffractometer at PSI, CH.
*
* Please note that this instrument file does not include the radial collimator of DMC. 
* To generate the full 800-bin angle coverage at DMC, you should combine two simulations
* using this instrumentfile, with SHIFT=0 and 0.1. This will displace the detector by half
* a bin-width, which is a standard procedure at the DMC diffractometer
*
* %Example: lambda=2.5666 Detector: Detector_I=7.5965E+02
*
* %Parameters
* lambda: [AA]     Choice of wavelength, affects both monochromator and source component
* R: [1]           Reflectivity of non-curved guides
* R_curve: [1]     Reflectivity of curved guides
* filename: [str]  Choice of reflection list file, e.g. from McStas data dir
* D_PHI: [deg]     Focusing 'd_phi' for PowderN, see mcdoc page
* SHIFT: [deg]     Rotation of detector, set to 0.1 to displace by half a bin
* PACK: [1]        Powder packing factor
* Dw: [1]          Powder Debye-Waller factor
* BARNS: [1]       Flag to define if powder reflection file |F2| is in Barns or fm
*
* %Link
* The <a href="http://lns00.psi.ch/montecarlo/">PSI Monte-Carlo</a> website.
* %Link
* P Willendrup, U Filges, L Keller, E Farhi, K Lefmann: Validation of a realistic powder sample using data from DMC at PSI, ICNS 2005 (Physica B, 386, (2006), 1032.)
*
* %End
***************************************************/
DEFINE INSTRUMENT PSI_DMC_simple(lambda=2.5666, R=0.87, R_curve=0.87, string filename="Na2Ca3Al2F14.laz",D_PHI=6, SHIFT=0, PACK=0.7, Dw=0.8, BARNS=1)

DECLARE
%{
  double mono_q = 1.8734;
  double OMA;  
  double RV;
  double y_mono = 0.025;
  double NV = 5;
  double d_phi_0;
  double TTM;
  double sample_radius = 0.008/2;
  double sample_height = 0.03;
  double can_radius = 0.0083/2;
  double can_height = 0.0303;
  double can_thick = 0.00015;
  
  /******Mirrorvalues*****/
  
  double alpha;
  double Qc=0.0217;
  double R0=0.995;
  double Mvalue=1.9;
  double W=1.0/250.0;
  
  double alpha_curve;
  double Qc_curve=0.0217;
  double R0_curve= 0.995;
  double Mvalue_curve=2.1;
  double W_curve=1.0/250.0;
  
  double ldiff=0.05;
  /* Curved guide element angle*/
  double angleGuideCurved;

%}

INITIALIZE
%{
  TTM = 2*asin(mono_q*lambda/(4*PI))*RAD2DEG;
  OMA = TTM/2;
  RV = fabs(2*2.82*sin(DEG2RAD*OMA));
  
  angleGuideCurved=-2.0*asin(0.4995 /2.0/3612)/PI*180;
  alpha=(R0-R)/Qc/(Mvalue-1);
  alpha_curve=(R0_curve-R_curve)/Qc_curve/(Mvalue_curve-1);
  
%}

TRACE

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

COMPONENT source = Source_Maxwell_3(
        yheight=0.156, xwidth=0.126,
	      Lmin=lambda-ldiff/2,Lmax=lambda+ldiff/2,
        dist=1.5, focus_xw = 0.02, focus_yh = 0.12,
        T1=296.16, I1=8.5E11,
        T2=40.68, I2=5.2E11)    
AT (0,0,0) RELATIVE source_arm  ROTATED (0,0,0) RELATIVE source_arm

COMPONENT PSDbefore_guides = PSD_monitor(
    nx = 128, ny = 128, filename = "PSDbefore_guides",
    xwidth = 0.02, yheight = 0.12)
AT (0, 0, 1.49999) RELATIVE source_arm

COMPONENT l_mon_source = L_monitor(
    nL = 101, filename = "lmonsource.dat", xwidth = 0.02,
    yheight = 0.12, Lmin = 0, Lmax = 20)
AT (0, 0, 1e-9) RELATIVE PREVIOUS


/* guide segment 1, m=2, 4.66 m */
COMPONENT guide1 = Guide_simple(w1= 0.02, h1=0.12, w2=0.02, h2=0.12,
                         l=4.66, R0=R0, Qc=Qc, alpha = alpha,
                         m = 1.8, W = W)
AT (0,0,1.50) RELATIVE source_arm ROTATED (0,0,0) RELATIVE source_arm

COMPONENT PSDbefore_curve = PSD_monitor(
    nx = 128, ny = 128, filename = "PSDbefore_curve",
    xwidth = 0.02, yheight = 0.12)
AT (0, 0, 4.664) RELATIVE guide1

COMPONENT guide2 = Bender(
    w = 0.02, h = 0.12, r = 3612, R0a = R0_curve, Qca = Qc_curve,
    alphaa = alpha_curve, ma = Mvalue_curve, Wa = W_curve, R0i = R0_curve, Qci = Qc_curve,
    alphai = alpha_curve, mi = 1, Wi = W_curve, R0s = R0_curve, Qcs = Qc_curve,
    alphas = alpha_curve, ms = Mvalue_curve, Ws = W_curve, l = 20)
AT (0, 0, 4.69) RELATIVE guide1

COMPONENT PSDafter_curve = PSD_monitor(
    nx = 128, ny = 128, filename = "PSDafter_curve",
    xwidth = 0.02, yheight = 0.12)
  AT (0, 0, 20.0001) RELATIVE guide2

/* bunker wall, m=2, 3.0 m */
COMPONENT bunker = Guide_simple(
        w1=0.02, h1=.12, w2=0.02, h2=.12,
        l=3.43, R0=R0,Qc=Qc, alpha = alpha,
        m = 1.6, W = W)
AT (0,0,20.1502) RELATIVE guide2 ROTATED (0,0,0) RELATIVE guide2

/* guide segment 3, m=2, 8.84 m */

COMPONENT guide3 = Guide_simple(
        w1=0.02, h1=.12, w2=0.02, h2=.12,
        l=12.275, R0=R0,Qc=Qc, alpha = alpha,
        m = 1.6, W = W)
AT (0,0,3.56) RELATIVE bunker ROTATED (0,0,0) RELATIVE bunker

COMPONENT guide4 = Guide_simple(
        w1=0.02, h1=.12, w2=0.02, h2=.12,
        l=5.66, R0=R0,Qc=Qc, alpha = alpha,
        m = 1.6, W = W)
AT (0,0,15.8555) RELATIVE bunker ROTATED (0,0,0) RELATIVE guide3

COMPONENT window1 = Al_window(
    thickness = 0.002)
AT (0, 0, 5.66+1e-9) RELATIVE PREVIOUS

/* guide exit at 42.5 m - measured flux 2.46 10e8 cm-2 s-1 mA-1 */

COMPONENT ydist_fluxpos = PSDlin_monitor(
    nbins = 11, filename = "ydist_fluxpos.dat", xwidth = 0.120,
    yheight = 0.02)
AT (0, 0, 5.66+1e-8+0.01) RELATIVE guide4 ROTATED (0,0,90) RELATIVE PREVIOUS

COMPONENT PSD_fluxpos = PSD_monitor(
    nx = 100, ny = 100, filename = "xdist_fluxposy.dat",
    xwidth = 0.02, yheight = 0.12)
AT (0, 0, 5.66+1e-7+0.01) RELATIVE guide4

COMPONENT xdist_flux_pos =  PSDlin_monitor(
    nbins = 11, filename = "xdist_fluxpos.dat", xwidth = 0.020,
    yheight = 0.12)
AT (0, 0, 1e-9) RELATIVE PREVIOUS

COMPONENT PSD_fluxposB = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_fluxposB.dat",
    xwidth = 0.02, yheight = 0.12)
AT (0, 0, 6.24-1e-7-0.01) RELATIVE guide4

COMPONENT window2 = Al_window(
    thickness = 0.002)
    AT (0, 0, 1e-9) RELATIVE PREVIOUS

/*0.580 m gap for DMC monochromator*/

COMPONENT in_slit = Slit(
        xmin = -0.01, xmax=0.01 , ymin = -0.06, ymax = 0.06)
  AT (0,0,0.0021) RELATIVE  PREVIOUS

COMPONENT lambda_in =  L_monitor(
    xmin=-0.011, xmax=0.011, ymin=-0.061, ymax=0.061, Lmin=0, Lmax=2*lambda, nL=128, filename="L_in.dat")
  AT ( 0, 0, 0.001) RELATIVE in_slit

/* Monochromator description */ 

COMPONENT sma = Arm() /* source - monochromator arm */ 
AT (0, 0, 0.65) RELATIVE in_slit ROTATED (0,OMA,0) RELATIVE in_slit

COMPONENT foc_mono = Monochromator_2foc(
    zwidth = 0.05, yheight = 0.025, gap = 0.0005, NH = 1, NV = 5,
    mosaich = 38, mosaicv = 38, r0 = 0.7, Q = mono_q, RV = RV,
    RH = 0)
AT (0, 0, 0) RELATIVE sma

COMPONENT msa = Arm() /* monochromator - sample arm */
   AT (0, 0, 0) RELATIVE sma ROTATED (0, TTM, 0) RELATIVE in_slit

COMPONENT out1_slit = Slit(
        xmin=-0.01, xmax=0.01, ymin=-0.06, ymax=0.06)
  AT (0,0,0.2) RELATIVE msa ROTATED (0,0,0) RELATIVE msa  

COMPONENT Amoin_slit = Slit(
        xmin=-0.01, xmax=0.01, ymin=-0.06, ymax=0.06)
  AT (0,0,0.325) RELATIVE  msa ROTATED (0,0,0) RELATIVE msa

COMPONENT Bmoin_slit = Slit(
        xmin=-0.01, xmax=0.01, ymin=-0.06, ymax=0.06) 
  AT (0,0,0.525) RELATIVE  msa ROTATED (0,0,0) RELATIVE msa

COMPONENT out2_slit = Slit(
        xmin=-0.01, xmax=0.01, ymin=-0.06, ymax=0.06)
  AT (0,0,0.65) RELATIVE msa ROTATED (0,0,0) RELATIVE msa  

COMPONENT PSD_sample =  PSD_monitor(
    xmin=-0.05, xmax=0.05, ymin=-0.07, ymax=0.07,
    nx=80, ny=80, filename="PSD_sample.dat")
  AT ( 0, 0, 2.77) RELATIVE msa

COMPONENT lambda_sample =  L_monitor(
    xmin=-sample_radius, xmax=sample_radius, ymin=-sample_height/2, ymax=sample_height/2, Lmin=lambda-0.2, Lmax=lambda+0.2, nL=128, filename="L_sample.dat")
  AT ( 0, 0, 2.81) RELATIVE msa

COMPONENT sa_arm = Arm()
   AT (0, 0, 2.82) RELATIVE msa 
   ROTATED (0, 0, 0) RELATIVE msa 

COMPONENT sample = PowderN (
	d_phi=D_PHI,radius=sample_radius, yheight=sample_height, DW=Dw,
  pack=PACK, reflections = filename, barns=BARNS, p_transmit=0, p_inc=0)
AT ( 0, 0, 0) RELATIVE sa_arm

COMPONENT STOP = Beamstop(radius=0.3)
AT (0,0,1.4) RELATIVE sa_arm
ROTATED (0,0,0) RELATIVE sa_arm


COMPONENT Detector = Monitor_nD(
  xwidth=3.0, yheight=0.09, filename="detector.dat", min=19.9+SHIFT, max=99.9+SHIFT, bins=400,
  options="banana, theta")
AT (0,0,0) RELATIVE sa_arm
ROTATED (0, 0, 180) RELATIVE sa_arm


END