/************************************************
*
* 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
* Release: McStas 1.12c
* Version: $Revision: 1.4 $
* %INSTRUMENT_SITE: PSI
*
* %Description
* McStas model of the DMC powder diffractometer at PSI, CH.
*
* Please note that this instrumentfile 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: 
* mcrun PSI_DMC.instr lambda=2.5666 FileName="Na2Ca3Al2F14.laz" -n1e8
*
* %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 filename, 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(lambda = 2.5666, R = 0.87, R_curve = 0.87, string FileName="Cd.laz",D_PHI=30, SHIFT=0, PACK=0.7, Dw=0.8, BARNS=1 , sample_radius = 0.002)

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_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(
        height=0.156, width=0.126,
	l_low=lambda-ldiff/2,l_high=lambda+ldiff/2,
        dist=1.5, xw = 0.025, 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(
    nchan = 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(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(
        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(
        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(
        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(
    win_thick = 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(
    nx = 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(
    nx = 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(
    win_thick = 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, nchan=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

SPLIT 10 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, nchan=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 

SPLIT 10 COMPONENT sample = PowderN (
	d_phi=D_PHI,radius=sample_radius, h=sample_height, DW=Dw,
        pack=PACK, reflections = FileName, barns=BARNS, tfrac=0, frac=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,
  options="banana, theta limits [19.9+SHIFT 99.9+SHIFT] bins=400, file=detector.dat")
AT (0,0,0) RELATIVE sa_arm
ROTATED (0, 0, 180) RELATIVE sa_arm


END