/************************************************
*
* 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: PSDbefore_guides_I=2.54332e+09
*
* %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(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 = Arm()
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 Mod_Face = Shape(xwidth=0.1,yheight=0.1)
AT (0,0,-0.16) RELATIVE ABSOLUTE
EXTEND %{
  /* Propagate back to a small rectangle in front of moderator */
  printf("(before) x=%g y=%g z=%g\n", x, y, z);
  
  ALLOW_BACKPROP;
  PROP_Z0;
  SCATTER;
  
  printf("(after) x=%g y=%g z=%g\n", x, y, z);  
  printf("(if) %d\n", (fabs(x)> xwidth/2 || fabs(y)> yheight/2));
  /* Remove neutrons that are not from around the moderators */
  if (fabs(x)> xwidth/2 || fabs(y)> yheight/2) 
  {
    printf("(absorbed) x=%g y=%g z=%g\n", x, y, z);   
    ABSORB;   
  }
%}


  
COMPONENT PSDbefore_guides = PSD_monitor(
    nx = 128, ny = 128, filename = "PSDbefore_guides",
    xwidth = 0.025, yheight = 0.125)
AT (0, 0, 1.5) RELATIVE source_arm

END