Fix feeder not beeing fully illuminated and add VS and MF divergence monitors

analysis/compress_h5.sh

@@ -0,0 +1,9 @@
+ #!/bin/bash
+# compress McStas hdf5 files to increase performance and save disk space
+
+for fi in $*
+do
+  echo $fi
+  h5repack -i $fi -o $fi.c -f /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:GZIP=5 -l /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:CHUNK=3072x8
+  mv $fi.c $fi
+done

analysis/mcstas_reader.py

@@ -11,7 +11,7 @@ from numpy import *
 try:
     import h5py
 except ImportError:
-    print "h5py not found, modern NeXuS format will not be readable."
+    print("h5py not found, modern NeXuS format will not be readable.")
 
 try:
   from IPython import display #@UnusedImport
@@ -49,7 +49,7 @@ class McSim(object):
       elif os.path.exists(os.path.join(path, 'mccode.sim')):
         self._init_old(os.path.join(path, 'mccode.sim'))
       else:
-        raise IOError, "Can't locate mccode.h5 of mccode.sim file in %s"%path
+        raise IOError("Can't locate mccode.h5 of mccode.sim file in %s"%path)
 
   def _init_hdf(self, path):
     self.hdf=h5py.File(path, 'r')
@@ -61,11 +61,11 @@ class McSim(object):
     self.data_loader=DataLoaderOld(self.info, os.path.dirname(path))
 
   def keys(self):
-    return self.info['data'].keys()
+    return list(self.info['data'].keys())
 
   def monitors(self):
     output=[]
-    for val in self.info['data'].values():
+    for val in list(self.info['data'].values()):
       xy=(val['xvar'], val['yvar'])
       if not xy in output:
         output.append(xy)
@@ -74,7 +74,7 @@ class McSim(object):
   
   def plot(self, monitors=None):
     graphs={}
-    for key, val in self.info['data'].items():
+    for key, val in list(self.info['data'].items()):
       xy=(val['xvar'], val['yvar'])
       if monitors is not None and xy!=monitors:
         continue
@@ -102,12 +102,12 @@ class McSim(object):
   def __getitem__(self, item):
     if item in self._data:
       return self._data[item]
-    elif item in self.keys():
+    elif item in list(self.keys()):
       data=self.data_loader.load_item(item)
       self._data[item]=data
       return data
     else:
-      raise KeyError, "Can't find dataset %s"%item
+      raise KeyError("Can't find dataset %s"%item)
 
 class HeaderFile(object):
   '''
@@ -124,7 +124,7 @@ class HeaderFile(object):
     data=open(path, 'r').read()
     data_lines=data.splitlines()
     if not data.startswith('McStas simulation description file'):
-      raise IOError, 'Not a valid McStas description file.'
+      raise IOError('Not a valid McStas description file.')
     self._data['start_time']=data_lines[1].split(':', 1)[1].strip()
     self._data['program_name']=data_lines[2].split(':', 1)[1].strip()
     self._data['data']=self.get_data(data)
@@ -167,13 +167,13 @@ class DataLoaderOld(object):
     fname=os.path.join(self.root, item_info['filename'])
     x_col=item_info['xvar']
     y_col=item_info['yvar']
-    if x_col=='Li' and y_col=='p': # Detector_nD     
+    if x_col in ['Li', 'List'] and y_col=='p': # Detector_nD     
       cols=item_info['variables'].split() 
       data=loadtxt(fname, dtype={'names': cols, 'formats': ['f4']*len(cols)})
       return TofData(data, item_info)
     else:
       raw=loadtxt(fname)
-      data=raw[:len(raw)/3]
+      data=raw[:len(raw)//3]
       return Dataset(data, item_info)
 
   def load_item_1d(self, item):
@@ -192,11 +192,11 @@ class DataLoaderHDF(object):
     self.hdf=hdf['entry1']
     self.info={}
     self.info['data']={}
-    for item in self.hdf['data'].keys():
+    for item in list(self.hdf['data'].keys()):
       node=self.hdf['data/'+item]
       info={}
-      for key, value in node.attrs.items():
-        info[key.strip()]=value.strip()
+      for key, value in list(node.attrs.items()):
+        info[key.strip()]=value.strip().decode('ascii')
       if info['filename'][-4:]=='.dat' or '_list.' in info['filename']:
         self.info['data'][info['component']]=info
       else:
@@ -211,7 +211,7 @@ class DataLoaderHDF(object):
     node=self.hdf[item_info['datapath']]
     x_col=item_info['xvar']
     y_col=item_info['yvar']
-    if x_col=='Li' and y_col=='p': # Detector_nD
+    if x_col in ['Li', 'List'] and y_col=='p': # Detector_nD
       cols=item_info['variables'].split()
       evds=node['events']
       if len(evds)<=MAX_EVTS_BATCH:
@@ -253,7 +253,7 @@ class Dataset1D(object):
       import pylab
       ax=pylab.gca()
 
-    limits=map(float, self.info['xlimits'].split())
+    limits=list(map(float, self.info['xlimits'].split()))
     x=linspace(limits[0], limits[1], len(self.data))
 
     ax.errorbar(x, self.data, yerr=self.errors)
@@ -291,12 +291,12 @@ class Dataset(object):
       import pylab
       ax=pylab.gca()
 
-    limits=map(float, self.info['xylimits'].split())
+    limits=list(map(float, self.info['xylimits'].split()))
 
     if log:
-      ax.imshow(self.data, extent=limits, aspect='auto', norm=LogNorm())
+      ax.imshow(self.data, origin='lower', extent=limits, aspect='auto', norm=LogNorm())
     else:
-      ax.imshow(self.data, extent=limits, aspect='auto')
+      ax.imshow(self.data, origin='lower', extent=limits, aspect='auto')
     ax.set_xlabel(self.info['xlabel'])
     ax.set_ylabel(self.info['ylabel'])
     ax.set_title(self.info['component'])
@@ -342,7 +342,7 @@ class TofData(Dataset):
     if fltr is None:
       I, x=histogram(columns[col], bins=bins, weights=w)
     else:
-      if isinstance(fltr, basestring):
+      if isinstance(fltr, str):
         fltr=eval(fltr, globals(), columns)
       I, x=histogram(columns[col][fltr], bins=bins, weights=w[fltr])
     return x, I
@@ -367,7 +367,7 @@ class TofData(Dataset):
       I, y, x=histogram2d(columns[ycol], columns[xcol],
                         bins=bins, weights=self.data['p'])
     else:
-      if isinstance(fltr, basestring):
+      if isinstance(fltr, str):
         fltr=eval(fltr, globals(), columns)
       I, y, x=histogram2d(columns[ycol][fltr], columns[xcol][fltr],
                         bins=bins, weights=columns['p'][fltr])

simulation/Estia_baseline.instr

@@ -53,8 +53,10 @@ DEFINE INSTRUMENT ESS_reflectometer_Estia
 int operationmode = 0, double over_illumination = 0.0, double theta_resolution = 0.04,
 double lambda_min      =  3.75, double lambda_start    =  3.0, double lambda_end      = 12.0,
 int enable_chopper = 0, int enable_gravity=0, int enable_windows=1,
- int enable_polarizer = 0, int enable_analyzer = 0,
-double source_power = 2, double frame_usage     =   0.98
+ int enable_polarizer = 0, int enable_analyzer = 0, 
+ double pol1_start=0.3, double pol1_angle=1.66, double pol2_start=0.3, double pol2_angle=1.66,
+double source_power = 2,
+double selene1_foot1y =0.0, double selene1_foot2y = 0.0
 )
 
 DECLARE
@@ -70,8 +72,14 @@ double iscs_z=0.0547095; // ISCS in McStas coordinates using mod-view-opt.instr
 double iscs_rot_y=0.4; // ISCS is at TCS-35.6 degree, McStas at -36 degree
 double iscs_rot_x=0.7; // downward tilt of Estia axis
 
-// Selene 1 (Will follow later with detailed parameters)
+// Selene 1 geometry parameters (optics parameters in Estia_selene1.instr)
+double selene1_foot1   =  4.20; // distance of first foot to VS focus
+double selene1_foot2   =  7.00; // distance of second foot to VS focus
+double selene1_center;
+double selene1_shift;
+double selene1_rot;
 
+int p_int=0; // a flag that gets incremented if a polarizer mirror scatters
 
 // Selene 2
 
@@ -97,9 +105,6 @@ double slit_distance   =   1.775; //m, distance slit to sample
 double selene_theta    =    1.25;
 double velocity_max           ;     // m/s  neutron velocity of lambda_min
 
-/* flags */
-int PP_small         = 1 ;   // plot all PSD with small area
-int PP_large         = 1 ;   // plot all PSD with large area
 
 double analyzer_max_height = 0.01; // Maximum sample height to be covered by te polarization analyzers
 double analyzer1_start  = 0.65; // Distance to sample to start the first analyzer
@@ -115,12 +120,12 @@ INITIALIZE
 // chopper coupling together
 velocity_max= 3.956034E3 / lambda_min; // h/m_n over lambda - ((3.9..e-7m^2/s))/(1e-10m)
 chopper_phase=360.0*(chopper_pos*chopper_freq/velocity_max+(pulse_zero-opening_time)*chopper_freq);
-chopper_open =360.0*(chopper_pos/(total_length+detector_arm))*frame_usage;
+chopper_open = 98.0;
 
 
 /* print out some calculated parameter for checking purposes */
 printf(" Chopper phase = %.1f deg\n", chopper_phase);
-printf(" Chopper open  = %.1f deg\n", chopper_open);
+//printf(" Chopper open  = %.1f deg\n", chopper_open);
 
 
 dist_ana_vfocus=analyzer_max_height/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
@@ -130,7 +135,11 @@ Theta2_analyzer1=atan((dist_ana_vfocus+analyzer1_start+analyzer1_length)/dist_an
 Theta1_analyzer2=atan((dist_ana_vfocus+analyzer2_start)/dist_ana_vfocus*analyzer_max_height/2.0/analyzer2_start)*180.0/PI;
 Theta2_analyzer2=atan((dist_ana_vfocus+analyzer2_start+analyzer2_length)/dist_ana_vfocus*analyzer_max_height/2.0/(analyzer2_start+analyzer2_length))*180.0/PI;
 
-
+selene1_center=(selene1_foot1+selene1_foot2)/2.;
+selene1_shift=(selene1_foot1y+selene1_foot2y)/1000.0;
+selene1_rot=atan((selene1_foot2y-selene1_foot1y)/(selene1_foot2-selene1_foot1)/1000.0)*180.0/PI;
+//printf(" Selene 1 rotation  = %.4f deg\n", selene1_rot);
+//printf(" Selene 1 shift  = %.1f mm\n", selene1_shift*1e3);
 
 %}
 TRACE
@@ -157,10 +166,12 @@ COMPONENT ISCS = Arm() // rotate around y-axis (slight downward tilt)
 COMPONENT arm_feeder = Arm()
     AT (0, 0, 0) RELATIVE ISCS
     ROTATED (0, 0, 0) RELATIVE ISCS
+COMPONENT arm_selene1_center = Arm()
+    AT (selene1_shift, 0, 2*NBOA_c+selene1_center) RELATIVE ISCS
+    ROTATED (0, selene1_rot, 0) RELATIVE ISCS
 // Axes starting at focus of feeder parallel to c-axis of Selene guides
 COMPONENT arm_selene1 = Arm()
-    AT (0, 0, 2*NBOA_c) RELATIVE ISCS
-    ROTATED (0, 0, 0) RELATIVE ISCS
+    AT (0, 0, -selene1_center) RELATIVE arm_selene1_center
 COMPONENT arm_selene2 = Arm()
     AT (0, 0, 2*NBOA_c+2*selene_c) RELATIVE ISCS
     ROTATED (0, 0, 0) RELATIVE ISCS
@@ -194,20 +205,11 @@ COMPONENT moderator = ESS_butterfly(
     sector = "E", beamline = 2, yheight = 0.03, cold_frac = 0.9,
         focus_xw = E02_01_01_Cu_in_xmax-E02_01_01_Cu_in_xmin+0.002, 
         focus_yh = E02_01_01_Cu_in_yheight+0.002,
-    target_index=4, Lmin = lambda_start, Lmax = lambda_end,
+    target_index=3, Lmin = lambda_start, Lmax = lambda_end,
     n_pulses = 1+enable_chopper, acc_power=source_power)
     AT (0, 0, 0) RELATIVE origin 
 
 
-// monitor used for brilliance transfer calculations
-COMPONENT tof_source = Monitor_nD(
-    filename = "tof_source",
-    options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
-    xwidth=0.01, yheight = 0.01)
-    WHEN sample==4
-    AT (0, 0, 0.03) RELATIVE ISCS
-        ROTATED (0,-selene_theta,0) RELATIVE ISCS
-
 /***************************************
  * Geometry of neutron feeder separate *
  ***************************************/
@@ -216,11 +218,14 @@ COMPONENT tof_source = Monitor_nD(
 /**************************************************** 
  * Beam manipulation area around the virtual source *
  ****************************************************/
-/* Absorber to reduce beam to needed size an for shielding purposes (FeNi in CAD model) */
-COMPONENT FeNi_in = Slit(xwidth=0.028, yheight=0.078)
-    AT (0, 0, -0.800) RELATIVE arm_virtual_source_beam
+/* Absorber to reduce beam to needed size an for shielding purposes (CPC1 in CAD model) */
+COMPONENT CPC1_in = Slit(xwidth=0.0303, yheight=0.0792)
+    AT (0, 0, -0.890) RELATIVE arm_virtual_source_beam
 
-COMPONENT FeNi_out = Slit(xwidth=0.013, yheight=0.038)
+COMPONENT CPC1_monitor = Slit(xwidth=0.016344, yheight=0.05547)
+    AT (0, 0, -0.3573) RELATIVE arm_virtual_source_beam
+
+COMPONENT CPC1_out = Slit(xwidth=0.013, yheight=0.038)
     AT (0, 0, -0.220) RELATIVE arm_virtual_source_beam
 
 
@@ -231,30 +236,17 @@ COMPONENT chopper = DiskChopper(radius=chopper_diameter/2.0, yheight=0.02,
     WHEN enable_chopper==1
     AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
     
-    // Beam monitoring before VS
-    COMPONENT DP_before_virtual_source = DivPos_monitor(
-        filename = "DP_before_virtual_source" ,
-        nh = 200, ndiv = 200,
-        xwidth = 0.2, yheight = 0.2,
-        maxdiv_h = 6,
-        restore_neutron = 1)
-        AT (0, 0, -0.022) RELATIVE arm_virtual_source_beam
 
-    COMPONENT PP_small_before_virtual_source = PSD_monitor(
-        filename = "PP_small_before_virtual_source",
-        nx = 200, ny = 200,
-        xwidth = 0.1, yheight = 0.1,
-        restore_neutron = 1)
-        WHEN PP_small ==1
-        AT (0, 0, -0.021) RELATIVE arm_virtual_source_beam
+COMPONENT vs_divergence_h = DivPos_monitor(nh=21, ndiv=41, filename="vs_hordiv",
+          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0,
+		  restore_neutron=1)
+		  AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
 
-    COMPONENT PP_large_before_virtual_source = PSD_monitor(
-        filename = "PP_large_before_virtual_source",
-        nx = 200, ny = 200,
-        xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
-        restore_neutron = 1)
-        WHEN PP_large ==1
-        AT (0, 0, 2*NBOA_c-0.02) RELATIVE ISCS
+COMPONENT vs_divergence_v = DivPos_monitor(nh=21, ndiv=41, filename="vs_verdiv",
+          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0,
+		  restore_neutron=1)
+		  AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
+		  ROTATED (0,0,90) RELATIVE arm_selene1
 
 
 /* The actual virtual source mask, two L-shaped absorbers (first top-right) */
@@ -263,66 +255,26 @@ COMPONENT virtual_source_TR = Slit(
     WHEN sample!=4
     AT (-over_illumination, 0, -0.5*sample_length) RELATIVE arm_virtual_source
 
-// window to absorb neutron not passing through heavy collimation
-COMPONENT virtual_source_HC = Slit(
-    xmin = -0.015, xmax = 0.015,
-    ymin = -0.015, ymax = 0.015)
+// window to cut down to defined size for test setting
+COMPONENT virtual_source_HC = Slit(xwidth=sample_length, yheight=sample_height)
     WHEN sample==4
     AT (0, 0, 0) RELATIVE arm_virtual_source
-    
+//     
 /* The actual virtual source mask, two L-shaped absorbers (second bottom-left) */
 COMPONENT virtual_source_BL = Slit(
     xmin = -1.0, xmax = 0.0, ymin = -sample_height/2-over_illumination*5, ymax = 1.0)
     WHEN sample!=4
     AT (over_illumination, 0, 0.5*sample_length) RELATIVE arm_virtual_source
 
-    
-    // Beam monitoring after VS
-    COMPONENT tof_virtual_source = Monitor_nD(
-        filename = "tof_virtual_source",
-        options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
-        xwidth=sample_length, yheight = sample_height)
-        WHEN sample==4
-        AT (0, 0, 0.019) RELATIVE arm_virtual_source_beam
-                
-    COMPONENT DL_behind_virtual_source = DivLambda_monitor(
-        filename = "DL_behind_virtual_source",
-        nL = 200, nh = 200,
-        xwidth = 0.02, yheight = 0.06,
-        maxdiv_h = 3,
-        Lmin = lambda_start, Lmax = lambda_end, 
-        restore_neutron = 1)
-        AT (0, 0, 0.02) RELATIVE arm_virtual_source_beam
-
-    COMPONENT DP_behind_virtual_source = DivPos_monitor(
-        filename = "DP_behind_virtual_source" ,
-        nh = 200, ndiv = 200,
-        xwidth = 0.2, yheight = 0.2,
-        maxdiv_h = 3,
-        restore_neutron = 1)
-        AT (0, 0, 0.021) RELATIVE arm_virtual_source_beam
-
-    COMPONENT PP_small_behind_virtual_source = PSD_monitor(
-        filename = "PP_small_behind_virtual_source",
-        nx = 200, ny = 200, 
-        xwidth = 0.02, yheight = 0.05,
-        restore_neutron = 1)
-        WHEN PP_small ==1
-        AT (0, 0, 0.022) RELATIVE arm_virtual_source_beam
-
-    COMPONENT PP_large_behind_virtual_source = PSD_monitor(
-        filename = "PP_large_behind_virtual_source",
-        nx = 200, ny = 200,
-        xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
-        restore_neutron = 1)
-        WHEN PP_large ==1
-        AT (0, 0, 2*NBOA_c+0.023) RELATIVE ISCS
 
 /*************************************
  * Geometry of Selene guide separate *
  *************************************/
-%include "Estia_selene.instr"
+%include "Estia_selene1.instr"
 
+%include "Estia_mf.instr"
+
+%include "Estia_selene2.instr"
 
 /**********************************
  * Optical components within cave *
@@ -343,60 +295,14 @@ COMPONENT ac_slit = Slit(
 /***************
  * Sample area *
  ***************/
-    // Beam monitoring before sample
-    COMPONENT PP_small_sample_focus = PSD_monitor(
-        filename = "PP_small_sample_focus",
-        nx = 200, ny = 200, 
-        xwidth = 0.02, yheight = 0.05,
-        restore_neutron = 1)
-        WHEN PP_small ==1
-        AT (0, 0, -0.001) RELATIVE arm_sample
-
-    COMPONENT DL_sample_focus = DivLambda_monitor(
-        filename = "DL_sample_focus",
-        nL = 200, nh = 200,
-        xwidth = 0.02, yheight = 0.06,
-        maxdiv_h = 3,
-        Lmin = lambda_start, Lmax = lambda_end, 
-        restore_neutron = 1)
-        AT (0, 0, -0.001) RELATIVE arm_sample
-
-    /* monitor at sample position for footprint */
-    COMPONENT PL_sample = Monitor_nD(
-        filename = "PL_sample",
-        options = "lambda limits=[0 35] bins=350 x limits=[-0.02 0.02] bins 160",
-        xwidth = 0.04, yheight = 0.02, restore_neutron=1)
-        WHEN sample!=4
-        AT (0, 0, -0.001) RELATIVE arm_sample
-        ROTATED (0, -90, 0) RELATIVE arm_sample
-
-    COMPONENT PSD_sample = PSD_monitor(
-        filename = "PSD_sample", 
-        nx = 100, ny = 100,
-        xwidth = 4*sample_length, yheight = 4*sample_height,
-        restore_neutron = 1)
-        WHEN sample!=4
-        AT (0, 0, -0.001) RELATIVE arm_sample
-        ROTATED (0, -90, 0) RELATIVE arm_sample
-
     COMPONENT tof_sample = Monitor_nD(
-        filename = "tof_sample",
-        options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
-        xwidth=sample_length, yheight = sample_height)
+        filename = "tof_sample", user1=p_int,
+        options = "x limits=[-0.025 0.025] bins=1000 y limits=[-0.025 0.025] bins=1000 xdiv limits=[-0.75 0.75] bins=150 ydiv limits=[-2.0 2.0] bins=400 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sy limits=[-1 1] bins=2000 user1 list all",
+        xwidth=0.05, yheight = 0.05)
         WHEN sample==4
-        AT (0, 0, -0.0005) RELATIVE arm_sample_beam
-        ROTATED (0, 0, 0) RELATIVE arm_sample_beam
-                
-    /* monitor for full beam at sample position ( n/cm²/s ) */
-    COMPONENT PP_small_sample = PSD_monitor(
-        filename = "PSD_sample_perp",
-        nx = 200, ny = 200, 
-        xwidth = 0.05, yheight = 0.05,
-        restore_neutron = 1)
-        WHEN sample==4 
         AT (0, 0, 0) RELATIVE arm_sample_beam
         ROTATED (0, 0, 0) RELATIVE arm_sample_beam
-
+                
 
 /* NiTi multilayer sample */
 COMPONENT sample = Mirror(
@@ -448,14 +354,14 @@ COMPONENT arm_analyzer2 = Arm()
     ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector
 
 /* polarization analyser */
-COMPONENT analyzer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer1_start, length=analyzer1_length,
+COMPONENT analyzer1 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer1_start, length=analyzer1_length,
                                 delta_theta=(Theta1_analyzer1+Theta2_analyzer1+0.05)*PI/180.0, h2=0.14, h1=0.05, abs_ref=0, 
                                 m_u=5.5, m_d=0.45, both_coated=0, alpha=2.3, W = 0.0014)
         WHEN enable_analyzer
         AT (0, 0.0, analyzer1_start) RELATIVE arm_analyzer
         ROTATED (0,0,0.0) RELATIVE arm_analyzer
 
-COMPONENT analyzer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer2_start, length=analyzer2_length,
+COMPONENT analyzer2 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer2_start, length=analyzer2_length,
                                 delta_theta=(Theta1_analyzer2+Theta2_analyzer2+0.05)*PI/180.0, h2=0.2, h1=0.05, abs_out=0, 
                                 m_u=5.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
         WHEN enable_analyzer==2
@@ -463,31 +369,13 @@ COMPONENT analyzer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, ref
         ROTATED (0,0,0.0) RELATIVE arm_analyzer
 
 /* detector */
-    // Beam monitoring at detector position
-    COMPONENT PP_large_detector = PSD_monitor(
-        filename = "PP_large_detector",
-        nx = 200, ny = 200, 
-        xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
-        restore_neutron = 1)
-        AT (0, 0, total_length+detector_arm-0.001) RELATIVE ISCS
-
-    // Detector with 0.5x1mm² resolution
-    COMPONENT PP_detector = PSD_monitor(
-        filename = "PP_detector",
-        nx = 1000, ny = 500,
-        xwidth = 0.5, yheight = 0.5,
-        restore_neutron = 1)
-        AT (0, 0, detector_arm) RELATIVE arm_detector
-
     COMPONENT tof_detector = Monitor_nD(
-        filename = "tof_detector",
-        options = "x limits=[-0.25 0.25] bins=1000 y limits=[-0.5 0.5] bins=1000 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sx limits=[-1 1] bins=1000 sy limits=[-1 1] bins=1000, list all",
+        filename = "tof_detector", user1=p_int,
+        options = "x limits=[-0.25 0.25] bins=1000 y limits=[-0.5 0.5] bins=1000 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sx sy sz user1, list all",
         xwidth = 0.5, yheight = 1.0)
         WHEN sample!=4
         AT (0, 0, detector_arm+0.00001) RELATIVE arm_detector
-    /*COMPONENT detector = PSD_TOF_monitor(xwidth=0.5, yheight=0.5,
-            nx=1000, ny=250, tmin=35000, tmax=115000, nt=160, filename="detector")
-        AT (0, 0, detector_arm) RELATIVE arm_detector
+        ROTATED (0, 0, 0) RELATIVE arm_detector
 
 /***********************************************************************/
 

simulation/Estia_mf.instr

@@ -0,0 +1,246 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: Estia_mf
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
+* Date: 01. 03. 2018
+* Origin: PSI
+* Release: McStas 2.4.1
+* Version: 1.0
+* %INSTRUMENT_SITE: ESS (E02)
+*
+* Estia is a vertical sample, focusing reflectometer for small sample
+*
+* %Description
+* These are the components near the middle focus between the Selene guide 1 and 2.
+* This file is not intended for stand alone use but is included in the Estia instrument model. 
+* The removable components, however, allow for a stand alone use for debugging purpose.
+*
+* %Parameters
+* enable_gravity: [0/1] Use gravity in elliptical guide model.
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
+
+DECLARE
+%{
+
+// Polarizer parameters used for the detailed design
+double pol1_xstart = 1.007; // [m] up-stream of focus point
+double pol1_ystart = 0.0064; // [m] (relative to c-axis)
+double pol1_xend = 0.365; // [m] up-stream of focus point
+double pol1_yend = 0.013; // [m]
+double pol1_gamma = 1.65; // [deg] optimal incident angle for beams passing through focus point
+double pol1_m = 5.08;
+double pol1_lin_xstart = 1.007; //[m]
+double pol1_lin_ystart = 0.0064; //[m]
+double pol1_lin_xend = 0.812; //[m]
+double pol1_lin_yend = -0.0014; //[m]
+double pol1_lin_m = 5.5;
+
+
+double pol2_xstart = 1.286; // [m] up-stream of focus point
+double pol2_ystart = 0.009; // [m] (relative to c-axis)
+double pol2_xend = 0.500; // [m] up-stream of focus point
+double pol2_yend = 0.0175; // [m]
+double pol2_gamma = 1.70; // [deg] optimal incident angle for beams passing through focus point
+double pol2_m = 5.08;
+double pol2_lin_xstart = 1.286; //[m]
+double pol2_lin_ystart = 0.009; //[m]
+double pol2_lin_xend = 1.093; //[m]
+double pol2_lin_yend = 0.0014; //[m]
+double pol2_lin_m = 5.08;
+
+double pol_hclose = 0.12; // [m] height of entrance
+double pol_hfar = 0.12; // [m] height of exit
+
+// parameters for polarizer component calculated from the above
+double pol1_length, pol2_length, pol1_lin_length, pol2_lin_length;
+double Theta_pol1, Theta_pol2, Theta_rot1, Theta_rot2;
+double Theta1_pol1, Theta1_pol2, Theta2_pol1, Theta2_pol2;
+double pol1_lin_rot, pol2_lin_rot;
+
+%}
+
+INITIALIZE
+%{
+
+pol1_length = pol1_xstart-pol1_xend;
+pol2_length = pol2_xstart-pol2_xend;
+pol1_lin_length = pol1_lin_xstart-pol1_lin_xend;
+pol2_lin_length = pol2_lin_xstart-pol2_lin_xend;
+
+Theta1_pol1=atan2(pol1_ystart, pol1_xstart); // [rad] angle between c-axis and entrance point
+Theta2_pol1=atan2(pol1_yend, pol1_xend); // [rad] angle between c-axis and exit point
+Theta1_pol2=atan2(pol2_ystart, pol2_xstart); // [rad] angle between c-axis and entrance point
+Theta2_pol2=atan2(pol2_yend, pol2_xend); // [rad] angle between c-axis and exit point
+
+Theta_pol1=(Theta2_pol1-Theta1_pol1)*180.0/PI; // [deg] full covered divergence angle
+Theta_pol2=(Theta2_pol2-Theta1_pol2)*180.0/PI; // [deg]
+Theta_rot1=(Theta1_pol1+Theta2_pol1)/2.0*180.0/PI; // [deg] rotation of center of polarizer
+Theta_rot2=(Theta1_pol2+Theta2_pol2)/2.0*180.0/PI; // [deg]
+
+pol1_lin_rot=atan2(pol1_lin_ystart-pol1_lin_yend, 
+                   pol1_lin_xstart-pol1_lin_xend)*180.0/PI; // [deg] rotation angle of polarizer
+pol2_lin_rot=atan2(pol2_lin_ystart-pol2_lin_yend, 
+                   pol2_lin_xstart-pol2_lin_xend)*180.0/PI; // [deg] rotation angle of polarizer
+   
+
+
+printf(" Polarizer 1 angle      = %.2f deg\n", Theta_rot1);
+printf(" Polarizer 1 divergence = %.2f deg\n", Theta_pol1);
+printf(" Polarizer 1 distance   = %.1f mm\n", pol1_xstart*1000.0);
+printf(" Polarizer 1 length     = %.1f mm\n", pol1_length*1000.0);
+printf(" PolStraight 1 angle      = %.2f deg\n", pol1_lin_rot);
+
+printf(" Polarizer 2 angle      = %.2f deg\n", Theta_rot1);
+printf(" Polarizer 2 divergence = %.2f deg\n", Theta_pol2);
+printf(" Polarizer 2 distance   = %.1f mm\n", pol2_xstart*1000.0);
+printf(" Polarizer 2 length     = %.1f mm\n", pol2_length*1000.0);
+printf(" PolStraight 2 angle      = %.2f deg\n", pol2_lin_rot);
+
+%}
+TRACE
+
+REMOVABLE COMPONENT origin = Progress_bar()
+    AT (0,0,0) ABSOLUTE
+
+REMOVABLE COMPONENT ISCS = Arm()
+    AT (0, 0, 0) RELATIVE origin    
+    
+REMOVABLE COMPONENT arm_selene1 = Arm()
+    AT (0, 0, 0) RELATIVE ISCS    
+    
+REMOVABLE COMPONENT arm_selene2 = Arm()
+    AT (0, 0, 2*selene_c) RELATIVE ISCS    
+
+REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002, 
+                                       target_index=3, Emin=1, Emax=15)
+    AT (0,0,0) RELATIVE ISCS
+    ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
+
+
+/**************************************
+* Middle focus between Selene guides *
+**************************************/
+COMPONENT arm_polref = Arm()
+    AT (0, 0, 2*selene_c) RELATIVE arm_selene1
+    ROTATED (selene_theta, -selene_theta, 0) RELATIVE ISCS
+
+COMPONENT arm_polarizer = Arm()
+    AT (0, 0, 2*selene_c) RELATIVE arm_selene1
+    ROTATED (0, 0, 0) RELATIVE arm_selene1
+
+COMPONENT arm_pol1 = Arm()
+    AT (0, 0, 0) RELATIVE arm_polarizer
+    ROTATED (0, -Theta_rot1, 0) RELATIVE arm_polarizer
+
+COMPONENT arm_pol2 = Arm()
+    AT (0, 0, 0) RELATIVE arm_polarizer
+    ROTATED (0, -Theta_rot2, 0) RELATIVE arm_polarizer
+
+
+COMPONENT polarizer2_lin = Pol_mirror(zwidth = pol2_lin_length, yheight = pol_hfar, p_reflect=0,
+                                      rUpFunc=TableReflecFunc, rUpPar="SNAG_m5p5_Tup.ref", 
+									  rDownFunc=TableReflecFunc, rDownPar="SNAG_m5p5_Tdown.ref", useTables=1)
+		WHEN (enable_polarizer>0)
+        AT ((pol2_lin_ystart+pol2_lin_yend)/2.0, 0, -(pol2_lin_xstart+pol2_lin_xend)/2.0) RELATIVE arm_polarizer
+        ROTATED (0,-pol2_lin_rot,0.0) RELATIVE arm_polarizer
+		GROUP polarizer2_set
+		EXTEND
+		%{
+			p_int=0; // reset for each trace
+			if (SCATTERED==2) {
+				p_int +=2;
+				}
+		%}
+
+COMPONENT polarizer2 = Polariser(lin=-pol2_xstart, length=pol2_length,
+                                 enable_ref=1, abs_ref=1, abs_out=1, both_coated=1,
+                                 d_substrate = 5e-4, T_loss=4.0e3,
+                                 m_u=pol2_m, m_d=0.6, m_residual=0.55,
+                                 alpha=2.3, W = 0.0014,  reflect_d=0, reflect_u=0, 
+                                 delta_theta=Theta_pol2*PI/180.0, 
+                                 h2=pol_hfar, h1=pol_hclose)
+		WHEN (enable_polarizer>0)
+        AT (0, 0, -pol2_xstart) RELATIVE arm_pol2
+        ROTATED (0,0,90.0) RELATIVE arm_pol2
+		GROUP polarizer2_set
+		EXTEND
+		%{
+			if (SCATTERED) {
+				ALLOW_BACKPROP;    
+				PROP_Z0;
+				p_int +=1;
+				} else ABSORB;
+		%}
+
+COMPONENT replacement_pol2_set = Slit(xwidth=0.2, yheight=0.2)
+	WHEN (enable_polarizer==0)
+    AT (0, 0, -0.2) RELATIVE arm_polarizer
+	GROUP polarizer2_set
+
+COMPONENT polarizer1_lin = Pol_mirror(zwidth = pol1_lin_length, yheight = pol_hfar, p_reflect=0,
+                                      rUpFunc=TableReflecFunc, rUpPar="SNAG_m5p5_Tup.ref", 
+									  rDownFunc=TableReflecFunc, rDownPar="SNAG_m5p5_Tdown.ref", useTables=1)
+		WHEN (enable_polarizer>0)
+        AT ((pol1_lin_ystart+pol1_lin_yend)/2.0, 0, -(pol1_lin_xstart+pol1_lin_xend)/2.0) RELATIVE arm_polarizer
+        ROTATED (0,-pol1_lin_rot,0.0) RELATIVE arm_polarizer
+		GROUP polarizer1_set
+		EXTEND
+		%{
+			if (SCATTERED==2) {
+				p_int +=8;
+			}
+		%}
+
+COMPONENT polarizer1 = Polariser(lin=-pol1_xstart, length=pol1_length,
+                                 enable_ref=1, abs_ref=1, abs_out=1, both_coated=1,
+                                 d_substrate = 5e-4, T_loss=4.0e3,
+                                 m_u=pol1_m, m_d=0.6, m_residual=0.55,
+                                 alpha=2.3, W = 0.0014,  reflect_d=0, reflect_u=0, 
+                                 delta_theta=Theta_pol1*PI/180.0, 
+                                 h2=pol_hfar, h1=pol_hclose)
+		WHEN (enable_polarizer>0)
+        AT (0, 0, -pol1_xstart) RELATIVE arm_pol1
+        ROTATED (0,0,90.0) RELATIVE arm_pol1
+		GROUP polarizer1_set
+		EXTEND
+		%{
+			if (SCATTERED) {
+				p_int +=4;
+			}
+		%}
+
+COMPONENT replacement_pol1_set = Slit(xwidth=0.2, yheight=0.2)
+	WHEN (enable_polarizer==0)
+    AT (0, 0, -0.1) RELATIVE arm_polarizer
+	GROUP polarizer1_set
+
+
+
+COMPONENT mf_divergence_h = DivPos_monitor(nh=21, ndiv=41, filename="mf_hordiv",
+          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0)
+		  AT (0, 0, 0) RELATIVE arm_polarizer
+
+COMPONENT mf_divergence_v = DivPos_monitor(nh=21, ndiv=41, filename="mf_verdiv",
+          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0)
+		  AT (0, 0, 0) RELATIVE arm_polarizer
+		  ROTATED (0,0,90) RELATIVE arm_polarizer
+
+
+
+COMPONENT virtual_flipper = Pol_SF_ideal(ny=1, xwidth=0.1, yheight=0.1, zdepth=0.001)
+        WHEN (enable_polarizer>2)
+        AT (0, 0, 0) RELATIVE arm_polarizer
+    
+
+FINALLY
+%{
+%}
+
+END

simulation/Estia_monitor.instr

@@ -0,0 +1,204 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: ESS_reflectometer_Estia
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
+* Date: 01. 03. 2018
+* Origin: PSI
+* Release: McStas 2.4.1
+* Version: 1.0
+* %INSTRUMENT_SITE: ESS (E02)
+*
+* Estia is a vertical sample, focusing reflectometer for small sample
+*
+* %Description
+* The instrument consists of a two part elliptical feeder that focuses onto 
+* a slit mask that defines the sample footprint, called virtual source.
+* The virtual source is imaged on the sample position with a Selene 
+* type neutron guide (two ellipses).
+* This version of the instrument, used in the ESS butterfly moderator, 
+* has two vertical Selene guide systems that share the same focal points.
+* The two guides are only implemented in the feeder as the selene mirrors
+* are part of the upgrade program.
+*
+* %Parameters
+* omegaa: [deg]  sample rotation omega
+* sample: flag to switch between     (0) sample Ni/Ti-multilayer
+*                                    (1) reference ( R(q_z) = 1 for all q_z )
+*                                    (2) Ni-film on silicon
+*                                    (3) Natural SiO2 on silicon
+*                                    (4) monitor instead of sample (normal to beam)
+* sample_length: [m] Size of sample in beam direction, also controls virtual source opening
+* sample_height: [m] Size of sample in vertical direction, also controls virtual source opening
+* 
+* operationmode: operation mode      (0) high-intensity specular reflectivity
+*                                    (1) almost conventional TOF
+* over_illumination: [m] Extra opening of virtual source compared to sample size
+* theta_resolution: Delta theta / theta  adjusted with the slit (operationmode=1)
+* 
+* lambda_start: [A] Beginning of simulated wavelength range
+* lambda_end:   [A] End of simulated wavelength range
+* enable_gravity: [0/1] Use gravity in elliptical guide model.
+* enable_chopper: [0/1/2/3] Activate chopper component if !=0. Numbers larger than
+*                           1 define the number of pulses per chopper opening.
+*                           (pulses skipped=enable_chopper-1)
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT ESS_reflectometer_Estia
+(double lambda_start    =  3.0, double lambda_end      = 12.0,
+int enable_gravity=1, int enable_windows=1, direct_beam=0,
+double source_power = 5, foil_thickness=0.00001)
+
+DECLARE
+%{
+/* Geometrical parameters from CAD model of Estia (ESS-0050413)
+* TCS coordinate and directional rotation first focus point 
+* refered to as focus_moderator_y_rot
+*/
+//TCS position of ISCS: (110,-105,137)
+double iscs_x=0.0199717;
+double iscs_y=0.0;
+double iscs_z=0.0547095; // ISCS in McStas coordinates using mod-view-opt.instr
+double iscs_rot_y=0.4; // ISCS is at TCS-35.6 degree, McStas at -36 degree
+double iscs_rot_x=0.7; // downward tilt of Estia axis
+
+// Selene 1 geometry parameters (optics parameters in Estia_selene1.instr)
+double selene1_foot1   =  4.20; // distance of first foot to VS focus
+double selene1_foot2   =  7.00; // distance of second foot to VS focus
+double selene1_center;
+double selene1_shift;
+double selene1_rot;
+
+// Selene 2
+
+
+    
+/* general instrument geometry parameters */
+double total_length    = 35.0 ;     // m    distance moderator-sample (2*c_feeder+4*c_Selene)
+double detector_arm    =  4.0 ;     // m    distance sample-detector
+
+
+/* chopper parameters */
+double chopper_diameter = 0.7;
+double chopper_pos     = 10.895;     // m  distance source-chopper
+double chopper_phase          ;       // deg phase between pulse and chopper opening
+double chopper_open           ;       // deg of opening angle in the chopper
+double pulse_zero      =   0.00175;   // ms  intensity weighted average time of emitted neutron pulse
+double opening_time    =   0.0015; // ms time to reach full intensity on detector, used to adjust phase to get full intensity at beginning of selected band
+double chopper_freq    =  14.0 ;     // Hz  chopper frequency
+
+double slit_distance   =   1.775; //m, distance slit to sample
+
+/* derived quantities */
+double selene_theta    =    1.25;
+
+%}
+
+INITIALIZE
+%{
+
+%}
+TRACE
+
+/********************************************************************************
+ * Initial geometry coordinate axes for important components of the simulation. *
+ ********************************************************************************/
+COMPONENT origin = Progress_bar()
+    AT (0,0,0) ABSOLUTE
+    
+/* The ISCS is the instrument coordinate system with directions where y points upwards and z lies on the instrument axes. */
+COMPONENT ISCS_rot1 = Arm() // position correctly and rotate around z-axis (x points in beam direction)
+    AT (iscs_x,iscs_y,iscs_z) RELATIVE origin
+    ROTATED (0,iscs_rot_y,0) RELATIVE origin
+COMPONENT ISCS = Arm() // rotate around y-axis (slight downward tilt)
+    AT (0,0,0) RELATIVE ISCS_rot1
+    ROTATED (iscs_rot_x,0,0) RELATIVE ISCS_rot1
+
+    
+/***********************
+* Instrument Skeleton *
+***********************/
+// Axes system parallel to the c-axis of the feeder
+COMPONENT arm_feeder = Arm()
+    AT (0, 0, 0) RELATIVE ISCS
+    ROTATED (0, 0, 0) RELATIVE ISCS
+COMPONENT arm_selene1_center = Arm()
+    AT (selene1_shift, 0, 2*NBOA_c+selene1_center) RELATIVE ISCS
+    ROTATED (0, selene1_rot, 0) RELATIVE ISCS
+// Axes starting at focus of feeder parallel to c-axis of Selene guides
+COMPONENT arm_selene1 = Arm()
+    AT (0, 0, -selene1_center) RELATIVE arm_selene1_center
+// Axes of the center of usable beam passing through the virtual source
+COMPONENT arm_virtual_source_beam = Arm()
+    AT (0, 0, 0) RELATIVE arm_selene1
+    ROTATED (0, selene_theta, 0) RELATIVE ISCS
+// 90 degree to main beam is the monitors detector position
+COMPONENT arm_monitor = Arm()
+    AT (0, 0, -0.320) RELATIVE arm_virtual_source_beam
+    ROTATED (-90, 0, 0) RELATIVE arm_virtual_source_beam
+
+
+
+/**********
+ * Source *
+ **********/
+COMPONENT moderator = ESS_butterfly(
+    sector = "E", beamline = 2, yheight = 0.03, cold_frac = 0.9,
+        focus_xw = E02_01_01_Cu_in_xmax-E02_01_01_Cu_in_xmin+0.002, 
+        focus_yh = E02_01_01_Cu_in_yheight+0.002,
+    target_index=4, Lmin = lambda_start, Lmax = lambda_end,
+    n_pulses = 1, acc_power=source_power)
+    AT (0, 0, 0) RELATIVE origin 
+
+
+/***************************************
+ * Geometry of neutron feeder separate *
+ ***************************************/
+%include "Estia_feeder.instr"
+
+/**************************************************** 
+ * Beam manipulation area around the virtual source *
+ ****************************************************/
+/* Absorber to reduce beam to needed size an for shielding purposes (CPC1 in CAD model) */
+COMPONENT CPC1_in = Slit(xwidth=0.0303, yheight=0.0792)
+    AT (0, 0, -0.890) RELATIVE arm_virtual_source_beam
+
+COMPONENT CPC1_monitor = Slit(xwidth=0.016344, yheight=0.05547)
+    AT (0, 0, -0.3573) RELATIVE arm_virtual_source_beam
+
+// As TOF detector is rectangular, use focus_r size to limit to actual
+// detector size of 0.5'' diameter cylinder
+COMPONENT Vanadium_Foil = Incoherent(focus_r=0.00635, p_interact=0.9,
+                                     xwidth=0.018, yheight=0.055, zdepth=foil_thickness,
+                                     sigma_abs=5.08, sigma_inc=5.08, Vc=13.827,
+                                     target_index=1)
+    WHEN direct_beam<2
+    AT (0, 0, -0.320) RELATIVE arm_virtual_source_beam
+    ROTATED (34.3, 0, 0) RELATIVE arm_virtual_source_beam
+
+COMPONENT tof_monitor = TOFLambda_monitor(
+    filename = "monitor",
+    tmin=0, tmax=120000, nt=1200,
+    Lmin=0,Lmax=35,nL=350,
+    xwidth = 0.025, yheight = 0.025)
+    WHEN direct_beam==0
+    AT (0, 0, 0.1) RELATIVE arm_monitor
+
+COMPONENT tof_direct = TOFLambda_monitor(
+    filename = "VS",
+    tmin=0, tmax=120000, nt=1200,
+    Lmin=0,Lmax=35,nL=350,
+    xwidth = 0.05, yheight = 0.05)
+    WHEN direct_beam>0
+    AT (0, 0, 0.08) RELATIVE arm_virtual_source_beam
+
+
+FINALLY
+%{
+%}
+
+END

simulation/Estia_selene.instr

@@ -1,240 +0,0 @@
-/*******************************************************************************
-*         McStas instrument definition URL=http://www.mcstas.org
-*
-* Instrument: Estia_selene
-*
-* %Identification
-* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
-* Date: 01. 03. 2018
-* Origin: PSI
-* Release: McStas 2.4.1
-* Version: 1.0
-* %INSTRUMENT_SITE: ESS (E02)
-*
-* Estia is a vertical sample, focusing reflectometer for small sample
-*
-* %Description
-* This is the beam extraction part of the instrument, it is not intended for use but is
-* included in the Estia instrument model. The removable components, however,
-* allow for a stand alone use for debugging purpose.
-*
-* %Parameters
-* enable_gravity: [0/1] Use gravity in elliptical guide model.
-*
-* %End
-*******************************************************************************/
-
-DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
-
-DECLARE
-%{
-    
-// Selene 1 (Will follow later with detailed parameters)
-
-
-// Selene 2
-
-    
-/* selene parameters */
-double selene_xi       =  0.6 ;     //      length of neutron guide over 2*c of ellipses, xi
-double selene_b        =  0.1047;     // m selene short half axis c
-double selene_c        =  6.0000;     // m selene focus half distance c
-double selene_coating  =  3.8;        //   m-value of selene coating
-
-
-/* other variables */
-double selene_entry ;
-double selene_length ;
-double selene_distance ;
-
-double polarizer_max_width = 0.0015; // Maximum sample height to be covered by te polarization analyzers
-double polarizer_start  = 0.3; // Distance to sample to start the first analyzer
-double polarizer_length = 0.6; // length of first analyzer
-double Theta1_polarizer, Theta2_polarizer, dist_pol_vfocus; // quantities calculated out of values above and 1.5 degree covered divergence
-
-%}
-
-INITIALIZE
-%{
-/* derived parameters for the selene guides */
-selene_length   = 2.0 * selene_c * selene_xi;
-selene_distance = ( 1.0 - selene_xi ) * selene_c;
-selene_entry    = selene_b * sqrt( 1.0 - pow(selene_xi,2.0) );
-
-dist_pol_vfocus=polarizer_max_width/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
-
-Theta1_polarizer=atan((dist_pol_vfocus+polarizer_start)/dist_pol_vfocus*polarizer_max_width/2.0/polarizer_start)*180.0/PI;
-Theta2_polarizer=atan((dist_pol_vfocus+polarizer_start+polarizer_length)/dist_pol_vfocus*polarizer_max_width/2.0/(polarizer_start+polarizer_length))*180.0/PI;
-
-%}
-TRACE
-
-REMOVABLE COMPONENT origin = Progress_bar()
-    AT (0,0,0) ABSOLUTE
-
-REMOVABLE COMPONENT ISCS = Arm()
-    AT (0, 0, 0) RELATIVE origin    
-    
-REMOVABLE COMPONENT arm_selene1 = Arm()
-    AT (0, 0, 0) RELATIVE ISCS    
-    
-REMOVABLE COMPONENT arm_selene2 = Arm()
-    AT (0, 0, 2*selene_c) RELATIVE ISCS    
-
-REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002, 
-                                       target_index=3, Emin=1, Emax=15)
-    AT (0,0,0) RELATIVE ISCS
-    ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
-
-/**************************
-* Selene 1 neutron guide *
-**************************/
-/* Absorber to cut direct view beam (Copper in CAD model) */
-COMPONENT slit_before_selene_guide_1 = Slit(
-    xmin = 0, xmax = selene_entry+0.001,
-    ymin = 0, ymax = selene_entry+0.005)
-    AT (0, 0, selene_distance-0.01) RELATIVE arm_selene1
-
-COMPONENT block_before_selene_guide_1 = Absorber(
-    xmin=-1, xmax=1, 
-    ymin=-selene_entry/4.0+0.005, 
-    ymax=selene_entry/4.0-0.005,
-    zmin=0.0, zmax=0.001)
-    AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
-
-/* Selene 1 elliptic guide */
-COMPONENT selene_guide_1 = Elliptic_guide_gravity(
-    l=selene_length, dimensionsAt = "mid",
-    linyh = selene_distance, loutyh= selene_distance,
-    linxw = selene_distance, loutxw= selene_distance,
-    xwidth=selene_b*2, yheight=selene_b*2,
-    mright=0, mleft=selene_coating, mtop=selene_coating, mbottom=0,
-    enableGravity=enable_gravity) 
-    AT (0, 0, selene_distance) RELATIVE arm_selene1
-
-/* Absorber to cut direct view beam (Copper in CAD model) */    
-COMPONENT slit_after_selene_guide_1 = Slit(
-    xmin = 0, xmax = selene_entry+0.005,
-    ymin = 0, ymax = selene_entry+0.01)
-    AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene1
-
-COMPONENT block_after_selene_guide_1 = Absorber(
-    xmin=-1, xmax=1, 
-    ymin=-selene_entry/4.0+0.015, 
-    ymax=selene_entry/4.0-0.015,
-    zmin=0.0, zmax=0.001)
-    AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene1
-        
-
-/**************************************
-* Middle focus between Selene guides *
-**************************************/
-COMPONENT arm_polarizer = Arm()
-    AT (0, 0, 2*selene_c) RELATIVE arm_selene1
-    ROTATED (selene_theta, -selene_theta, 0) RELATIVE ISCS
-
-COMPONENT polarizer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=-(polarizer_start+polarizer_length), length=polarizer_length,
-                                 delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
-        WHEN enable_polarizer
-        AT (0, 0, -(polarizer_start+polarizer_length)) RELATIVE arm_polarizer
-        ROTATED (0,0,90.0) RELATIVE arm_polarizer
-
-    COMPONENT polmon_mfocus=PolLambda_monitor(Lmin=lambda_start, Lmax=lambda_end, nL=30, npol=101, my=1, filename="polmon_mfocus")
-        WHEN enable_polarizer
-            AT (0, 0, 2*selene_c) RELATIVE arm_selene1
-            ROTATED (0.0,0.0,0.0) RELATIVE arm_selene1
-
-  
-    COMPONENT PP_small_before_middle_focus = PSD_monitor(
-        filename = "PP_small_before_middle_focus",
-        nx = 200, ny = 200, 
-        xwidth = 0.1, yheight = 0.1,
-        restore_neutron = 1)
-        AT (0, 0, 2*selene_c) RELATIVE arm_selene1
-
-/* The proximal polariser comes next */
-COMPONENT polarizer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=polarizer_start, length=polarizer_length,
-                                 delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
-        WHEN (enable_polarizer>1)
-        AT (0, 0, polarizer_start) RELATIVE arm_polarizer
-        ROTATED (0.0,0,90.0) RELATIVE arm_polarizer
-        
-
-
-/**************************
-* Selene 2 neutron guide *
-**************************/                
-/* Absorber to cut direct view beam (Copper in CAD model) */    
-COMPONENT slit_before_selene_guide_2 = Slit(
-        xmin = -selene_entry-0.005, xmax=0,
-        ymin = -selene_entry-0.01, ymax = 0)
-        AT (0, 0, selene_distance-0.002) RELATIVE arm_selene2
-
-COMPONENT block_before_selene_guide_2 = Absorber(
-    xmin=-1, xmax=1, 
-    ymin=-selene_entry/4.0+0.015, 
-    ymax=selene_entry/4.0-0.015,
-    zmin=0.0, zmax=0.001)
-    AT (0, 0, selene_distance-0.0015) RELATIVE arm_selene2
-
-    
-/* Selene 2 elliptic guide first half */
-COMPONENT selene_guide_21 = Elliptic_guide_gravity(
-    l=0.5*selene_length-0.001, dimensionsAt = "mid",
-    linyh = selene_distance, loutyh= selene_c+0.001,
-    linxw = selene_distance, loutxw= selene_c+0.001,
-    xwidth=selene_b*2, yheight=selene_b*2,
-    mright=selene_coating, mleft=0, mtop=0, mbottom=selene_coating,
-    enableGravity=enable_gravity) 
-    AT (0, 0, selene_distance) RELATIVE arm_selene2
-
-/* Absorber to cut direct view beam (Copper in CAD model) */    
-COMPONENT slit_within_selene_guide_2 = Slit(
-    xmin = -selene_b, xmax=-selene_b*0.45+0.005,
-    ymin = -selene_b, ymax=-selene_b*0.45+0.005)
-    AT (0, 0, 1.0*selene_c) RELATIVE arm_selene2
-
-/*COMPONENT block_within_selene_guide_2 = Absorber(
-    xmin =-1, xmax=1,
-    ymin =-selene_b*0.25+0.001, 
-    ymax = selene_b*0.25-0.001,
-    zmin=0.0, zmax=0.0001)
-    AT (0, 0, 3.0*selene_c+0.0001) RELATIVE arm_selene2*/
-
-/* Selene 2 elliptic guide first half */
-COMPONENT selene_guide_22 = Elliptic_guide_gravity(
-    l=0.5*selene_length-0.001, dimensionsAt = "mid",
-    linyh = selene_c+0.001, loutyh= selene_distance,
-    linxw = selene_c+0.001, loutxw= selene_distance,
-    xwidth=selene_b*2, yheight=selene_b*2,
-    mright=selene_coating, mleft=0, mtop=0, mbottom=selene_coating,
-    enableGravity=enable_gravity) 
-    AT (0, 0, 1.0*selene_c+0.001) RELATIVE arm_selene2
-
-/* Absorber to cut direct view beam (Copper in CAD model) */    
-COMPONENT slit_after_selene_guide_2 = Slit(
-    xmin = -selene_entry-0.001, xmax=0, xmax = 0.0,
-    ymin = -selene_entry-0.005, ymax = 0)
-    AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene2
-
-COMPONENT block_after_selene_guide_2 = Absorber(
-    xmin=-1, xmax=1, 
-    ymin=-selene_entry/4.0+0.015, 
-    ymax=selene_entry/4.0-0.015,
-    zmin=0.0, zmax=0.001)
-    AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene2
-
-REMOVABLE COMPONENT FocusMonitor = PSD_monitor(
-        filename = "Focus", 
-        nx = 100, ny = 100,
-        xwidth = 0.005, yheight = 0.005,
-        restore_neutron = 1)
-    AT (0, 0, 4*selene_c) RELATIVE ISCS
-    ROTATED (0, 1.25, 0) RELATIVE ISCS
-    
-
-FINALLY
-%{
-%}
-
-END

simulation/Estia_selene1.instr

@@ -0,0 +1,270 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: Estia_selene1
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
+* Date: 01. 03. 2018
+* Origin: PSI
+* Release: McStas 2.4.1
+* Version: 1.0
+* %INSTRUMENT_SITE: ESS (E02)
+*
+* Estia is a vertical sample, focusing reflectometer for small sample
+*
+* %Description
+* This is the Selene guide 1 part of the instrument, it is not intended for use but is
+* included in the Estia instrument model. The removable components, however,
+* allow for a stand alone use for debugging purpose.
+*
+* %Parameters
+* enable_gravity: [0/1] Use gravity in elliptical guide model.
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
+
+DECLARE
+%{
+    
+// Selene 1
+
+
+// Selene 2
+
+    
+/* selene parameters */
+double selene_xi       =  0.6 ;     //      length of neutron guide over 2*c of ellipses, xi
+double selene_b        =  0.1047;     // m selene short half axis c
+double selene_c        =  6.0000;     // m selene focus half distance c
+
+/* Coating along the elliptic guide to allow sharp wavelenth cut-off */
+double scoating_01     =  3.58;  // coating at Selene 1 start/ Selene 2 end
+double scoating_02     =  3.36;
+double scoating_03     =  3.19;
+double scoating_04     =  3.07;
+double scoating_05     =  2.98;
+double scoating_06     =  2.92;
+double scoating_07     =  2.88;
+double scoating_08     =  2.86;  // coating in center
+double scoating_09     =  2.86;  
+double scoating_10     =  2.87;  
+double scoating_11     =  2.91;  
+double scoating_12     =  2.98;  
+double scoating_13     =  3.07;  
+double scoating_14     =  3.20;  
+double scoating_15     =  3.38;  // coating at Selene 1 end/ Selene 2 start
+
+
+/* other variables */
+double selene_entry ;
+double selene_length ;
+double selene_distance ;
+
+double selene_segment ; 
+
+%}
+
+INITIALIZE
+%{
+/* derived parameters for the selene guides */
+selene_length   = 2.0 * selene_c * selene_xi;
+selene_distance = ( 1.0 - selene_xi ) * selene_c;
+selene_entry    = selene_b * sqrt( 1.0 - pow(selene_xi,2.0) );
+selene_segment  = selene_length / 15.;
+
+%}
+TRACE
+
+REMOVABLE COMPONENT origin = Progress_bar()
+    AT (0,0,0) ABSOLUTE
+
+REMOVABLE COMPONENT ISCS = Arm()
+    AT (0, 0, 0) RELATIVE origin    
+    
+REMOVABLE COMPONENT arm_selene1 = Arm()
+    AT (0, 0, 0) RELATIVE ISCS
+    ROTATED (0, 0, 0) RELATIVE ISCS
+    
+REMOVABLE COMPONENT arm_selene2 = Arm()
+    AT (0, 0, 2*selene_c) RELATIVE ISCS    
+
+REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002, 
+                                       target_index=3, Emin=1, Emax=15)
+    AT (0,0,0) RELATIVE ISCS
+    ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
+
+/**************************
+* Selene 1 neutron guide *
+**************************/
+/* Absorber to cut direct view beam (Copper in CAD model) */
+COMPONENT slit_before_selene_guide_1 = Slit(
+    xmin = 0, xmax = selene_entry+0.001,
+    ymin = -selene_entry-0.005, ymax = 0)
+    AT (0, 0, selene_distance-0.01) RELATIVE arm_selene1
+
+COMPONENT block_before_selene_guide_1 = Absorber(
+    xmin=-1, xmax=1, 
+    ymin=-selene_entry/4.0+0.005, 
+    ymax=selene_entry/4.0-0.005,
+    zmin=0.0, zmax=0.001)
+    AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
+
+/* Selene 1 elliptic guide */
+COMPONENT E02_03_01 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
+    linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_01, mbottom=scoating_01, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance) RELATIVE arm_selene1
+
+COMPONENT E02_03_02 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
+    linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_02, mbottom=scoating_02, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_03 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
+    linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_03, mbottom=scoating_03, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_04 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
+    linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_04, mbottom=scoating_04, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_05 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
+    linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_05, mbottom=scoating_05, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_06 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
+    linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_06, mbottom=scoating_06, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_07 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
+    linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_07, mbottom=scoating_07, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_08 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
+    linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_08, mbottom=scoating_08, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_09 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
+    linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_09, mbottom=scoating_09, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_10 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
+    linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_10, mbottom=scoating_10, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_11 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
+    linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_11, mbottom=scoating_11, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_12 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
+    linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_12, mbottom=scoating_12, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_13 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
+    linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_13, mbottom=scoating_13, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_14 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
+    linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_14, mbottom=scoating_14, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene1
+
+COMPONENT E02_03_15 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
+    linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=0, mleft=scoating_15, mbottom=scoating_15, mtop=0,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene1
+
+/* Absorber to cut direct view beam (Copper in CAD model) */    
+COMPONENT slit_after_selene_guide_1 = Slit(
+    xmin = 0, xmax = selene_entry+0.005,
+    ymin = -selene_entry-0.005, ymax = 0)
+    AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene1
+
+COMPONENT block_after_selene_guide_1 = Absorber(
+    xmin=-1, xmax=1, 
+    ymin=-selene_entry/4.0+0.015, 
+    ymax=selene_entry/4.0-0.015,
+    zmin=0.0, zmax=0.001)
+    AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene1
+        
+    
+
+FINALLY
+%{
+%}
+
+END

simulation/Estia_selene2.instr

@@ -0,0 +1,210 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: Estia_selene2
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
+* Date: 01. 03. 2018
+* Origin: PSI
+* Release: McStas 2.4.1
+* Version: 1.0
+* %INSTRUMENT_SITE: ESS (E02)
+*
+* Estia is a vertical sample, focusing reflectometer for small sample
+*
+* %Description
+* This is the second Selene guide part of the instrument, it is not intended for use but is
+* included in the Estia instrument model. This file does not work no its own, as it
+* requires the definitions specified in the Estia_selene1.instr DECLAR and INITIALIZE.
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT Estia_selene2()
+
+DECLARE
+%{
+    
+
+%}
+
+INITIALIZE
+%{
+
+%}
+TRACE
+
+
+/**************************
+* Selene 2 neutron guide *
+**************************/                
+/* Absorber to cut direct view beam (Bor-Al in CAD model) */
+COMPONENT slit_before_selene_guide_2 = Slit(
+        xmin = -selene_entry-0.005, xmax=0,
+        ymin = 0, ymax = selene_entry+0.01)
+        AT (0, 0, selene_distance-0.002) RELATIVE arm_selene2
+
+COMPONENT block_before_selene_guide_2 = Absorber(
+    xmin=-1, xmax=1, 
+    ymin=-selene_entry/4.0+0.015, 
+    ymax=selene_entry/4.0-0.015,
+    zmin=0.0, zmax=0.001)
+    AT (0, 0, selene_distance-0.0015) RELATIVE arm_selene2
+
+    
+/* Selene 2 elliptic guide first half */
+COMPONENT E02_04_01 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
+    linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_15, mleft=0, mbottom=0, mtop=scoating_15,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance) RELATIVE arm_selene2
+
+COMPONENT E02_04_02 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
+    linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_14, mleft=0, mbottom=0, mtop=scoating_14,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_03 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
+    linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_13, mleft=0, mbottom=0, mtop=scoating_13,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_04 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
+    linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_12, mleft=0, mbottom=0, mtop=scoating_12,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_05 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
+    linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_11, mleft=0, mbottom=0, mtop=scoating_11,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_06 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
+    linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_10, mleft=0, mbottom=0, mtop=scoating_10,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_07 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
+    linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_09, mleft=0, mbottom=0, mtop=scoating_09,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_08 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
+    linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_08, mleft=0, mbottom=0, mtop=scoating_08,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_09 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
+    linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_07, mleft=0, mbottom=0, mtop=scoating_07,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_10 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
+    linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_06, mleft=0, mbottom=0, mtop=scoating_06,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_11 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
+    linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_05, mleft=0, mbottom=0, mtop=scoating_05,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_12 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
+    linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_04, mleft=0, mbottom=0, mtop=scoating_04,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_13 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
+    linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_03, mleft=0, mbottom=0, mtop=scoating_03,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_14 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
+    linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_02, mleft=0, mbottom=0, mtop=scoating_02,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene2
+
+COMPONENT E02_04_15 = Elliptic_guide_gravity(
+    l=selene_segment, dimensionsAt = "mid",
+    linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
+    linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
+    xwidth=selene_b*2, yheight=selene_b*2,
+    mright=scoating_01, mleft=0, mbottom=0, mtop=scoating_01,
+    enableGravity=enable_gravity) 
+    AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene2
+
+/* Absorber to cut direct view beam (Bor-Al in CAD model) */
+COMPONENT slit_after_selene_guide_2 = Slit(
+    xmin = -selene_entry-0.001, xmax=0, xmax = 0.0,
+    ymin = 0, ymax = selene_entry+0.01)
+    AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene2
+
+COMPONENT block_after_selene_guide_2 = Absorber(
+    xmin=-1, xmax=1, 
+    ymin=-selene_entry/4.0+0.015, 
+    ymax=selene_entry/4.0-0.015,
+    zmin=0.0, zmax=0.001)
+    AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene2
+    
+
+FINALLY
+%{
+%}
+
+END

simulation/Pol_SF_ideal.comp

@@ -0,0 +1,91 @@
+/*******************************************************************************
+*
+* McStas, neutron ray-tracing package
+*         Copyright 1997-2002, All rights reserved
+*         Risoe National Laboratory, Roskilde, Denmark
+*         Institut Laue Langevin, Grenoble, France
+*
+* Written by: Erik B Knudsen
+* Date: May 2017
+* Version: $Revision$
+* Release: McStas 2.4
+* Origin: DTU Physics
+*
+* Component: Pol_RFSF_ideal
+*
+*
+* %I
+*
+* Ideal model of a spin flipper
+*
+* %D
+* This component simply mirrors the polarization vector of the neutron
+* ray in the plane through (0,0,0) with normal nx,ny,nz.
+* The flipper is surrounded by a perfectly absorbing box. Neutron rays not hitting
+* the box are left untouched.
+*
+* %P
+* Input parameters:
+* nx:    [ ] x-component of the normal vector of the flipping plane.
+* ny:    [ ] y-component of the normal vector of the flipping plane.
+* nz:    [ ] z-component of the normal vector of the flipping plane.
+* xwidth: [m] width of the spin flipper.
+* yheight: [m] height of the spin flipper.
+* zdepth: [m] length of the spin flipper.
+*
+*
+* %E
+*******************************************************************************/
+
+DEFINE COMPONENT Pol_SF_ideal
+DEFINITION PARAMETERS ()
+SETTING PARAMETERS (nx=0, ny=1, nz=0, xwidth=0.1, yheight=0.1, zdepth=0.1)
+OUTPUT PARAMETERS ()
+/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
+TRACE
+%{
+    int hit;
+    double t0,t1;
+    hit=box_intersect(&t0,&t1, x,y,z,vx,vy,vz, xwidth,yheight,zdepth);
+    if(hit){
+        PROP_DT(t0);
+        if(fabs(z- -zdepth*0.5)>DBL_EPSILON){
+            /*neutron must have hit the side walls*/
+            ABSORB;
+        }
+        /*move to center of box and flip*/
+        PROP_Z0;
+        SCATTER;
+        double s=scalar_prod(sx,sy,sz,nx,ny,nz);
+        if (s!=0){
+            sx-=2*s*nx;
+            sy-=2*s*ny;
+            sz-=2*s*nz;
+        }
+        PROP_DT((t1-t0)/2);/*propagate the remaining distance to the box exit*/
+        if(fabs(z-zdepth*0.5)>DBL_EPSILON){
+            ABSORB;
+        }
+    }
+
+
+%}
+
+MCDISPLAY
+%{
+    double dx=xwidth/16;
+    double dy=yheight/8;
+  
+  box(0,0,0,xwidth,yheight,zdepth);
+  
+  line(dx,-dy,0,dx,-dy+yheight/2.0,0);
+  line(-dx,dy,0,-dx,dy-yheight/2.0,0);
+  line(dx,-dy+yheight/2.0,0, dx+xwidth/16,-dy+yheight-yheight/16,0);
+  line(dx,-dy+yheight/2.0,0, dx-xwidth/16,-dy+yheight-yheight/16,0);
+  
+  line(-dx,dy-yheight/2.0,0, -dx+xwidth/16,dy-yheight+yheight/16,0);
+  line(-dx,dy-yheight/2.0,0, -dx-xwidth/16,dy-yheight+yheight/16,0);
+
+%}
+
+END

simulation/PolarizerTest.instr

@@ -0,0 +1,108 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: PolarizerTest
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch)
+* Date: 25. 07. 2019
+* Origin: PSI
+* Release: McStas 2.5
+* Version: 1.0
+*
+* Test for polarization parameters of curved transmission polarizer
+*
+* %Description
+* 
+*
+* %Parameters
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT PolarizerTest(m_up=5.08, m_down=0.6, R0=1.00,
+                                m_residual=0.55, Theta=1.5,
+                                lambda_min=2.0, lambda_max=32.0,
+                                gamma=1.66, focus_height=0.00001,
+                                phi=0.0, polarizer_start=0.2
+                               )
+
+DECLARE
+%{
+double lambda_mean,delta_lambda,polarizer_length,zero_position;
+double source_div, detector_height;
+
+int use_streight=0;
+  
+%}
+
+INITIALIZE
+%{
+lambda_mean=(lambda_min+lambda_max)/2.0;
+delta_lambda=(lambda_max-lambda_min)/2.0;
+
+polarizer_length=(exp(Theta/gamma)-1.0)*polarizer_start;
+zero_position=exp(0.5*Theta/gamma)*polarizer_start;
+
+source_div=0.1*Theta/180.*PI;
+detector_height=source_div/0.1*(polarizer_length+polarizer_start+0.5);
+
+%}
+
+TRACE
+
+COMPONENT origin = Progress_bar()
+    AT (0,0,0) ABSOLUTE
+
+COMPONENT focus  = Arm()
+    AT (0,0,0) RELATIVE origin
+    ROTATED (phi,0,0) RELATIVE origin
+
+COMPONENT Source = Source_simple(xwidth = 0.00001, yheight = focus_height, 
+                                 dist = 0.1, focus_xw = 0.0001, 
+                                 focus_yh = source_div, 
+                                 lambda0 = lambda_mean, dlambda = delta_lambda, 
+                                 gauss=0)
+        AT (0, 0, 0) RELATIVE origin
+
+COMPONENT FluxIn = L_monitor(xwidth = 0.5, yheight = 0.5, filename="fluxin", 
+          nL=601, Lmin=lambda_min, Lmax=lambda_max)
+         AT (0, 0, 0.101) RELATIVE origin
+
+COMPONENT polarizer = Polariser(enable_ref=1, abs_ref=0, 
+                                reflect_d=0, reflect_u=0, T_loss=4.0e3,
+                                m_substrate=0.469522, d_substrate = 5e-4, 
+                                m_u=m_up, m_d=m_down, R0=R0,
+                                m_residual=m_residual,
+                                lin=polarizer_start, length=polarizer_length,
+                                delta_theta=(Theta)*PI/180.0, h2=0.1, h1=0.05,
+                                both_coated=1, alpha=2.3, W = 0.0014)
+        WHEN use_streight==0
+        AT (0, -5e-4, polarizer_start) RELATIVE focus
+
+COMPONENT streight = Mirror(xwidth = 0.2, yheight = 0.2, 
+                            m=5.0, center=1, transmit=0)
+        WHEN use_streight==1
+        AT (0, 0, zero_position) RELATIVE origin
+        ROTATED (90.0-gamma,0,0) RELATIVE origin
+
+
+COMPONENT FluxOut = L_monitor(xwidth = 0.5, yheight = detector_height, 
+                              filename="fluxout", 
+          nL=601, Lmin=lambda_min, Lmax=lambda_max)
+        AT (0, 0, polarizer_length+polarizer_start+0.5) RELATIVE origin
+COMPONENT PolMon = MeanPolLambda_monitor(xwidth=0.5, yheight=detector_height, 
+                                         mx=-1, 
+                      nL=601, Lmin=lambda_min, Lmax=lambda_max, 
+                      filename="polarization")
+        AT (0, 0, polarizer_length+polarizer_start+0.5) RELATIVE origin
+
+
+    
+FINALLY
+%{
+%}
+
+END
+
+

simulation/SNAG_m5p5_Tdown.ref

@@ -0,0 +1,58 @@
+4.37173E-05	0.253838498
+0.001967279	0.548126562
+0.004240578	0.781912929
+0.006513878	0.831781609
+0.008656026	0.935568791
+0.010841891	0.996046699
+0.012984039	0.994254859
+0.015038752	0.780745054
+0.017443204	0.571166135
+0.019497917	0.421992388
+0.021771217	0.347878292
+0.023913365	0.314863976
+0.026317816	0.273648917
+0.028241378	0.248910367
+0.030514678	0.229613667
+0.032700543	0.210315136
+0.034711539	0.201079179
+0.03711599	0.179153501
+0.039301855	0.174572119
+0.041444003	0.161570874
+0.043629868	0.159330731
+0.045772016	0.158627627
+0.047957881	0.133760075
+0.050231181	0.128230852
+0.052373329	0.135460425
+0.054559194	0.125901964
+0.056701342	0.120064418
+0.058974642	0.122922579
+0.061160507	0.127084935
+0.063171503	0.12641248
+0.06548852	0.117637177
+0.067630668	0.105706827
+0.069903968	0.112415362
+0.072089833	0.119429932
+0.074231981	0.105676178
+0.076417846	0.107716872
+0.078428842	0.107129046
+0.080833293	0.10262909
+0.082888007	0.094164854
+0.085161306	0.092239896
+0.087303454	0.085655788
+0.089620471	0.081615114
+0.091631467	0.076984785
+0.093817332	0.076157254
+0.096090632	0.070972035
+0.09823278	0.060551277
+0.100287493	0.071366358
+0.102691945	0.065426801
+0.104834093	0.059045344
+0.107019958	0.056229269
+0.109162106	0.061290519
+0.111347971	0.051160701
+0.113490119	0.053836338
+0.115763418	0.05839439
+0.117949284	0.057514252
+0.120091431	0.047536765
+0.122364731	0.04877829
+0.300000000 0.000000000		

simulation/SNAG_m5p5_Tup.ref

@@ -0,0 +1,58 @@
+4.37173E-05	0.224586133
+0.001967279	0.518057252
+0.004240578	0.775672976
+0.006513878	0.825879049
+0.008656026	0.899985882
+0.010841891	0.967789786
+0.012984039	0.99196296
+0.015038752	0.992909754
+0.017443204	0.994177124
+0.019497917	0.991424846
+0.021771217	0.993887344
+0.023913365	0.998648926
+0.026317816	0.999361915
+0.028241378	0.998678679
+0.030514678	0.998109885
+0.032700543	0.997772638
+0.034711539	0.997392554
+0.03711599	0.999600524
+0.039301855	0.994779551
+0.041444003	0.996274647
+0.043629868	0.996149207
+0.045772016	0.998305774
+0.047957881	0.99658353
+0.050231181	0.995950077
+0.052373329	0.995563046
+0.054559194	0.995479284
+0.056701342	0.993347667
+0.058974642	0.990508095
+0.061160507	0.990476431
+0.063171503	0.988754593
+0.06548852	0.987959462
+0.067630668	0.987058572
+0.069903968	0.988939946
+0.072089833	0.986821236
+0.074231981	0.985149296
+0.076417846	0.987362592
+0.078428842	0.986507086
+0.080833293	0.98295765
+0.082888007	0.983313364
+0.085161306	0.983817788
+0.087303454	0.983554288
+0.089620471	0.979763465
+0.091631467	0.975348682
+0.093817332	0.975887375
+0.096090632	0.972236927
+0.09823278	0.965515882
+0.100287493	0.962852402
+0.102691945	0.956266464
+0.104834093	0.953509902
+0.107019958	0.945236222
+0.109162106	0.808306574
+0.111347971	0.186001557
+0.113490119	0.078929784
+0.115763418	0.051976893
+0.117949284	0.05338487
+0.120091431	0.04844847
+0.122364731	0.053399168
+0.300000000 0.000000000

simulation/compile_if_needed.sh

@@ -1,9 +1,13 @@
 #!/bin/bash
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/afs/psi.ch/project/sinq/sl6-64/mcstas2.4/mcstas/2.4/libs
+
 
 if [ Estia_baseline.instr -nt Estia_baseline.out ] || [ ! -f Estia_baseline.out ] \
    || [ Estia_feeder.instr -nt Estia_baseline.out ] \
-   || [ Estia_selene.instr -nt Estia_baseline.out ]; then
+   || [ Estia_selene1.instr -nt Estia_baseline.out ] \
+   || [ Estia_mf.instr -nt Estia_baseline.out ] \
+   || [ Estia_selene2.instr -nt Estia_baseline.out ]; then
 	rm Estia_baseline.c Estia_baseline.out
-	mcstas -o Estia_baseline.c Estia_baseline.instr
-	mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
+	mcstas -o Estia_baseline.c Estia_baseline.instr --trace
+	mpicc -O2 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
 fi

simulation/hostnames

@@ -0,0 +1,3 @@
+mcc05:72
+mcc06:72
+

simulation/pol_sweep_analyze.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+'''
+Run a simulation through various values of polarizer 1/2 distance to focus
+and generate a datafile with analyzed values for polarization and transmission
+for the short and long wavelength area.
+'''
+
+import sys, os
+from numpy import *
+sys.path.append('../analysis')
+from mcstas_reader import McSim
+
+Iref=1.0
+L=linspace(1.95, 32.05, 302)
+Lc=(L[:-1]+L[1:])/2.
+
+gamma1=1.8
+gamma2=1.66
+command_s_ref="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=1e10 --dir=../results/polarizer_s_ref --format=NeXus sample=1 omegaa=2.0 over_illumination=0.001  lambda_start=2.0 lambda_end=32.0 sample_length=0.01 enable_polarizer=0"
+command_l_ref="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=2e9 --dir=../results/polarizer_l_ref --format=NeXus sample=1 omegaa=6.0 over_illumination=0.001  lambda_start=2.0 lambda_end=32.0 sample_length=0.048 enable_polarizer=0"
+command_s="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=5e9 --dir=../results/polarizer_s --format=NeXus sample=1 omegaa=2.0 over_illumination=0.001 pol1_start=%%f pol1_angle=%f pol2_start=%%f pol2_angle=%f lambda_start=2.0 lambda_end=32.0 sample_length=0.01 enable_polarizer=3"%(gamma1, gamma2)
+command_l="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=1e9 --dir=../results/polarizer_l --format=NeXus sample=1 omegaa=6.0 over_illumination=0.001 pol1_start=%%f pol1_angle=%f pol2_start=%%f pol2_angle=%f lambda_start=2.0 lambda_end=32.0 sample_length=0.048 enable_polarizer=3"%(gamma1, gamma2)
+
+lengths=linspace(0.1, 0.6, 25)  
+
+def get_polarization(ds):
+    N,ignore=histogram(ds.data['L'], bins=L)
+    PN,ignore=histogram(ds.data['L'], bins=L, weights=sqrt(ds.data['sx']**2+ds.data['sy']**2+ds.data['sz']**2))
+    return PN/maximum(N,1)
+                       
+
+def get_intensity(ds):
+    I, ignore=histogram(ds.data['L'], bins=L, weights=ds.data['p'])
+    return I
+
+def analyze(ds):
+    P=get_polarization(ds)
+    I=get_intensity(ds)
+    T=I/maximum(Iref, Iref[Iref>0].min())
+    return P,T
+
+def get_data(fname):
+    det=McSim(fname)['tof_detector']
+    P,T=analyze(det)
+    return [P[17], P[230], T[17], T[230]]
+
+def simulate_and_analyze(pol1_start, pol2_start):
+    print pol1_start, pol2_start
+    os.system('rm -rf ../results/polarizer_?')
+    os.system(command_s%(pol1_start, pol2_start))
+    os.system(command_l%(pol1_start, pol2_start))
+    
+    global Iref
+    Iref=Iref_s
+    ds=get_data('../results/polarizer_s')
+    Iref=Iref_l
+    dl=get_data('../results/polarizer_l')
+    with open('../results/polarizer_data.dat', 'a') as fh:
+        fh.write((10*'%.6f '+'\n')%tuple([pol1_start, pol2_start]+ds+dl))
+
+if __name__=='__main__':
+    if not (os.path.exists('../results/polarizer_s_ref') and 
+            os.path.exists('../results/polarizer_l_ref')):
+        os.system('rm -rf ../results/polarizer_?_ref')
+        print 'Simulate reference intensities'
+        os.system(command_s_ref)
+        os.system(command_l_ref)
+      
+    
+    ref=McSim('../results/polarizer_s_ref')['tof_detector']
+    Iref_s=get_intensity(ref)
+    ref=McSim('../results/polarizer_l_ref')['tof_detector']
+    Iref_l=get_intensity(ref)
+    
+    with open('../results/polarizer_data.dat', 'w') as fh:
+        fh.write('# pol1   pol2     P3.7_s   P25_s    T3.7_s   T25_s    P3.7_l   P25_l    T3.7_l   T25_l\n')
+    
+    for pol1 in lengths:
+        for pol2 in lengths:
+            simulate_and_analyze(pol1, pol2)
+        with open('../results/polarizer_data.dat', 'a') as fh:
+            fh.write('\n')
+            

simulation/polarizer_coating.sh

@@ -0,0 +1,7 @@
+!#/bin/sh
+rm PolarizerTest.c PolarizerTest.out
+rm -r ../results/polarizer_coating
+mcstas -o PolarizerTest.c PolarizerTest.instr --trace
+#mpicc -O2 -o PolarizerTest.out PolarizerTest.c -lm -DUSE_MPI
+mpicc -O2 -o PolarizerTest.out PolarizerTest.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus -I/usr/local/include/nexus
+LD_LIBRARY_PATH=/usr/local/lib/ mpirun -np 6 ./PolarizerTest.out --ncount=1e8 --format=NeXus --dir=../results/polarizer_coating R0=1.0 gamma=1.66

simulation/run_monitor.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+
+mcstas -o Estia_monitor.c Estia_monitor.instr
+mpicc -O3 -o Estia_monitor.out Estia_monitor.c -lm -DUSE_MPI
+
+
+mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0
+mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_ref --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=2
+mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_trans --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=1
+
+mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor100 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.0001
+mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor001 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.000001

simulation/run_pol.sh

@@ -1,213 +1,77 @@
 #!/bin/bash
 
 DEST=../results
-ncount=2e7
-use_cores=4
+ncount=1e8
+use_cores=16
 sample=4
-omega=0.8
-sample_length=0.01
+sample_length=0.005
 sample_height=0.01
 lambda_start=3.0
 lambda_end=32.0
 
-###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
-if [ Estia_baseline.instr -nt Estia_baseline.out ]; then
-	rm Estia_baseline.c Estia_baseline.out
-	mcstas -o Estia_baseline.c Estia_baseline.instr
-	mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
-fi
-# 
-# if [ Estia_baseline_ana1.instr -nt Estia_baseline_ana1.out ]; then
-# 	rm Estia_baseline_ana1.c Estia_baseline_ana1.out
-# 	mcstas -o Estia_baseline_ana1.c Estia_baseline_ana1.instr
-# 	mpicc -O3 -o Estia_baseline_ana1.out Estia_baseline_ana1.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
-# fi
-# 
-# if [ Estia_baseline_ana2.instr -nt Estia_baseline_ana2.out ]; then
-# 	rm Estia_baseline_ana2.c Estia_baseline_ana2.out
-# 	mcstas -o Estia_baseline_ana2.c Estia_baseline_ana2.instr
-# 	mpicc -O3 -o Estia_baseline_ana2.out Estia_baseline_ana2.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
-# fi
+bash compile_if_needed.sh
 
+###################### Reference measurement 10x5mm² and 10x0.5mm² VS ####################
 
-###################### Reference and Ni-layer measurement 10x10mm² sample ####################
-# ncount=1e10
-sample_length=0.01
+ncount=1e8
+sample_length=0.005
 sample_height=0.01
-lambda_start=3.5
-lambda_end=30.0
-sample=1
 
-# omega=1.0
-# DESTi=$DEST/pol_ref_10x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=0
-# 
-# DESTi=$DEST/pol1_10x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=1 enable_analyzer=0
-# 
-# DESTi=$DEST/pol2_10x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=2 enable_analyzer=0
-# # 
-# omega=7.0
-# DESTi=$DEST/pol1_10x10_70
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=1 enable_analyzer=0
-# 
-# DESTi=$DEST/pol2_10x10_70
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=2 enable_analyzer=0
-# 
-
-# ncount=1e8
-omega=0.0
-# DESTi=$DEST/pol1_10x50_70
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=1 enable_analyzer=0
-
-DESTi=$DEST/pol2_10x50_70
+DESTi=$DEST/pol_ref_10x5
 if [ -e "$DESTi" ]; then
 	rm -r "$DESTi"
 fi
 
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --ncount=$ncount --gravitation \
-        omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
+mpirun  -np $use_cores --hostfile hostnames Estia_baseline.out \
+        --dir="$DESTi" --format=NeXus --ncount=$ncount \
         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=2 enable_analyzer=0
-# 
-# ncount=5e9
-# sample_length=0.01
-# omega=1.0
-# sample_length=0.003
-# DESTi=$DEST/pol_ref_3x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=0
-# 
-# DESTi=$DEST/pol1_3x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=1 enable_analyzer=0
-# 
-# DESTi=$DEST/pol2_3x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=2 enable_analyzer=0
-# 
-# 
-# sample_height=0.01
-# sample_length=0.01
-# DESTi=$DEST/ana1_10x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=1
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
+		enable_polarizer=0 enable_analyzer=0
+       
+ncount=5e8
+sample_length=0.0005
+sample_height=0.01
 
-# sample_height=0.001
-# DESTi=$DEST/ana1_10x1_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline_ana1.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=1
-#        
-# sample_height=0.01
-# DESTi=$DEST/ana2_10x10_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=2
+DESTi=$DEST/pol_ref_10x0p5
+if [ -e "$DESTi" ]; then
+	rm -r "$DESTi"
+fi
 
-# sample_height=0.001
-# DESTi=$DEST/ana2_10x1_10
-# if [ -e "$DESTi" ]; then
-# 	rm -r "$DESTi"
-# fi
-# 
-# mpirun  -np $use_cores Estia_baseline_ana2.out \
-#         --dir="$DESTi" --ncount=$ncount --gravitation \
-#         omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
-#         sample=$sample sample_length=$sample_length sample_height=$sample_height \
-#         lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=0 enable_analyzer=2
+mpirun  -np $use_cores --hostfile hostnames Estia_baseline.out \
+        --dir="$DESTi" --format=NeXus --ncount=$ncount \
+        sample=$sample sample_length=$sample_length sample_height=$sample_height \
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
+		enable_polarizer=0 enable_analyzer=0
        
 
-###################### Reference and Ni-layer measurement 1x1mm² sample ####################
+###################### Polarizers 10x5mm² and 10x0.5mm² VS ####################
 
+ncount=2e8
+sample_length=0.005
+sample_height=0.01
 
+DESTi=$DEST/pol_10x5
+if [ -e "$DESTi" ]; then
+	rm -r "$DESTi"
+fi
+
+mpirun  -np $use_cores --hostfile hostnames Estia_baseline.out \
+        --dir="$DESTi" --format=NeXus --ncount=$ncount \
+        sample=$sample sample_length=$sample_length sample_height=$sample_height \
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
+		enable_polarizer=1 enable_analyzer=0
+       
+ncount=1e9
+sample_length=0.0005
+sample_height=0.01
+
+DESTi=$DEST/pol_10x0p5
+if [ -e "$DESTi" ]; then
+	rm -r "$DESTi"
+fi
+
+mpirun  -np $use_cores --hostfile hostnames Estia_baseline.out \
+        --dir="$DESTi" --format=NeXus --ncount=$ncount \
+        sample=$sample sample_length=$sample_length sample_height=$sample_height \
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
+		enable_polarizer=1 enable_analyzer=0

simulation/run_selene.py

@@ -2,45 +2,51 @@
 #-*- coding: utf8 -*-
 
 import sys, os
-import numpy as np
 from subprocess import call
+from numpy import *
+from scipy.optimize import leastsq
 
-FOLDER='../results2/'
-PREFIX='selene_geo_'
-SEED=3798
-ITEMS=1000
 
-CALL='mpirun -np 6 Selene_geometry.out --dir=%s --ncount=4e8 enable_gravity=0 sample_length=0.003 '
-CALL+='tx_1=%f tz_1=%f ry_1=%f rz_1=%f tx_2=%f tz_2=%f ry_2=%f rz_2=%f '
+CALL='/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out '
+CALL=+'--dir=../results/selene_geo --ncount=1e9 '
+CALL=+'omegaa=1.0 sample=4 sample_length=0.001 sample_height=0.01 '
+CALL=+'lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 '
+CALL=+'lambda_min=3.75 selene1_foot1y=%.4f selene1_foot2y=%.4f '
 
-TX=[-0.001, 0.001]
-TZ=[-0.005, 0.005]
-RY=[-0.001, 0.001]
-RZ=[-2.5, 2.5]
+def B(x,w):
+    # box function with full width w
+    return float32(abs(x)<=(w/2.))
 
+def G(x,I0,x0,sigma):
+    # Gaussian with intensity I0, center x0 and standard deviation sigma
+    return I0*exp(-0.5*(x-x0)**2/sigma**2)
+
+def Intensity(x, p):
+    I0, x0, sigma, w=p
+    return convolve(B(xc,w), G(xc, I0, x0, sigma), mode='same')
+
+def Beam(x,I0,x0,sigma,w):
+    I=I0*where((x-x0)<(-w/2.), exp(-0.5*(x-x0+w/2.)**2/sigma**2),
+            where((x-x0)>(w/2.), exp(-0.5*(x-x0-w/2.)**2/sigma**2), 1.))
+    return I
+
+def residuals(p, x, y):
+    I0, x0, sigma, w=p
+    return y-Beam(x, I0, x0, sigma, w)
+
+def FWHM(pi):
+    rng=xc[where(Beam(xc, *pi)>=(pi[0]/2.))[0]]
+    return rng[-1]-rng[0]
+  
+x=linspace(-10.0005, 10.0005, 20001)
+xc=(x[1:]+x[:-1])/2.
+
+def analyze(fi):
+  sim=mr.McSim(fi)
+  det=sim['tof_sample']
+  ignore,y=det.project1d('x', bins=x/1000.)
+  p,res=leastsq(residuals, (y.max(), (xc*y).sum()/y.sum(), 0.01, 0.1), (xc, y))
+  return p,det
 
 if __name__=='__main__':
-    call('mcstas -t -o Selene_geometry.c Selene_geometry.instr'.split())
-    call('mpicc -O2 -o Selene_geometry.out Selene_geometry.c -lm -DUSE_MPI'.split())
-    
-    np.random.seed(SEED)
-    fh=open(os.path.join(FOLDER, PREFIX+'sims.dat'), 'w')
-    fh.write('# index tx_1 tz_1 ry_1 rz_1 tx_2 tz_2 ry_2 rz_2\n')
-    fh.write('# seed = %i\n'%SEED)
-    
-    for i in range(ITEMS):
-        tx_1=np.random.random()*(TX[1]-TX[0])+TX[0]
-        tx_2=np.random.random()*(TX[1]-TX[0])+TX[0]
-        tz_1=np.random.random()*(TZ[1]-TZ[0])+TZ[0]
-        tz_2=np.random.random()*(TZ[1]-TZ[0])+TZ[0]
-        ry_1=np.random.random()*(RY[1]-RY[0])+RY[0]
-        ry_2=np.random.random()*(RY[1]-RY[0])+RY[0]
-        rz_1=np.random.random()*(RZ[1]-RZ[0])+RZ[0]
-        rz_2=np.random.random()*(RZ[1]-RZ[0])+RZ[0]
-        ln=' '.join(['%f'%val for val in [tx_1, tz_1, ry_1, rz_1, tx_2, tz_2, ry_2, rz_2]])
-        fh.write(ln+'\n')
-        print i, ':', ln
-        call((CALL%(os.path.join(FOLDER, PREFIX+'%05i/'%i), tx_1, tz_1, ry_1, rz_1, tx_2, tz_2, ry_2, rz_2)).split())
-    
-    fh.close()
-    
+    pass

simulation/sbatch.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH -J mcEstia_1
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=1-00:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout.log
+#SBATCH -e stderr.log
+
+#SBATCH --partition=ll_long
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+/usr/bin/modulecmd tcsh load mcstas
+bash compile_if_needed.sh
+
+for y1 in $(seq -0.009 0.001 0.009)
+  do
+  for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
+    do
+    printf -v Y1 %.3f $y1
+    printf -v Y2 %.3f $y2
+    echo $Y1 $Y2
+    /afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
+            --dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
+            omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
+            lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
+            lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
+  done
+done
+
+        
+echo "Program finished with exit code $? at: `date`"
+

simulation/sbatch.tsh

@@ -0,0 +1,38 @@
+#!/bin/tcsh
+#SBATCH -J mcEstia_S
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=2:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout.log
+#SBATCH -e stderr.log
+
+#SBATCH --partition=ll_short
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+#module unload intel
+#module load intel
+module load mcstas
+bash compile_if_needed.sh
+
+mpirun -np $SLURM_NPROCS Estia_baseline.out --format=NeXus \
+        --dir=$*
+#        --dir=../results/$1 --ncount=1e9 \
+#        sample=1 enable_polarizer=1 \
+#        omegaa=6.0 sample_length=0.048 sample_height=0.01
+
+echo "Program finished with exit code $? at: `date`"
+
+
+echo "Compressing hdf5 file..."
+h5repack -i ../results/$1/mccode.h5 -o ../results/$1/mccode.h5.c -f /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:GZIP=5 -l /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:CHUNK=3072x8
+mv ../results/$1/mccode.h5.c ../results/$1/mccode.h5
+
+echo "Scipt finished."

simulation/sbatch_2.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH -J mcEstia_2
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=1-00:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout2.log
+#SBATCH -e stderr2.log
+
+#SBATCH --partition=ll_long
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+/usr/bin/modulecmd tcsh load mcstas
+bash compile_if_needed.sh
+
+for y1 in $(seq -0.09 0.01 -0.01)
+  do
+  for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
+    do
+    printf -v Y1 %.3f $y1
+    printf -v Y2 %.3f $y2
+    echo $Y1 $Y2
+    /afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
+            --dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
+            omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
+            lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
+            lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
+  done
+done
+
+        
+echo "Program finished with exit code $? at: `date`"
+

simulation/sbatch_3.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH -J mcEstia_3
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=1-00:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout3.log
+#SBATCH -e stderr3.log
+
+#SBATCH --partition=ll_long
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+/usr/bin/modulecmd tcsh load mcstas
+bash compile_if_needed.sh
+
+for y1 in $(seq 0.01 0.01 0.09)
+  do
+  for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
+    do
+    printf -v Y1 %.3f $y1
+    printf -v Y2 %.3f $y2
+    echo $Y1 $Y2
+    /afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
+            --dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
+            omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
+            lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
+            lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
+  done
+done
+
+        
+echo "Program finished with exit code $? at: `date`"
+

simulation/sbatch_4.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH -J mcEstia_4
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=1-00:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout4.log
+#SBATCH -e stderr4.log
+
+#SBATCH --partition=ll_long
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+/usr/bin/modulecmd tcsh load mcstas
+bash compile_if_needed.sh
+
+for y1 in $(seq -0.5 0.1 -0.1)
+  do
+  for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
+    do
+    printf -v Y1 %.3f $y1
+    printf -v Y2 %.3f $y2
+    echo $Y1 $Y2
+    /afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
+            --dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
+            omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
+            lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
+            lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
+  done
+done
+
+        
+echo "Program finished with exit code $? at: `date`"
+

simulation/sbatch_5.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH -J mcEstia_5
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --time=1-00:00:00
+#SBATCH --mail-type=fail
+#SBATCH --mail-user=artur.glavic@psi.ch
+
+#SBATCH -o stdout5.log
+#SBATCH -e stderr5.log
+
+#SBATCH --partition=ll_long
+
+echo "Starting at `date`"
+echo "Running on hosts: $SLURM_NODELIST"
+echo "Running on $SLURM_NNODES nodes."
+echo "Running on $SLURM_NPROCS processors."
+echo "Current working directory is `pwd`"
+
+/usr/bin/modulecmd tcsh load mcstas
+bash compile_if_needed.sh
+
+for y1 in $(seq 0.1 0.1 0.5)
+  do
+  for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
+    do
+    printf -v Y1 %.3f $y1
+    printf -v Y2 %.3f $y2
+    echo $Y1 $Y2
+    /afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
+            --dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
+            omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
+            lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
+            lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
+  done
+done
+
+        
+echo "Program finished with exit code $? at: `date`"
+