/********************************************************************************
* 
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
* 
* This file was written by McStasScript, which is a 
* python based McStas instrument generator written by 
* Mads Bertelsen in 2019 while employed at the 
* European Spallation Source Data Management and 
* Software Centre
* 
* Instrument SPHERES_POL
* 
* %Identification
* Written by: Python McStas Instrument Generator
* Date: 10:19:15 on April 15, 2025
* Origin: ESS DMSC
* %INSTRUMENT_SITE: Generated_instruments
* 
* 
* %Parameters
* 
* %End 
********************************************************************************/

DEFINE INSTRUMENT SPHERES_POL (
double source_num = 1, 
double inc_div_h = 3.6, 
double inc_div_v = 5.6, 
double inc_lambda0 = 6.2708, 
double inc_energy = 2.0803, 
double inc_lambda_fwhm = 12.5, 
double use_input_file = 0, 
string inc_file = "DNS_1E8.dat", 
double mosaic_v = 2.2, 
double mosaic_h = 5.7, 
double graphite_tilt = 2.5, 
double graphite_distance = 0.005, 
double r0 = 0.7, 
double speed = 225, 
double radius = 0.43, 
double pst_type = 0, 
double mono_distance = 2, 
double mono_width = 0.5, 
double mono_height = 0.25, 
double mono_dE = 0.032, 
double Vlength = 0.3, 
double V_nvs = 1, 
string shape_file = "shape.dat", 
double Glength = 7.2, 
double sample_num = 1, 
double Sample_radius = 0.012, 
double Sample_thickness = 0.001, 
double Sample_height = 0.045, 
double Sample_split = 5 
)

DECLARE 
%{
double gspeedoverR;
double gradius;
double ggraphite_tilt;
int gpst_type;
int guse_input_file;
int scattered_on_crystal1;
int scattered_on_crystal2;
int scattered_on_crystal3;
double E0;
char E_mono_opts[200];
int scattered_on_sample;
int scattered_on_analyzer = 0;
int scattered_largeanalyzer = 0;
int scattered_smallanalyzer = 0;
double xi_analyzer;
int scattered_twice = 0;
double vd_value;
%}

INITIALIZE 
%{
// Start of initialize for generated SPHERES_POL
	
   if (use_input_file)
      MPI_MASTER(printf("Using input file '%s'\n",inc_file));
      
   guse_input_file = use_input_file;
      
   E0 = 81.8042 / inc_lambda0 / inc_lambda0;
   gspeedoverR = speed / radius;
   gradius = radius;
   gpst_type = pst_type;
   ggraphite_tilt = graphite_tilt;
   
   sprintf(E_mono_opts,"energy limits=[%g %g] bins=100, noborder",E0-mono_dE,E0+mono_dE);
   
  if (mcgravitation) {
     MPI_MASTER(fprintf(stderr,"ERROR: This will not work with gravity on!\n"));
     exit(-1);
  }
%}

TRACE 
COMPONENT origin = Progress_bar()
WHEN (source_num == 1)
AT (0,0,-0.06) ABSOLUTE

COMPONENT source = Source_div(
 xwidth = 0.04, yheight = 0.04,
 focus_aw = 2*inc_div_h, focus_ah = 2*inc_div_v,
 lambda0 = inc_lambda0, dlambda = inc_lambda0*inc_lambda_fwhm/100*FWHM2RMS,
 gauss = 1, flux = 1)
WHEN (source_num == 1)
AT (0,0,-0.05) ABSOLUTE

COMPONENT CompInstanceName = Virtual_input(
 filename = inc_file)
WHEN (source_num == 2)
AT (0, 0, -0.001) RELATIVE PREVIOUS

COMPONENT makeSpin = Spin_random()
AT (0,0,0) ABSOLUTE

COMPONENT hdiv_source = Monitor_nD(
 xwidth = 0.5, yheight = 0.5,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, hdiv limits=[-5 5] bins=200 border")
AT (0,0,0) ABSOLUTE

COMPONENT vdiv_source = Monitor_nD(
 xwidth = 0.5, yheight = 0.5,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, vdiv limits=[-5 5] bins=200 border")
AT (0,0,0) ABSOLUTE

COMPONENT l_source = Monitor_nD(
 options = "previous lambda limits=[4 8] bins=200")
AT (0,0,0) ABSOLUTE

COMPONENT guide_bsel = Guide_gravity(
 w1 = 0.06, h1 = 0.12,
 l = 9.4, m = 2,
 mleft = 2, mright = 2,
 mtop = 2, mbottom = 2)
AT (0, 0, 0) ABSOLUTE

COMPONENT hdiv_bsel = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, hdiv limits=[-5 5] bins=200 border")
AT (0,0,9.4) RELATIVE PREVIOUS

COMPONENT vdiv_bsel = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, vdiv limits=[-5 5] bins=200 border")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT l_bsel = Monitor_nD(
 options = "previous lambda limits=[4.6 8] bins=200")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT vel_sel = MySelector(
 xmin = -0.06, xmax = 0.06,
 ymin = -0.03, ymax = 0.03,
 length = 0.25, nslit = 64,
 d = 0.0004, radius = 0.16,
 alpha = 48.24, nu = 341)
AT (0, 0,  9.4+0.14) RELATIVE guide_bsel
ROTATED (0,0,-90) ABSOLUTE

COMPONENT after_sel = Arm()
AT (0,0,0) RELATIVE vel_sel
ROTATED (0,0,90) RELATIVE vel_sel

COMPONENT hdiv_asel = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, hdiv limits=[-5 5] bins=200 border")
AT (0,0,0.25) RELATIVE PREVIOUS

COMPONENT vdiv_asel = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, vdiv limits=[-5 5] bins=200 border")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT l_asel = Monitor_nD(
 options = "previous lambda limits=[4.6 8] bins=200")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT E_asel = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1.3 3.2] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT focus_guide = myGuide_tapering2(
 file = shape_file, l = Glength,
 R0 = 0.99, Qcx = 0.021,
 Qcy = 0.021, alphax = 2.95,
 alphay = 2.95, W = 0.0015,
 curvature = 0, curvature_v = 0,
 verbose = 0)
AT (0,0,0.251) RELATIVE after_sel

COMPONENT chopper_window = Slit(
 xmin = -0.03, xmax = 0.03,
 ymin = -0.03, ymax = 0.03)
AT (0,0,Glength+0.01) RELATIVE focus_guide

COMPONENT PSD_guideEnd = Monitor_nD(
 xwidth = 0.05, yheight = 0.05,
 bins = 200, restore_neutron = 1,
 options = "x y")
AT (0,0,Glength+0.01) RELATIVE focus_guide

COMPONENT hdiv_PST = Monitor_nD(
 xwidth = 0.15, yheight = 0.15,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, hdiv limits=[-5 5] bins=200 border")
AT (0,0,0.12) RELATIVE PREVIOUS

COMPONENT vdiv_PST = Monitor_nD(
 xwidth = 0.15, yheight = 0.15,
 restore_neutron = 1, options = "energy limits=[1 3.5] bins=200, vdiv limits=[-5 5] bins=200 border")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT PSD_PST = Monitor_nD(
 xwidth = 0.1, yheight = 0.1,
 bins = 200, restore_neutron = 1,
 options = "x y")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT PST_arm = Arm()
AT (0,0,0.12) RELATIVE chopper_window
ROTATED (0,20.84,0) ABSOLUTE

COMPONENT move_in_crystal1 = Arm()
WHEN (gpst_type==0)
AT (0,0,0) RELATIVE PST_arm
EXTEND %{
{
PROP_Z0;
vx -= gspeedoverR * (gradius + y);
vy += gspeedoverR * x;
}
%}

COMPONENT PST0_crystal1 = Monochromator_flat(
 zmin = -0.03, zmax = 0.03,
 ymin = -0.016, ymax = 0.016,
 mosaich = 60*mosaic_v, mosaicv = 60*mosaic_v,
 r0 = r0, DM = 3.355)
WHEN (gpst_type==0)
AT (0,0,1e-5) RELATIVE move_in_crystal1
ROTATED (0,90,0) RELATIVE move_in_crystal1
EXTEND %{
{
scattered_on_crystal1 = SCATTERED;
}
%}

COMPONENT move_out_crystal1 = Arm()
WHEN (gpst_type==0)
AT (0,0,0) RELATIVE move_in_crystal1
EXTEND %{
{
PROP_Z0;
vx += gspeedoverR * (gradius + y);
vy -= gspeedoverR * x;
}
%}

COMPONENT move_in_crystal2 = Arm()
WHEN (gpst_type==0 && !scattered_on_crystal1)
AT (0,0,graphite_distance) RELATIVE PST_arm
ROTATED (0,graphite_tilt,0) RELATIVE PST_arm
EXTEND %{
{
PROP_Z0;
vx -= gspeedoverR * (gradius + y) * cos(ggraphite_tilt*DEG2RAD);
vz -= gspeedoverR * (gradius + y) * sin(ggraphite_tilt*DEG2RAD);
vy += gspeedoverR * x;
}
%}

COMPONENT PST0_crystal2 = Monochromator_flat(
 zmin = -0.03, zmax = 0.03,
 ymin = -0.016, ymax = 0.016,
 mosaich = 60*mosaic_v, mosaicv = 60*mosaic_v,
 r0 = r0, DM = 3.355)
WHEN (gpst_type==0 && !scattered_on_crystal1)
AT (0,0,1e-5) RELATIVE move_in_crystal2
ROTATED (0,90,0) RELATIVE move_in_crystal2
EXTEND %{
{
scattered_on_crystal1 = SCATTERED;
}
{
scattered_on_crystal2 = SCATTERED;
}
%}

COMPONENT move_out_crystal2 = Arm()
WHEN (gpst_type==0 && !scattered_on_crystal1)
AT (0,0,0) RELATIVE move_in_crystal2
EXTEND %{
{
vx += gspeedoverR * (gradius + y) * cos(ggraphite_tilt*DEG2RAD);
vz += gspeedoverR * (gradius + y) * sin(ggraphite_tilt*DEG2RAD);
vy -= gspeedoverR * x;
}
%}

COMPONENT move_in_crystal3 = Arm()
WHEN (gpst_type==0 && !scattered_on_crystal1 && !scattered_on_crystal2)
AT (0,0,3*graphite_distance) RELATIVE PST_arm
ROTATED (0,-graphite_tilt,0) RELATIVE PST_arm
EXTEND %{
{
PROP_Z0;
vx -= gspeedoverR * (gradius + y) * cos(-ggraphite_tilt*DEG2RAD);
vz -= gspeedoverR * (gradius + y) * sin(-ggraphite_tilt*DEG2RAD);
vy += gspeedoverR * x;
}
%}

COMPONENT PST0_crystal3 = Monochromator_flat(
 zmin = -0.03, zmax = 0.03,
 ymin = -0.016, ymax = 0.016,
 mosaich = 60*mosaic_v, mosaicv = 60*mosaic_v,
 r0 = r0, DM = 3.355)
WHEN (gpst_type==0 && !scattered_on_crystal1 && !scattered_on_crystal2)
AT (0,0,1e-5) RELATIVE move_in_crystal3
ROTATED (0,90,0) RELATIVE move_in_crystal3
EXTEND %{
{
scattered_on_crystal1 = SCATTERED;
}
{
scattered_on_crystal3 = SCATTERED;
}
%}

COMPONENT move_out_crystal3 = Arm()
WHEN (gpst_type==0 && !scattered_on_crystal1 && !scattered_on_crystal2)
AT (0,0,0) RELATIVE move_in_crystal3
EXTEND %{
{
vx += gspeedoverR * (gradius + y) * cos(-ggraphite_tilt*DEG2RAD);
vz += gspeedoverR * (gradius + y) * sin(-ggraphite_tilt*DEG2RAD);
vy -= gspeedoverR * x;
}
%}

COMPONENT PST_out = Arm()
AT (0,0,0) RELATIVE PST_arm
ROTATED (0,180.0+20.84,0) RELATIVE PST_arm

COMPONENT Monochromator_position = Arm()
AT (0,0,mono_distance) RELATIVE PST_out

COMPONENT PSD_out = Monitor_nD(
 xwidth = 2, yheight = 1,
 bins = 200, restore_neutron = 1,
 options = "x y")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT l_out = Monitor_nD(
 xwidth = 1, yheight = 1,
 restore_neutron = 1, options = "lambda limits=[5.6 6.8] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT E_out = Monitor_nD(
 xwidth = 1, yheight = 1,
 restore_neutron = 1, options = "energy limits=[1.9 2.3] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT v_out = Monitor_nD(
 xwidth = 1, yheight = 1,
 restore_neutron = 1, options = "user1 limits = [600 650] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT vdiv_E_out = Monitor_nD(
 xwidth = 0.7, yheight = 0.7,
 restore_neutron = 1, options = "energy limits=[1.7 2.5] bins=200, vdiv limits=[-11 11] bins=200 border")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT hdiv_E_out = Monitor_nD(
 xwidth = 0.7, yheight = 0.7,
 restore_neutron = 1, options = "energy limits=[1.7 2.5] bins=200, hdiv limits=[-20 20] bins=200")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT E_mono = Monitor_nD(
 xwidth = mono_width, yheight = mono_height,
 restore_neutron = 1, options = E_mono_opts)
AT (0,0,0) RELATIVE Monochromator_position

COMPONENT beamstop_fc = Beamstop(
 xmin = -0.03, xmax = 0.03,
 ymin = -0.02, ymax = 0.02)
AT (0,0,0.2) RELATIVE chopper_window

COMPONENT Monochromator = Perfect_Crystal(
 width = 0.5, height = 0.25,
 centerfocus = 0, radius = 2.180,
 lambda = 6.2708, R0 = 0.98,
 sigma = 0.12e-3, ismirror = 0,
 speed = 4.7, amplitude = 0.075,
 smartphase = 1, exclusive = 1,
 verbose = 1)
AT (0,0,mono_distance) RELATIVE PST_out

COMPONENT Sample_position = Arm()
AT (0,0,-0.36) RELATIVE PST_out
ROTATED (0,180.0,0) RELATIVE PST_out

COMPONENT PSD_sample = Monitor_nD(
 xwidth = 0.1, yheight = 0.1,
 bins = 200, restore_neutron = 1,
 options = "x y")
AT (0,0,0) RELATIVE Sample_position

COMPONENT PSD_sample_lim = Monitor_nD(
 xwidth = 0.03, yheight = 0.04,
 bins = 200, restore_neutron = 1,
 options = "x y")
AT (0,0,0) RELATIVE Sample_position

COMPONENT l_sample = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "lambda limits=[6.1 6.4] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT E_sample = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[2.078 2.083] bins=100")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT vdiv_E_sample = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1.7 2.5] bins=200, vdiv limits=[-11 11] bins=200 border")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT hdiv_E_sample = Monitor_nD(
 xwidth = 0.25, yheight = 0.25,
 restore_neutron = 1, options = "energy limits=[1.7 2.5] bins=200, hdiv limits=[-20 20] bins=200")
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT MeanPol_Det = MeanPolLambda_monitor(
 xwidth = 0.1, yheight = 0.1,
 nL = 100, restore_neutron = 1,
 filename = "MeanPol_LA_center.csv", mx = 0,
 my = 1, mz = 0,
 Lmin = 6.2, Lmax = 6.35)
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT divlambda_Down = DivLambda_monitor_Down(
 nL = 100, nh = 100,
 filename = "divlambda_Down.csv", xwidth = 0.1,
 yheight = 0.1, maxdiv_h = 20,
 Lmin = 6.2, Lmax = 6.35,
 restore_neutron = 1)
AT (0,0,0) RELATIVE PREVIOUS

COMPONENT div_lambda_hor_Up_back = DivLambda_monitor_Up(
 nL = 100, nh = 100,
 filename = "divlambda_Up.csv", xwidth = 0.1,
 yheight = 0.1, maxdiv_h = 20,
 Lmin = 6.2, Lmax = 6.35,
 restore_neutron = 1)
AT (0,0,0) RELATIVE PREVIOUS

FINALLY 
%{
// Start of finally for generated SPHERES_POL
%}

END