/*******************************************************************************
*   McStas simulation of the Powder Instrument for Transition in Structure Investigations (PITSI)
*
* Instrument: SAFARI_PITSI
*
* %Identification
* Written by: Deon Marais (deon.marais@necsa.co.za)
* Date: September 2013
* Origin: Necsa
* %INSTRUMENT_SITE: Necsa
*
* Powder Instrument for Transition in Structure Investigations
*
* %Description
* Necsa Neutron Powder Diffractometer located at beam port 5 of the SAFARI-1 research reactor, South Africa
*
* %Parameters
* source_lam_min: [Angs]  Minimum wavelenth of source
* source_lam_max: [Angs]  Maximum wavelenth of source
* hi_res: []                 Selects hi-resolution(1) or hi-intensity(0) reactor beam through primary shutter
* mono_Si_type: []           Monochromator Silicon type options: 422 400 311 511 111 331 (Mostly used: 331 and 551)
* mono_mosh: [arc min]    Monochromator horizontal mosaicity
* mono_mosv: [arc min]    Monochromator vertical mosaicity
* mono_dx: [m]            Monochromator delta x - positive to left of reactor beam
* mono_dy: [m]            Monochromator delta y - positive upward of reactor beam
* mono_dz: [m]            Monochromator delta z - positive along reactor beam
* mono_takeoff: [deg]     Monochromator takeoff angle - positive anti-clockwise from reactor beam
* mono_dtilt: [deg]       Monochromator tilt angle - not implemented yet
* mono_r_h: [m]           Monochromator horizontal focus radius (min 5.0 m)
* mono_r_v: [m]           Monochromator vertical focus radius (min 0.9 m)
* port_takeoff: [deg]     Port takeoff angle - positive anti clockwise from reactor beam (70 or 90 deg)
* inc_slit_rot: [deg]     Incident slit delta rotation - not implemented yet
* inc_slit_dx: [m]        Incident slit delta x position - positive to left of incident beam
* inc_slit_to_cor: [m]    Incident slit to sample stage center of rotation
* inc_slit_width: [m]     Incident slit width 0.00013m to 0.006m
* inc_slit_height: [m]    Incident slit height 0m to 0.05m
* inc_slit_sep: []           Incident slit separation between width and height. <0:use emperical calc, >=0:distance in m
* mono_to_cor: [m]        Distance between monochromator and center of rotation
* sample_dx: [m]          Sample delta x - positive to left of incident beam if sample_dom=0
* sample_dy: [m]          Sample delta y - positive upword of incident beam
* sample_dz: [m]          Sample delta x - positive along of incident beam if sample_dom=0
* sample_dom: [deg]       Sample delta omega - positive anti-clockwise from incident beam
* det_takeoff: [deg]      Detector takeoff angle - positive anti-clockwise from incident beam
* cor_to_det: [m]         Distance between sample centre of rotation and detector
* dangle_interest: [deg]  Delta angle of interenterest.
* hi_res: []
* mono_Si_type: []
* inc_slit_sep: []
* full_instrument: []
* full_instrument: []        When 1, simulates the complete instrument. When 0, only simulate from the outlet collimator
*
* %Link
* The South African Nuclear Energy Corporation <a href="http://www.necsa.co.za"> website</a>
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT SAFARI_PITSI(
  source_lam_min=0.5,
  source_lam_max=2.0,
  hi_res=0, mono_Si_type=551,
  mono_mosh=30,
  mono_mosv=30,
  mono_dx=0,
  mono_dy=0,
  mono_dz=0,
  mono_takeoff=90.0,
  mono_dtilt=0,
  mono_r_h=5.0,
  mono_r_v=3.572,
  port_takeoff=90.0,
  inc_slit_rot=0,
  inc_slit_dx=0,
  inc_slit_to_cor=0.01,
  inc_slit_width=0.006,
  inc_slit_height=0.05,
  inc_slit_sep=0,
  mono_to_cor=2.5,
  sample_dx=0,
  sample_dy=0,
  sample_dz=0,
  sample_dom=0,
  det_takeoff=-114.375,
  cor_to_det=1.179,
  dangle_interest=125,
  full_instrument=1)

DECLARE
%{
  double hi_res, port_takeoff;

  double mono_Si_type, mono_mosh, mono_mosv;
  double mono_dx, mono_dy, mono_dz;
  double mono_takeoff, mono_dtilt, mono_r_h, mono_r_v;
  double mono_r_req_v;
  double mono_pts_v, mono_pts_h, focal_dist;
  double a,b,c,d,y;      //Gaussian fit params

  double inc_slit_rot, inc_slit_dx, inc_slit_to_cor;
  double inc_slit_width, inc_slit_height, inc_slit_sep;
  double mono_to_cor;
  double sample_dx, sample_dy, sample_dz, sample_dom;

  double det_takeoff, cor_to_det, dangle_interest, ndet, det_cover_angle;
  double det_width, det_height;

  double diff_slit_dx, diff_slit_to_cor;
  double diff_slit_width, diff_slit_height;

  double inc_slit_xmin, inc_slit_xmax, inc_slit_ymin, inc_slit_ymax;
  double inc_slit_xmin_h, inc_slit_xmax_h, inc_slit_ymin_w, inc_slit_ymax_w;
  double diff_slit_xmin, diff_slit_xmax, diff_slit_ymin, diff_slit_ymax;
  double wafer_d, start_wafer_pos, mono_turns, mono_Rh_req, mono_d, mono_q;
  double lam;
  double as;
  int msw;

  double chamber_col_start, chamber_col_length, outside_chamber_collimator_w, outside_chamber_collimator_h;

  double from_col=1;
  double full_instrument;
%}

INITIALIZE
%{
  printf ("\n------------------\n");

  /* Constants */
  if (port_takeoff == 70.0) {
    chamber_col_start=0.44057;   //from the monochromator diffraction center
    chamber_col_length=1.464;
    outside_chamber_collimator_w=0.01954;
    outside_chamber_collimator_h=0.04099;
  } else {        //90 degrees
    chamber_col_start=0.414;   //from the monochromator diffraction center
    chamber_col_length=1.3458;
    outside_chamber_collimator_w=0.02186;
    outside_chamber_collimator_h=0.04853;
  }

  /* Incident slit */
  inc_slit_xmin = -inc_slit_width/2.0;
  inc_slit_xmax = inc_slit_width/2.0;
  inc_slit_ymin = -inc_slit_height/2.0;
  inc_slit_ymax = inc_slit_height/2.0;
  if (inc_slit_sep < 0) {      //Emperical formula. Linear dependance derived from measurements: sep=13.18mm when width=0.13mm; sep=4.18mm when width=4.95mm
    inc_slit_sep= -1.86721992 * inc_slit_width + 0.0134227385;
  }
  printf ("Incident slit separation = %.4lf mm\n",inc_slit_sep*1000);

  if (inc_slit_sep == 0) {      //No gap, so horisonal and vertical slits lie ontop of each other and can have same dimensions
    inc_slit_xmin_h = inc_slit_xmin;
    inc_slit_xmax_h = inc_slit_xmax;
    inc_slit_ymin_w = inc_slit_ymin;
    inc_slit_ymax_w = inc_slit_ymax;
  } else {        //Perform some trig to be sure that none of the neutrons are accidentally cut of
    double col_to_slit = (mono_to_cor-inc_slit_to_cor) - (mono_to_cor-chamber_col_start+chamber_col_length);
    inc_slit_ymax_w = 1.01*(inc_slit_ymax + (0.023745 + inc_slit_ymax)*inc_slit_sep/col_to_slit);
    if (inc_slit_ymax_w < inc_slit_ymax) inc_slit_ymax_w = inc_slit_ymax;
    inc_slit_ymin_w = -inc_slit_ymax_w;
    inc_slit_xmax_h = -1.01*(((0.01077+inc_slit_xmax)*inc_slit_sep/col_to_slit)-inc_slit_xmax);
    if (inc_slit_xmax_h < inc_slit_xmax) inc_slit_xmax_h = inc_slit_xmax;
    inc_slit_xmin_h = - inc_slit_xmax_h;
  }
  printf ("Height slit: width=%.4lfmm, height=%.4lfmm\n",2*inc_slit_xmax_h*1000, 2*inc_slit_ymax*1000);
  printf ("Width  slit: width=%.4lfmm, height=%.4lfmm\n",2*inc_slit_xmax*1000, 2*inc_slit_ymax_w*1000);


  /*
  PoiConst = 0.2;
  miu = 0.00004;
  F = 4.1534*0.01;
  D = 0.000725;*/

  /* Monochromator */
  wafer_d = (6.0/13.0/1000.0);    //Wafer diameter. Total blade diameter (6mm) made up from 13? wafers
  start_wafer_pos = 6.0 * wafer_d;  //In order to have the middle waver situated at the centre of the port takeoff

  msw=(int)mono_Si_type;
  mono_d=0.0;
  if (msw==422) mono_d = 1.10858;
  if (msw==400) mono_d = 1.35773;
  if (msw==311) mono_d = 1.63748;
  if (msw==511) mono_d = 1.04518;
  if (msw==111) mono_d = 3.135;
  if (msw==331) mono_d = 1.24594;
  if (msw==551) mono_d = 0.76049;

  mono_q = 2*PI/mono_d;
  double mono_omega = fabs(mono_takeoff/2.0);
  lam=2.0*mono_d*sin(DEG2RAD*mono_omega);
  printf ("mono_Si_type = %i, mono_d=%.4lfA, mono_omega = %.2lfdeg, lambda = %.4lfA\n",msw, mono_d, mono_omega, lam);

  /* Monochromator-sample distance Vertical*/
  mono_pts_v=55.55+(805.55-55.55)/(1.0/0.9)*(1.0/mono_r_v); //emperical formula that equates number of Control Pts to the radius. Derived from '13-09-23 Necsa Si Monochromator Data Sheet'
  as = -tan(2*DEG2RAD*mono_omega)/(tan(DEG2RAD*mono_omega));
  focal_dist = mono_r_v*(1-(0.5/as))*sin(DEG2RAD*mono_omega);
  printf ("Monochromator Vertical to focal point = %.4lfm (%.2lf Control Pts)\n",focal_dist, mono_pts_v);
  double mono_r_req_v = mono_to_cor/((1-(0.5/as))*sin(DEG2RAD*mono_omega));
  mono_pts_v=55.55+(805.55-55.55)/(1.0/0.9)*(1.0/mono_r_req_v); //emperical formula that equates number of Control Pts to the radius. Derived from '13-09-23 Necsa Si Monochromator Data Sheet'
  printf ("To Vertical focus on COR at %.3lfm, let mono_r_v = %.4lfm (%.2lf Pts)\n",mono_to_cor,mono_r_req_v, mono_pts_v);

  /* Monochromator-sample distance Horisontal*/
  a = 0.227642476966334;  //Gausian fitted parameters for Control Pts against curvature. y=a*exp((x-b)^2/(-2*c^2))+d
  b = 427.435873967018;  //
  c = 141.018775963938;
  d = 0.000855170045959246;
  y = 1.0/mono_r_h;    //curvature
  mono_pts_h = b + c * sqrt(2*log((y-d)/a));
  as = -tan(2*DEG2RAD*mono_omega)/(tan(DEG2RAD*mono_omega));
  focal_dist = mono_r_h*(1-(0.5/as))*sin(DEG2RAD*mono_omega);
  printf ("Monochromator Horizontan to focal point = %.4lfm (%.2lf Control Pts)\n",focal_dist, mono_pts_h);

  //mono_Rh_req = mono_to_cor/((1-(0.5/as))*sin(DEG2RAD*mono_omega));
  //mono_turns=(1005.5/mono_Rh_req)-13;    //emperical formula that equates the number of turns of the motor to the curvature. Taken from Multi-wafer silicon monochromator for stress machine at SAFARI, South Africa. Mihai Popovici
  //printf ("To focus on COR at %.3lfm, let mono_r_h = %.4lfm (%.2lf turns)\n",mono_to_cor,mono_Rh_req, mono_turns);

  /* Number of detectors to construct and calculates the angle offset between them */
  det_width=0.66;
  det_height=0.38;
  det_cover_angle = 2.0*RAD2DEG*atan((det_width/2.0)/cor_to_det);
  ndet = ceil(dangle_interest/det_cover_angle);
  printf ("Single detector coverage:%.2lf deg. %i detectors coverage angle: %.2f deg\n",\
    det_cover_angle, (int)ndet, ndet*det_cover_angle);


  printf ("------------------\n\n");
%}

TRACE
//**************************************************************************************************
//Source
COMPONENT Progress = Progress_bar(percent=1,flag_save=0)
  AT (0,0,0) ABSOLUTE

COMPONENT Reactorbeam = Arm()
  AT (0,0,0) ABSOLUTE

COMPONENT Prim_axes = Arm()
  AT (0, 0, 5.140) RELATIVE Reactorbeam
  ROTATED (0, port_takeoff, 0) RELATIVE Reactorbeam

COMPONENT  Source = Source_gen(
    radius = 0.0905, dist = 2.86805, focus_xw = 0.1, focus_yh = 0.05,
    Lmin = source_lam_min, Lmax = source_lam_max, I1 = 0)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Reactorbeam

COMPONENT PSD_Source = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_Source",
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Reactorbeam

COMPONENT LAM_Source = L_monitor(
    nL = 100, Lmin = source_lam_min, Lmax = source_lam_max,
    filename = "LAM_Source.out",
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Reactorbeam


COMPONENT Window_before_filter = Slit(
    //xmin = -0.05, xmax = 0.05, ymin = -0.02922, ymax = 0.02922)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 3.24045) RELATIVE Reactorbeam

//**************************************************************************************************
//Filter
// The Al2O3_sapphire.trm filter is specified for a 0.0508m thickness,
// ours is 0.15876m thick, therefore thickness will be defined as 3.1752
COMPONENT Sapphire_filter = Filter_gen(
    filename = "Al2O3_sapphire.trm", options = "multiply",
    xmin = -0.053975, xmax = 0.053975 , ymin = -0.0381, ymax = 0.0381, thickness=3.125)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Window_before_filter

COMPONENT Window_after_filter = Slit(
    xmin = -0.05, xmax = 0.05, ymin = -0.032005, ymax = 0.032005)
  WHEN (full_instrument==1)
  AT (0, 0, 0.15876) RELATIVE Sapphire_filter

COMPONENT LAM_After_sapphire = L_monitor(
    nL = 100, Lmin = source_lam_min, Lmax = source_lam_max,
    filename = "LAM_After_sapphire.out",
    //xmin = -0.05, xmax = 0.05, ymin = -0.032005, ymax = 0.032005)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Window_after_filter

COMPONENT PSD_After_sapphire = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_After_sapphire",
    //xmin = -0.05*1.1, xmax = 0.05*1.1, ymin = -0.032005*1.1, ymax = 0.032005*1.1)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Window_after_filter


COMPONENT HighResOutlet = Slit(
    xmin = -0.0325, xmax = 0.0325, ymin = -0.05, ymax = 0.05)
  WHEN ((hi_res==1) && (full_instrument==1))
  AT (0, 0, 4.84055) RELATIVE Reactorbeam

COMPONENT HighIntensityOutlet = Slit(
    xmin = -0.05, xmax = 0.05, ymin = -0.05, ymax = 0.05)
  WHEN ((hi_res==0) && (full_instrument==1))
  AT (0, 0, 4.84055) RELATIVE Reactorbeam

COMPONENT PSD_After_Outlet = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_After_Outlet",
    //xmin = -0.05*1.1, xmax = 0.05*1.1, ymin = -0.06*1.1, ymax = 0.06*1.1)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 4.84055) RELATIVE Reactorbeam

//**************************************************************************************************
//Monochromator
COMPONENT Mono_axis = Arm()
  AT (0+mono_dx, 0+mono_dy, 5.140+mono_dz) RELATIVE Reactorbeam
  ROTATED (0, mono_takeoff/2.0, 0) RELATIVE Reactorbeam


COMPONENT Blade_1 = Monochromator_curved(
   zwidth = 0.22/51.0, yheight = 0.1395/9.0, gap = 0.0, NH = 51, NV = 9,
   mosaich = mono_mosh,  mosaicv = mono_mosv, r0 = 1.0, t0 = 1.0, RV = mono_r_v, RH = mono_r_h,
   DM = mono_d )
  WHEN (full_instrument==1)
  AT ( -start_wafer_pos, 0,0) RELATIVE Mono_axis
  GROUP Monochro

COMPONENT Blade_2 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_3 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_4 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_5 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_6 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_7 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_8 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_9 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_10 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_11 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_12 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

COMPONENT Blade_13 = COPY(Blade_1)
  WHEN (full_instrument==1)
  AT (wafer_d,0,0) RELATIVE PREVIOUS
  GROUP Monochro

/*COMPONENT PSD_BehindMono = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_BehindMono.out",
    //xmin = -0.13, xmax = 0.13, ymin = -0.1125, ymax = 0.1125)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  AT (0, 0, 5.8) RELATIVE Reactorbeam
  GROUP Monochro

COMPONENT LAM_BehindMono = L_monitor(
    nL = 100, Lmin = source_lam_min, Lmax = source_lam_max,
    filename = "LAM_BehindMono.out",
    //xmin = -0.13, xmax = 0.13, ymin = -0.1125, ymax = 0.1125)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  AT (0, 0, 5.8) RELATIVE Reactorbeam
  GROUP Monochro

COMPONENT ReactorBeamStop = Beamstop(
    xmin = -0.1,xmax = 0.1, ymin=-0.1, ymax=0.1)
  AT (0, 0, 5.8) RELATIVE Reactorbeam*/





//**************************************************************************************************
//Chamber Collimator
//COMPONENT Prim_axes = Arm()
//  AT (0, 0, 5.55177) RELATIVE Reactorbeam
//  ROTATED (0, port_takeoff, 0) RELATIVE Reactorbeam

  //About at the inside surface of the secondary shutter
COMPONENT PSD_At_sec_shutter = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_At_sec_shutter.out",
    //xmin = -0.13, xmax = 0.13, ymin = -0.1125, ymax = 0.1125)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, chamber_col_start-0.1736) RELATIVE Prim_axes

COMPONENT LAM_At_sec_shutter = L_monitor(
    nL = 100, Lmin = source_lam_min, Lmax = source_lam_max,
    filename = "LAM_At_sec_shutter.out",
    //xmin = -0.13, xmax = 0.13, ymin = -0.1125, ymax = 0.1125)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT Inside_chamber_collimator = Slit(
      xmin = -0.021, xmax = 0.021, ymin = -0.0425, ymax = 0.0425)
  WHEN (full_instrument==1)
  AT (0, 0, chamber_col_start) RELATIVE Prim_axes

COMPONENT PSD_After_inside_chamber_collimator = PSD_monitor(
      nx = 100, ny = 100, filename = "PSD_After_inside_chamber_collimator.out",
      //xmin = -0.021*1.1, xmax = 0.021*1.1, ymin = -0.0425*1.1, ymax = 0.0425*1.1)
      xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Inside_chamber_collimator

COMPONENT Outside_chamber_collimator = Slit(
      xmin = -1*outside_chamber_collimator_w/2.0, xmax = outside_chamber_collimator_w/2.0,
      ymin = -1*outside_chamber_collimator_h/2.0, ymax = outside_chamber_collimator_h/2.0)
  WHEN (full_instrument==1)
  AT (0, 0, chamber_col_length) RELATIVE Inside_chamber_collimator

COMPONENT PSD_Outside_chamber_collimator_1 = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_Outside_chamber_collimator_1.out",
    //xmin = -0.012545*1.1, xmax = 0.012545*1.1, ymin = -0.02012*1.1, ymax = 0.02012*1.1)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  WHEN (full_instrument==1)
  AT (0, 0, 0) RELATIVE Outside_chamber_collimator

/*COMPONENT Source_capture = Virtual_output(
  filename="D:\\McStas\\PITSI_Source_capture.dat", bufsize=1e7)
   WHEN (full_instrument==1)
   AT (0, 0, chamber_col_start+chamber_col_length) RELATIVE Prim_axes
*/


/*COMPONENT Colimator_source = Virtual_input(
  filename="C:\\mcstas-2.0\\workspace\\Source_capture.dat",verbose=1, repeat_count=1, smooth=1)
   WHEN (full_instrument==0)
   AT (0, 0, chamber_col_start+chamber_col_length) RELATIVE Prim_axes */


COMPONENT PSD_Outside_chamber_collimator = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_Outside_chamber_collimator.out",
    //xmin = -0.012545*1.1, xmax = 0.012545*1.1, ymin = -0.02012*1.1, ymax = 0.02012*1.1)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  AT (0, 0, chamber_col_start+chamber_col_length) RELATIVE Prim_axes

//**************************************************************************************************
//Incident slit

COMPONENT Incident_slit_h = Slit(
      xmin = inc_slit_xmin_h, xmax = inc_slit_xmax_h, ymin = inc_slit_ymin, ymax = inc_slit_ymax)
  AT (inc_slit_dx, 0, mono_to_cor-inc_slit_to_cor-inc_slit_sep) RELATIVE Prim_axes

COMPONENT Incident_slit_w = Slit(
      xmin = inc_slit_xmin, xmax = inc_slit_xmax, ymin = inc_slit_ymin_w, ymax = inc_slit_ymax_w)
  AT (inc_slit_dx, 0, mono_to_cor-inc_slit_to_cor) RELATIVE Prim_axes

COMPONENT PSD_After_Incident_slit_w = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_After_Incident_slit_w.out",
    //xmin = inc_slit_xmin*1.5, xmax = inc_slit_xmax*1.5, ymin = inc_slit_ymin*1.5, ymax = inc_slit_ymax*1.5)
    xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
  AT (0, 0, 0) RELATIVE Incident_slit_w

COMPONENT Center_of_rotation = Arm()
  AT (0, 0, mono_to_cor) RELATIVE Prim_axes


//**************************************************************************************************
//Sample

COMPONENT Sample_rotation = Arm()
  AT (0, 0, 0) RELATIVE Center_of_rotation
  ROTATED (0, sample_dom, 0) RELATIVE Center_of_rotation

COMPONENT Sample_location = Arm()
  AT (sample_dx, sample_dy, sample_dz) RELATIVE Sample_rotation


COMPONENT PSD_Center_of_rotation = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_Center_of_rotation",
    //xmin = inc_slit_xmin*1.5, xmax = inc_slit_xmax*1.5, ymin = inc_slit_ymin*1.5, ymax = inc_slit_ymax*1.5)
    xmin = -0.035, xmax = 0.035, ymin = -0.035, ymax = 0.035)
  AT (0, 0, 0) RELATIVE Center_of_rotation

COMPONENT DIV_Center_of_rotation = Divergence_monitor(
    nh = 100, nv = 100, filename = "DIV_Center_of_rotation",
    //xmin = inc_slit_xmin*1.5, xmax = inc_slit_xmax*1.5, ymin = inc_slit_ymin*1.5, ymax = inc_slit_ymax*1.5)
    xmin = -0.035, xmax = 0.035, ymin = -0.035, ymax = 0.035)
  AT (0, 0, 0) RELATIVE Center_of_rotation

/*COMPONENT Sample =  Incoherent(radius = 0.005,
  yheight = 0.05, focus_r = 0, pack = 1,
  target_x = 0, target_y = 0, target_z = 1)
  AT (0,0,0) RELATIVE Sample_location
  EXTEND
  %{
    if (!SCATTERED) ABSORB;
  %}*/

SPLIT COMPONENT Sample = PowderN(
    reflections = "Fe.laz", radius = 0.006,
    yheight = 0.05)
  AT (0,0,0) RELATIVE Sample_location
  EXTEND
  %{
    if (!SCATTERED) ABSORB;
  %}

/*
SPLIT COMPONENT Sample = Single_crystal(xwidth=0.01, yheight=0.01, zdepth=0.01,
mosaic = 5, reflections="YBaCuO.lau")
  AT (0,0,0) RELATIVE Sample_location
  EXTEND
  %{
    if (!SCATTERED) ABSORB;
  %}
*/

/*
COMPONENT PSD_After_Sample_location = PSD_monitor(
    nx = 100, ny = 100, filename = "PSD_After_Sample_location",
    xmin = -0.2, xmax = 0.2, ymin = -0.2, ymax = 0.2)
  AT (0, 0, 1) RELATIVE Sample_location
*/

/*
COMPONENT DiffBeamStop = Beamstop(
    xmin = -0.2,xmax = 0.2, ymin=-0.2, ymax=0.2)
  AT (0, 0, mono_to_cor+1.35) RELATIVE Prim_axes
*/


//**************************************************************************************************
//Detector Axes
COMPONENT Det_axis = Arm()
  AT (0, 0, 0) RELATIVE Center_of_rotation
  ROTATED (0, det_takeoff, 0) RELATIVE Center_of_rotation

COMPONENT Det_axis_2 = COPY(Det_axis)
  AT (0, 0, 0) RELATIVE Det_axis
  ROTATED (0, det_cover_angle, 0) RELATIVE Det_axis

COMPONENT Det_axis_3 = COPY(Det_axis)
  AT (0, 0, 0) RELATIVE Det_axis_2
  ROTATED (0, det_cover_angle, 0) RELATIVE Det_axis_2

COMPONENT Det_axis_4 = COPY(Det_axis)
  AT (0, 0, 0) RELATIVE Det_axis_3
  ROTATED (0, det_cover_angle, 0) RELATIVE Det_axis_3

//**************************************************************************************************
//ROC
//Single ROC: yheight=0.147, length=0.0863 , 45degrees total, 35 blades, radius at window=0.20066
//Cannot make three seperate collimators as it complains about redefinition of the widths and theta ??
COMPONENT RadColl = Collimator_radial(
    xwidth=0.006096, yheight=0.147, length=0.0863, theta_min=-15.625, theta_max=15.625*7, nchan=35*4, radius=0.20066, nslit=35*4)
AT (0, 0, 0) RELATIVE Det_axis

/*COMPONENT RadColl = Collimator_radial(
    xwidth=0.006096, yheight=0.147, length=0.0863, theta_min=-22.5, theta_max=22.5, nchan=0, radius=0.20066, nslit=35)
AT (0, 0, -0.0002) RELATIVE Det_axis*/


//**************************************************************************************************
//Detector

/*COMPONENT PSD_4pi = PSD_monitor_4PI(radius=0.7, nx=101, ny=51,
  filename="vanadium.psd")
  AT (0,0,0) RELATIVE Det_axis*/

COMPONENT PSD_Detector = PSD_monitor(
    nx = 330, ny = 15, filename = "PSD_Detector",
    xmin = -1*det_width/2.0, xmax = det_width/2.0, ymin = -1*det_height/2.0, ymax = det_height/2.0)
  AT (0, 0, cor_to_det) RELATIVE Det_axis
  GROUP Detectors

COMPONENT PSD_Detector_2 = COPY(PSD_Detector)(
    filename = "PSD_Detector_2")
  AT (0, 0, cor_to_det) RELATIVE Det_axis_2
  GROUP Detectors

COMPONENT PSD_Detector_3 = COPY(PSD_Detector)(
    filename = "PSD_Detector_3")
  AT (0, 0, cor_to_det) RELATIVE Det_axis_3
  GROUP Detectors

COMPONENT PSD_Detector_4 = COPY(PSD_Detector)(
    filename = "PSD_Detector_4")
  AT (0, 0, cor_to_det) RELATIVE Det_axis_4
  GROUP Detectors

/*COMPONENT Det_axis_multi = Arm(
    )
  AT (0, 0, 0) RELATIVE PREVIOUS(2)
  ROTATED (0, det_cover_angle, 0) RELATIVE PREVIOUS(2)

COMPONENT COPY(PSD_Detector) = COPY(PSD_Detector)
  AT (0, 0, cor_to_det) RELATIVE PREVIOUS
  JUMP Det_axis_multi ITERATE ndet */


/*COMPONENT Det_axis_multi = Arm(
    )
  AT (0, 0, 0) RELATIVE PREVIOUS
  ROTATED (0, det_cover_angle, 0) RELATIVE PREVIOUS

COMPONENT PSD_Decttor_multi = */



FINALLY
%{
%}
END