Add pulse skipping simulation to run script

.gitignore

@@ -8,7 +8,6 @@ old
 *.backup
 .project
 .settings
-.idea
 *.pyc
 *.pyo
 hosts

analysis/compress_h5.sh

@@ -4,6 +4,6 @@
 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
+  h5repack -i $fi -o $fi.c -f /entry1/data/of_detector_list_p_x_y_t_L_sx_sy/events:GZIP=5 -l /entry1/data/of_detector_list_p_x_y_t_L_sx_sy/events:CHUNK=3072x7
   mv $fi.c $fi
 done

analysis/estia_analyze.py

@@ -8,149 +8,38 @@ import sys, os
 from numpy import *
 seterr(all='ignore')
 import estia_help as eh
-from aglib import mcstas_reader as mr
-
-
-# Scaling factor of source monitor after normalizing by area to get to brilliance.
-# Converts from 1.5x4.0 degree² acceptance range to 1/steradian
-# Factor 10 due to 0.1 Angstrom binning.
-BT_SOURCE_SCALE=10.*(180.**2/pi**2)/1.5/4.0
-BT_SAMPLE_SCALE=10.*(180.**2/pi**2)/1.5/1.5
+import mcstas_reader as mr
 
+# work around for wrong wavelength calculation in pulse skipping mode
+eh.l_code='((t-0.0335)%%(%f/14.)+0.0335-0.0015)*101.436605'
 
 if __name__=='__main__':
-  if not os.path.exists('../results/analyzed'):
-    os.mkdir('../results/analyzed')
+  if not os.path.exists('../results/mooc_analyzed'):
+    os.mkdir('../results/mooc_analyzed')
   
-  if len(sys.argv)>1:
-    only_items=map(int, sys.argv[1:])
-  else:
-    only_items=range(100)
 
-  if 1 in only_items:
-    print "Brilliance Transfer 5x10"
-    bt_10=mr.McSim('../results/brilliance_5x10/mccode.h5')
-    d=bt_10['tof_sample.L']
-    vs=bt_10['tof_virtual_source.L']
-    r=bt_10['tof_source.L']
-    x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*2., r, BT_SOURCE_SCALE)
-    savetxt('../results/analyzed/brilliance_transfer_5x10.dat', array([x, y, e]).T)
-    savetxt('../results/analyzed/brilliance_5x10.dat', array([x, d.data*BT_SAMPLE_SCALE*2.,
-                                           d.errors*BT_SAMPLE_SCALE*2.]).T)
-    x, y, e=eh.calc_brilliance_transfer(vs, BT_SAMPLE_SCALE, r, BT_SOURCE_SCALE)
-    savetxt('../results/analyzed/brilliance_transfer_vs_5x10.dat', array([x, y, e]).T)
-    
-    btn_10=mr.McSim('../results/brilliance_nowindow_5x10/mccode.h5')
-    d=btn_10['tof_sample.L']
-    vs=btn_10['tof_virtual_source.L']
-    r=btn_10['tof_source.L']
-    x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*2., r, BT_SOURCE_SCALE)
-    savetxt('../results/analyzed/brilliance_transfer_nowindow_5x10.dat', array([x, y, e]).T)
-    x, y, e=eh.calc_brilliance_transfer(vs, BT_SAMPLE_SCALE, r, BT_SOURCE_SCALE)
-    savetxt('../results/analyzed/brilliance_transfer_nowindow_vs_5x10.dat', array([x, y, e]).T)
+  print "pulse skipping mode"
+  opts=dict(qres=0.02, qmin=0.007, qmax=0.2, mindq=2e-5, crop_overlap=False)
+  names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)', 
+         '1s-Stat(ToF)', '1pulse-Stat(ToF)', '10s-Statistics(ToF)', 
+         'Intensity(λ)', 'Reflectivity(λ)']
+  units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
+  info='Reflectivity simulation for 50x10mm² sample, 2-pulse skipping'
 
-
-  if 2 in only_items:
-    print "Brilliance Transfer 1x1"
-    bt_1=mr.McSim('../results/brilliance_1x1/mccode.h5')
-    d=bt_1['tof_sample.L']
-    r=bt_1['tof_source.L']
-    x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*100, r, BT_SOURCE_SCALE)
-    savetxt('../results/analyzed/brilliance_transfer_1x1.dat', array([x, y, e]).T)
-    savetxt('../results/analyzed/brilliance_1x1.dat', array([x, d.data*BT_SAMPLE_SCALE*100,
-                                           d.errors*BT_SAMPLE_SCALE*100]).T)
-
-  if 3 in only_items:
-    print "Refelctiviy 10x10 sample"
-    opts=dict(qres=0.01, qmin=0.003, qmax=0.5, mindq=5e-4)
-    names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)', 'dR (1s ToF)', 'Intensity(λ)', 'Reflectivity(λ)']
-    units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
-    info='Reflectivity simulation for 10x10mm² sample, omega=%.1f deg'
-
-    print "0.8 degree"
-    ref=mr.McSim('../results/reference_10x10_10/mccode.h5')['tof_detector']
-    ni_08=mr.McSim('../results/nickle_10x10_08/mccode.h5')['tof_detector']
-    q, I, R=eh.calcR(ni_08, ref, 0.8, 1.0, detcorr=True, **opts)
-    q, I_real, R_real=eh.calcR(ni_08, ref, 0.8, 1.0, use_tof=False, **opts)
-    idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
-    eh.save_w_header('../results/analyzed/reflectivity_10x10_08.dat',
-                     [q[idx_start:idx_end+1],
-                      I[idx_start:idx_end+1], R[idx_start:idx_end+1],
-                      (R/sqrt(I))[idx_start:idx_end+1],
-                      I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
-                     info, names, units)
-    del(ni_08)
-
-    print "3.0 degree"
-    ni_30=mr.McSim('../results/nickle_10x10_30/mccode.h5')['tof_detector']
-    q, I, R=eh.calcR(ni_30, ref, 3.0, 1.0, detcorr=True, **opts)
-    q, I_real, R_real=eh.calcR(ni_30, ref, 3.0, 1.0, use_tof=False, **opts)
-    idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
-    eh.save_w_header('../results/analyzed/reflectivity_10x10_30.dat',
-                     [q[idx_start:idx_end+1],
-                      I[idx_start:idx_end+1], R[idx_start:idx_end+1],
-                      (R/sqrt(I))[idx_start:idx_end+1],
-                      I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
-                     info, names, units)
-    del(ni_30)
-
-    print "8.0 degree"
-    ni_80=mr.McSim('../results/nickle_10x10_80/mccode.h5')['tof_detector']
-    q, I, R=eh.calcR(ni_80, ref, 8.0, 1.0, detcorr=True, **opts)
-    q, I_real, R_real=eh.calcR(ni_80, ref, 8.0, 1.0, use_tof=False, **opts)
-    idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
-    eh.save_w_header('../results/analyzed/reflectivity_10x10_80.dat',
-                     [q[idx_start:idx_end+1],
-                      I[idx_start:idx_end+1], R[idx_start:idx_end+1],
-                      (R/sqrt(I))[idx_start:idx_end+1],
-                      I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
-                     info, names, units)
-
-  if 4 in only_items:
-    print "10x10 sample pulse skipping mode."
-    opts=dict(qres=0.04, qmin=0.003, qmax=0.2, mindq=2e-4)
-    names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)', 'Intensity(λ)', 'Reflectivity(λ)']
-    units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
-    info='Reflectivity simulation for 10x10mm² sample, 2-pulse skipping'
-
-    ref=mr.McSim('../results/single_skip_reference/mccode.h5')['tof_detector']
-    ni=mr.McSim('../results/single_skip_nickle/mccode.h5')['tof_detector']
-    q, I, R=eh.calcR(ni, ref, 2.0, skip_pulses=2, **opts)
-    q, I_real, R_real=eh.calcR(ni, ref, 2.0, skip_pulses=2, use_tof=False, **opts)
+  ref=mr.McSim('../results/mooc_reference_ps')['tof_detector']
+  
+  for i in range(6):
+    print "dataset",i
+    sample=mr.McSim('../results/mooc_model_ps_%i'%i)['tof_detector']
+    q, I, R=eh.calcR(sample, ref, 1.5, 1.5, skip_pulses=2, **opts)
+    q, I_real, R_real=eh.calcR(sample, ref, 1.5, 1.5, skip_pulses=2, use_tof=False, **opts)
     idx_start, idx_end=where((I>0.))[0][[0,-1]]
-    eh.save_w_header('../results/analyzed/single_skip_reflectivity.dat',
-                     [q[idx_start:idx_end+1],
-                      I[idx_start:idx_end+1], R[idx_start:idx_end+1],
+    eh.save_w_header('../results/mooc_analyzed/model%i_skip_reflectivity.dat'%i,
+                      [q[idx_start:idx_end+1],
+                      I[idx_start:idx_end+1], R[idx_start:idx_end+1], 
+                      random.poisson(I[idx_start:idx_end+1])/I[idx_start:idx_end+1]*R[idx_start:idx_end+1],
+                      random.poisson(I[idx_start:idx_end+1]*3./14.)/I[idx_start:idx_end+1]*14./3.*R[idx_start:idx_end+1],
+                      random.poisson(I[idx_start:idx_end+1]*10.)/I[idx_start:idx_end+1]/10.*R[idx_start:idx_end+1],
                       I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
-                     info, names, units)
+                      info, names, units)
 
-  if 5 in only_items:
-    print "Extracting event statistics."
-    fh=open('../results/analyzed/event_stats.dat', 'w')
-    fh.write('# event statistics for Estia\n')
-    fh.write('# item      avg. total     peak total     avg. ROI       peak ROI\n')
-    fh.write('#           [cps]          [cps]          [cps/mm²]      [cps/mm²]\n')
-    line='%(name)s  %(total) 12e  %(peak) 12e  %(roi) 12e  %(roi_peak) 12e\n'
-
-    ds=mr.McSim('../results/reference_10x10_10/mccode.h5')['tof_detector']
-    res=eh.calcStat(ds)
-    res['name']='Reference'
-    fh.write(line%res)
-
-    ds=mr.McSim('../results/nickle_10x10_08/mccode.h5')['tof_detector']
-    res=eh.calcStat(ds)
-    res['name']='Ni-0.8deg'
-    fh.write(line%res)
-
-    ds=mr.McSim('../results/nickle_10x10_30/mccode.h5')['tof_detector']
-    res=eh.calcStat(ds)
-    res['name']='Ni-3.0deg'
-    fh.write(line%res)
-
-    ds=mr.McSim('../results/nickle_10x10_80/mccode.h5')['tof_detector']
-    res=eh.calcStat(ds)
-    res['name']='Ni-8.0deg'
-    fh.write(line%res)
-
-    fh.close()
-    

analysis/mcstas_reader.py

@@ -0,0 +1,422 @@
+#-*- coding: utf-8 -*-
+'''
+Simple support library for reading, analyzing and plotting of McStas results
+from the Estia instrument simulations.
+Meant to be used from IPython Notebook or QtConsole, but can also be run stand alone.
+'''
+
+import os, sys
+from numpy import *
+
+try:
+    import h5py
+except ImportError:
+    print "h5py not found, modern NeXuS format will not be readable."
+
+try:
+  from IPython import display #@UnusedImport
+  from IPython.core.pylabtools import print_figure
+  from matplotlib.figure import Figure
+  from matplotlib.backends.backend_agg import FigureCanvasAgg
+  from matplotlib.colors import LogNorm
+except ImportError:
+  # IPython and/or matplotlib not available, no interactive functionalities
+  display=None
+
+MAX_EVTS_BATCH=50000
+    
+class McSim(object):
+  '''
+  Object representing complete simulation from McStas.
+  Supports either the old style McStas output with separate ASCII files 
+  or HDF5 single file output.
+  
+  Different monitors can be accessed as keys like in a dictionary. The data is only loaded when the monitor is first accessed.
+  '''
+
+  def __init__(self, path):
+    '''
+    Initialize the simulation. path should be the file name to either the NeXuS or the mccode.sim file.
+    '''
+    self._data={}
+    if path.endswith('.h5'):
+      self._init_hdf(path)
+    elif path.endswith('.sim'):
+      self._init_old(path)
+    else:
+      if os.path.exists(os.path.join(path,'mccode.h5')):
+        self._init_hdf(os.path.join(path,'mccode.h5'))
+      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
+
+  def _init_hdf(self, path):
+    self.hdf=h5py.File(path, 'r')
+    self.data_loader=DataLoaderHDF(self.hdf)
+    self.info=self.data_loader.info
+
+  def _init_old(self, path):
+    self.info=HeaderFile(path)
+    self.data_loader=DataLoaderOld(self.info, os.path.dirname(path))
+
+  def keys(self):
+    return self.info['data'].keys()
+
+  def monitors(self):
+    output=[]
+    for val in self.info['data'].values():
+      xy=(val['xvar'], val['yvar'])
+      if not xy in output:
+        output.append(xy)
+    output.sort()
+    return output
+  
+  def plot(self, monitors=None):
+    graphs={}
+    for key, val in self.info['data'].items():
+      xy=(val['xvar'], val['yvar'])
+      if monitors is not None and xy!=monitors:
+        continue
+      data=self[key]
+      if xy in graphs:
+        graphs[xy].append(data)
+      else:
+        graphs[xy]=[data]
+    for xy, datasets in sorted(graphs.items()):
+      cols=min(len(datasets), 3)
+      rows=len(datasets)/3+1
+
+      fig=Figure(figsize=(12, 5*rows), dpi=300, facecolor='#FFFFFF')
+      FigureCanvasAgg(fig)
+
+      for i, data in enumerate(datasets):
+        ax=fig.add_subplot(rows, cols, i+1)
+        data.plot(ax=ax)
+      graphs[xy]=fig
+    if monitors is None:
+      return graphs
+    else:
+      return fig
+
+  def __getitem__(self, item):
+    if item in self._data:
+      return self._data[item]
+    elif item in self.keys():
+      data=self.data_loader.load_item(item)
+      self._data[item]=data
+      return data
+    else:
+      raise KeyError, "Can't find dataset %s"%item
+
+class HeaderFile(object):
+  '''
+  Analyze McStas mccode.sim header files.
+  '''
+  _data=None
+
+  @property
+  def data(self):
+    return self._data['data']
+
+  def __init__(self, path):
+    self._data={}
+    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.'
+    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)
+
+  def get_data(self, data):
+    output={}
+    start_idx=0
+    while start_idx<len(data):
+      start_idx=data.find('begin data', start_idx)
+      if start_idx==-1:
+        break
+      end_idx=data.find('end data', start_idx)
+      block=data[start_idx:end_idx].splitlines()[1:]
+      block_info={}
+      for bline in block:
+        item, value=bline.strip().split(':', 1)
+        block_info[item]=value.strip()
+        if item=='component':
+          output[value.strip()]=block_info
+
+      start_idx=end_idx
+    return output
+
+  def __getitem__(self, item):
+    return self._data[item]
+
+class DataLoaderOld(object):
+  '''
+    Load and analyze a old style McStas format with a simulation and a set of data text files.
+  '''
+  def __init__(self, info, root):
+    self.info=info
+    self.root=root
+
+  def load_item(self, item):
+    item_info=self.info['data'][item]
+    if not item_info['type'].startswith('array_2d'):
+      return self.load_item_1d(item)
+
+    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     
+      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]
+      return Dataset(data, item_info)
+
+  def load_item_1d(self, item):
+    item_info=self.info['data'][item]
+    fname=os.path.join(self.root, item_info['filename'])
+    raw=loadtxt(fname).T
+    data=raw[1]
+    errors=raw[2]
+    return Dataset1D(data, errors, item_info)
+
+class DataLoaderHDF(object):
+  '''
+    Load and analyze a new style NeXuS file format.
+  '''
+  def __init__(self, hdf):
+    self.hdf=hdf['entry1']
+    self.info={}
+    self.info['data']={}
+    for item in self.hdf['data'].keys():
+      node=self.hdf['data/'+item]
+      info={}
+      for key, value in node.attrs.items():
+        info[key.strip()]=value.strip()
+      if info['filename'][-4:]=='.dat' or '_list.' in info['filename']:
+        self.info['data'][info['component']]=info
+      else:
+        self.info['data'][info['filename']]=info
+      info['datapath']='data/'+item
+
+  def load_item(self, item):
+    item_info=self.info['data'][item]
+    if not item_info['type'].startswith('array_2d'):
+      return self.load_item_1d(item)
+
+    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
+      cols=item_info['variables'].split()
+      evds=node['events']
+      if len(evds)<=MAX_EVTS_BATCH:
+        data=evds.value.astype(float32).view(
+            dtype={'names': cols, 'formats': ['f4']*len(cols)}).flatten()
+      else:
+        ds=[]
+        sys.stdout.write('Reading large dataset:\n')
+        for i in range(len(evds)//MAX_EVTS_BATCH+1):
+          sys.stdout.write('\r%i/%i'%(i*MAX_EVTS_BATCH, len(evds)))
+          sys.stdout.flush()
+          ds.append(evds[i*MAX_EVTS_BATCH:(i+1)*MAX_EVTS_BATCH])
+        sys.stdout.write('\r%i/%i\n'%(len(evds), len(evds)))
+        data=vstack(ds).astype(float32).view(
+            dtype={'names': cols, 'formats': ['f4']*len(cols)}).flatten()
+      return TofData(data, item_info)
+    else:
+      data=node['data'].value.T
+      return Dataset(data, item_info)
+  
+  def load_item_1d(self, item):
+    item_info=self.info['data'][item]
+    node=self.hdf[item_info['datapath']]
+    data=node['data'].value
+    errors=node['errors'].value
+    return Dataset1D(data, errors, item_info)
+
+class Dataset1D(object):
+  '''
+    Representation of a standard McStas dataset of one variable.
+  '''
+  def __init__(self, data, errors, info):
+    self.data=data
+    self.errors=errors
+    self.info=info
+
+  def plot(self, log=False, ax=None):
+    if ax is None:
+      import pylab
+      ax=pylab.gca()
+
+    limits=map(float, self.info['xlimits'].split())
+    x=linspace(limits[0], limits[1], len(self.data))
+
+    ax.errorbar(x, self.data, yerr=self.errors)
+    if log:
+      ax.set_yscale('log')
+    else:
+      ax.set_yscale('linear')
+    ax.set_xlabel(self.info['xlabel'])
+    ax.set_ylabel(self.info['ylabel'])
+    ax.set_title(self.info['component'])
+
+  def _repr_png_(self):
+    '''
+      Image representation form ipython console/notebook. 
+    '''
+    fig=Figure(figsize=(8, 5), dpi=300, facecolor='#FFFFFF')
+    FigureCanvasAgg(fig)
+    ax=fig.add_subplot(111)
+
+    self.plot(ax=ax)
+
+    return print_figure(fig, dpi=72)
+
+
+class Dataset(object):
+  '''
+    Representation of a standard McStas dataset with 2D view.
+  '''
+  def __init__(self, data, info):
+    self.data=data
+    self.info=info
+
+  def plot(self, log=False, ax=None):
+    if ax is None:
+      import pylab
+      ax=pylab.gca()
+
+    limits=map(float, self.info['xylimits'].split())
+
+    if log:
+      ax.imshow(self.data, extent=limits, aspect='auto', norm=LogNorm())
+    else:
+      ax.imshow(self.data, extent=limits, aspect='auto')
+    ax.set_xlabel(self.info['xlabel'])
+    ax.set_ylabel(self.info['ylabel'])
+    ax.set_title(self.info['component'])
+
+  def _repr_png_(self):
+    '''
+      Image representation form ipython console/notebook. 
+    '''
+    fig=Figure(figsize=(8, 5), dpi=300, facecolor='#FFFFFF')
+    FigureCanvasAgg(fig)
+    ax=fig.add_subplot(111)
+
+    self.plot(ax=ax)
+
+    return print_figure(fig, dpi=72)
+
+class TofData(Dataset):
+  '''
+    Representation of a dataset collected with Monitor_nD
+  '''
+  
+  def project1d(self, col, bins=50, fltr=None, newcols=None, norm=None):
+    ''' 
+      Generate binned data for arbitrary binning and columns.
+      
+      The fltr argument can be a string with a filtering condition on the available columns.
+      Syntax for this filtering follows numpy array convetions like (x>5)&(abs(L)<0.3) would
+      be a valid statement.
+      
+      newcols can be a list of (name, code) tuples, that calculate new columns from
+      existing ones, like newcols=[('q', '4*pi/L*sin(theta)')].
+    '''
+    columns=dict([(coli, self.data[coli]) for coli in self.data.dtype.names])
+    if newcols is not None:
+      for name, code in newcols:
+        columns[name]=eval(code, globals(), columns)
+
+    if norm is None:
+      w=self.data['p']
+    else:
+      w=eval(norm+'*p', globals(), columns)
+
+    if fltr is None:
+      I, x=histogram(columns[col], bins=bins, weights=w)
+    else:
+      if isinstance(fltr, basestring):
+        fltr=eval(fltr, globals(), columns)
+      I, x=histogram(columns[col][fltr], bins=bins, weights=w[fltr])
+    return x, I
+
+  def project2d(self, xcol, ycol, bins=50, fltr=None, newcols=None):
+    ''' 
+      Generate 2D binned data for arbitrary binning and columns.
+      
+      The fltr argument can be a string with a filtering condition on the available columns.
+      Syntax for this filtering follows numpy array convetions like (x>5)&(abs(L)<0.3) would
+      be a valid statement.
+      
+      newcols can be a list of (name, code) tuples, that calculate new columns from
+      existing ones, like newcols=[('q', '4*pi/L*sin(theta)')].
+    '''
+    columns=dict([(coli, self.data[coli]) for coli in self.data.dtype.names])
+    if newcols is not None:
+      for name, code in newcols:
+        columns[name]=eval(code, globals(), columns)
+
+    if fltr is None:
+      I, y, x=histogram2d(columns[ycol], columns[xcol],
+                        bins=bins, weights=self.data['p'])
+    else:
+      if isinstance(fltr, basestring):
+        fltr=eval(fltr, globals(), columns)
+      I, y, x=histogram2d(columns[ycol][fltr], columns[xcol][fltr],
+                        bins=bins, weights=columns['p'][fltr])
+    return x, y, I
+
+  def plot(self, xcol='x', ycol='y', log=False, ax=None, bins=50, fltr=None, newcols=None):
+    if ax is None:
+      import pylab
+      ax=pylab.gca()
+
+    x, y, I=self.project2d(xcol, ycol, bins=bins, fltr=fltr, newcols=newcols)
+
+    if log:
+      ax.pcolormesh(x, y, I, norm=LogNorm())
+    else:
+      ax.pcolormesh(x, y, I)
+
+    col_names=self.info['title'].split()
+    cols=self.info['variables'].split()
+
+    try:
+      ax.set_xlabel(xcol+'-'+col_names[cols.index(xcol)])
+    except ValueError:
+      ax.set_xlabel(xcol)
+    try:
+      ax.set_ylabel(ycol+'-'+col_names[cols.index(ycol)])
+    except ValueError:
+      ax.set_ylabel(ycol)
+    ax.set_title(self.info['component'])
+
+
+  def plot1d(self, col='x', log=False, ax=None, bins=50, fltr=None, newcols=None):
+    if ax is None:
+      import pylab
+      ax=pylab.gca()
+
+    x, I=self.project1d(col, bins=bins, fltr=fltr, newcols=newcols)
+
+    if log:
+      ax.semilogy((x[:-1]+x[1:])/2., I)
+    else:
+      ax.plot((x[:-1]+x[1:])/2., I)
+    col_names=self.info['title'].split()
+    cols=self.info['variables'].split()
+
+    try:
+      ax.set_xlabel(col+'-'+col_names[cols.index(col)])
+    except ValueError:
+      ax.set_xlabel(col)
+    ax.set_ylabel('Intensity')
+    ax.set_title(self.info['component'])
+

simulation/Estia_baseline.instr

@@ -49,14 +49,13 @@
 *******************************************************************************/
 
 DEFINE INSTRUMENT ESS_reflectometer_Estia
-(double omegaa = 1.2,  int sample = 0,  double sample_length   = 0.01, double sample_height   = 0.01,
+(double omegaa = 1.2,  int sample = 1,  string sample_file="alessandra_model_2_nores_000.dat", 
+double sample_length   = 0.01, double sample_height   = 0.01,
 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 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
+ int enable_polarizer = 0, int enable_analyzer = 0,
+double source_power = 2, double frame_usage     =   0.98
 )
 
 DECLARE
@@ -72,12 +71,8 @@ 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 1 (Will follow later with detailed parameters)
+
 
 // Selene 2
 
@@ -95,7 +90,7 @@ double chopper_phase          ;       // deg phase between pulse and chopper ope
 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 chopper_freq = 14.0;     // Hz  chopper frequency
 
 double slit_distance   =   1.775; //m, distance slit to sample
 
@@ -103,35 +98,32 @@ 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 analyzer_flipper_start  = 0.4; // Distance from sample to analyzer flipper
-
-
 double analyzer1_start  = 0.65; // Distance to sample to start the first analyzer
 double analyzer1_length = 1.3; // length of first analyzer
 double analyzer2_start  = 0.7; // Distance to sample to start the first analyzer
 double analyzer2_length = 1.6; // length of first analyzer
-
 double Theta1_analyzer1, Theta1_analyzer2, Theta2_analyzer1, Theta2_analyzer2, dist_ana_vfocus; // quantities calculated out of values above and 1.5 degree covered divergence
-int p_int;
-%}
 
-USERVARS %{
-  int p_int; // a flag that gets incremented if a polarizer mirror scatters
 %}
 
 INITIALIZE
 %{
+
+chopper_freq/=fmax(1.0, enable_chopper);
 // 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 = 98.0;
-p_int=0;
+chopper_open =360.0*(chopper_pos/(total_length+detector_arm))*frame_usage;
+
 
 /* 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
@@ -141,11 +133,7 @@ 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
@@ -172,12 +160,10 @@ 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, -selene1_center) RELATIVE arm_selene1_center
+    AT (0, 0, 2*NBOA_c) RELATIVE ISCS
+    ROTATED (0, 0, 0) RELATIVE ISCS
 COMPONENT arm_selene2 = Arm()
     AT (0, 0, 2*NBOA_c+2*selene_c) RELATIVE ISCS
     ROTATED (0, 0, 0) RELATIVE ISCS
@@ -211,11 +197,10 @@ 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=3, Lmin = lambda_start, Lmax = lambda_end,
+    target_index=4, Lmin = lambda_start, Lmax = lambda_end,
     n_pulses = 1+enable_chopper, acc_power=source_power)
     AT (0, 0, 0) RELATIVE origin 
 
-
 /***************************************
  * Geometry of neutron feeder separate *
  ***************************************/
@@ -224,14 +209,11 @@ COMPONENT moderator = ESS_butterfly(
 /**************************************************** 
  * 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
+/* 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
 
-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)
+COMPONENT FeNi_out = Slit(xwidth=0.013, yheight=0.038)
     AT (0, 0, -0.220) RELATIVE arm_virtual_source_beam
 
 
@@ -239,46 +221,34 @@ COMPONENT CPC1_out = Slit(xwidth=0.013, yheight=0.038)
 COMPONENT chopper = DiskChopper(radius=chopper_diameter/2.0, yheight=0.02,
     theta_0=chopper_open, phase=chopper_phase+chopper_open/2.0,
     nu=chopper_freq,  nslit=1)
-    WHEN enable_chopper==1
+    WHEN enable_chopper!=0
     AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
     
 
-COMPONENT vs_divergence_h = DivPos_monitor(nb=21, ndiv=41, filename="vs_hordiv",
-          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0,
-		  restore_neutron=1)
-		  AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
-
-COMPONENT vs_divergence_v = DivPos_monitor(nb=21, ndiv=41, filename="vs_verdiv",
-          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0,
-		  restore_neutron=1, vertical=1)
-		  AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
-
 /* The actual virtual source mask, two L-shaped absorbers (first top-right) */
 COMPONENT virtual_source_TR = Slit(
     xmin = 0.0, xmax = 1.0, ymin = -1.0, ymax = sample_height/2+over_illumination*5)
     WHEN sample!=4
     AT (-over_illumination, 0, -0.5*sample_length) RELATIVE arm_virtual_source
 
-// window to cut down to defined size for test setting
-COMPONENT virtual_source_HC = Slit(xwidth=sample_length, yheight=sample_height)
+// 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)
     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
 
-
 /*************************************
  * Geometry of Selene guide separate *
  *************************************/
-%include "Estia_selene1.instr"
+%include "Estia_selene.instr"
 
-%include "Estia_mf.instr"
-
-%include "Estia_selene2.instr"
 
 /**********************************
  * Optical components within cave *
@@ -299,24 +269,6 @@ COMPONENT ac_slit = Slit(
 /***************
  * Sample area *
  ***************/
-    COMPONENT tof_sample = Monitor_nD(
-        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) RELATIVE arm_sample_beam
-        ROTATED (0, 0, 0) RELATIVE arm_sample_beam
-                
-
-/* NiTi multilayer sample */
-COMPONENT sample = Mirror(
-    xwidth = sample_length, yheight = sample_height,
-    center = 1, transmit = 0,
-    reflect = "NiTiML.ref"
-    )
-    WHEN sample==0
-    AT (0, 0, 0) RELATIVE arm_sample
-    ROTATED (0, 90, 0) RELATIVE arm_sample
 
 /* ideal reflector as reference */
 COMPONENT reference_sample = Mirror(
@@ -328,88 +280,28 @@ COMPONENT reference_sample = Mirror(
     AT (0, 0, 0) RELATIVE arm_sample
     ROTATED (0, 90, 0) RELATIVE arm_sample
 
-/* Nickel film on silicon */
-COMPONENT ni_sample = Mirror(
+/* User defined sample file */
+COMPONENT mooc_sample = Mirror(
     xwidth = sample_length, yheight = sample_height,
     center = 1, transmit = 0,
-    reflect = "Si-Ni.ref"
+    reflect = sample_file
     )
     WHEN sample==2
     AT (0, 0, 0) RELATIVE arm_sample
     ROTATED (0, 90, 0) RELATIVE arm_sample
 
-/* Silicon with natural oxide */
-COMPONENT si_sample = Mirror(
-    xwidth = sample_length, yheight = sample_height,
-    center = 1, transmit = 0,
-    reflect = "Si-SiO2.ref"
-    )
-    WHEN sample==3
-    AT (0, 0, 0) RELATIVE arm_sample
-    ROTATED (0, 90, 0) RELATIVE arm_sample
 
 
-
-
-
-/* Rotate spin 90 degrees*/
-COMPONENT arm_hack = Arm()
-       AT (0, 0, 0) RELATIVE arm_sample
-       EXTEND %{
-         double tmp = sy;
-         sy = sx;
-         sx = -tmp;
-       %}
-
-
-
-
-/* flipped arm detector to position analyzer correctly for beam path 1*/
-COMPONENT arm_detector2 = Arm()
-    AT (0, 0, 0) RELATIVE arm_detector
-  ROTATED (0,0,180) RELATIVE arm_detector
-
-
-
-COMPONENT arm_analyzer = Arm()
-    AT (0, 0, 0) RELATIVE arm_detector2
-    ROTATED (-selene_theta+(Theta1_analyzer1-Theta2_analyzer1)/2.0, 0, 0) RELATIVE arm_detector2
-
-COMPONENT virtual_analyzer_flipper = Arm() // Gone -> Pol_SF_ideal(ny=1, xwidth=1, yheight=1, zdepth=0.001)
-        WHEN (enable_analyzer>2)
-        AT (0, 0, analyzer_flipper_start) RELATIVE arm_analyzer
-        ROTATED (0,0,0.0) RELATIVE arm_analyzer
-    EXTEND %{
-        if(INSTRUMENT_GETPAR(enable_analyzer)>2) sx *=-1;
-    %}
-
-COMPONENT arm_analyzer2 = Arm()
-    AT (0, 0, 0) RELATIVE arm_detector2
-    ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector2
-
-/* polarization analyser */
-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(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
-        AT (0, 0.0, analyzer2_start) RELATIVE arm_analyzer2
-        ROTATED (0,0,0.0) RELATIVE arm_analyzer2
-
 /* 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 sy sz user1, list all",
+        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",
         xwidth = 0.5, yheight = 1.0)
         WHEN sample!=4
         AT (0, 0, detector_arm+0.00001) RELATIVE arm_detector
-        ROTATED (0, 0, 0) 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
 
 /***********************************************************************/
 

simulation/Estia_feeder.instr

@@ -144,7 +144,7 @@ COMPONENT NBOA_window=Al_window(thickness=NBOA_Al_entrance_length)
     AT (0,0,NBOA_Al_entrance_start) RELATIVE ISCS
 
     
-COMPONENT NBOA_side = Absorber(xmin=NBOA_side_x, xmax=NBOA_side_x+0.012,
+COMPONENT NBOA_side = AbsorberFixed(xmin=NBOA_side_x, xmax=NBOA_side_x+0.012, 
                                       ymin=-0.2, ymax=0.2, 
                                       zmin=0.0, zmax=E02_01_01_Cu_length+E02_01_01_length+E02_01_02_length+E02_01_03_length)
     AT (0, 0, E02_01_01_Cu_start) RELATIVE ISCS
@@ -181,7 +181,7 @@ COMPONENT NBOA_Cu_collimator=Guide_gravity(
 * feeder neutron guide *
 ************************/
 
-COMPONENT E02_01_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_011_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_011_wh/2.0, ymax=E02_01_011_wh/2.0, 
                                       zmin=0.0, zmax=E02_01_01_length-0.5)
     AT (0, 0, E02_01_01_start) RELATIVE ISCS
@@ -190,7 +190,7 @@ COMPONENT E02_01_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_01_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_012_wh/2.0, ymax=E02_01_012_wh/2.0, 
                                       zmin=0.5-E02_01_01_Cu_length, zmax=E02_01_01_length)
     AT (0, 0, E02_01_01_start) RELATIVE ISCS
@@ -210,7 +210,7 @@ COMPONENT E02_01_01 = Elliptic_guide_gravity(
     enableGravity=enable_gravity)
     AT (0, 0, E02_01_01_start) RELATIVE ISCS
 
-COMPONENT E02_01_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_021_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_021_wh/2.0, ymax=E02_01_021_wh/2.0, 
                                       zmin=0.0, zmax=E02_01_02_length/2.0)
     AT (0, 0, E02_01_02_start) RELATIVE ISCS
@@ -219,7 +219,7 @@ COMPONENT E02_01_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_01_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_022_wh/2.0, ymax=E02_01_022_wh/2.0, 
                                       zmin=E02_01_02_length/2.0, zmax=E02_01_02_length)
     AT (0, 0, E02_01_02_start) RELATIVE ISCS
@@ -238,7 +238,7 @@ COMPONENT E02_01_02 = Elliptic_guide_gravity(
     enableGravity=enable_gravity)
     AT (0, 0, E02_01_02_start) RELATIVE ISCS
 
-COMPONENT E02_01_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_031_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_031_wh/2.0, ymax=E02_01_031_wh/2.0, 
                                       zmin=0.0, zmax=0.5)
     AT (0, 0, E02_01_03_start) RELATIVE ISCS
@@ -247,7 +247,7 @@ COMPONENT E02_01_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_01_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_032_wh/2.0, ymax=E02_01_032_wh/2.0, 
                                       zmin=0.5, zmax=1.0)
     AT (0, 0, E02_01_03_start) RELATIVE ISCS
@@ -256,7 +256,7 @@ COMPONENT E02_01_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_01_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_01_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_01_033_wh/2.0, ymax=E02_01_033_wh/2.0, 
                                       zmin=1.0, zmax=E02_01_03_length)
     AT (0, 0, E02_01_03_start) RELATIVE ISCS
@@ -299,7 +299,7 @@ COMPONENT NFGA_Al_window_in=Al_window(
     WHEN enable_windows
     AT (0,0,NFGA_Al_start) RELATIVE ISCS
 
-COMPONENT NFGA_side1 = Absorber(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
+COMPONENT NFGA_side1 = AbsorberFixed(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008, 
                                       ymin=-0.2, ymax=0.2, 
                                       zmin=0.002, zmax=0.502)
     AT (0, 0, E02_02_01_start-0.002) RELATIVE ISCS
@@ -309,7 +309,7 @@ COMPONENT NFGA_side1 = Absorber(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT NFGA_side2 = Absorber(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
+COMPONENT NFGA_side2 = AbsorberFixed(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008, 
                                       ymin=-0.2, ymax=0.2, 
                                       zmin=0.502, zmax=E02_02_01_length+E02_02_02_length+E02_02_03_length+0.002)
     AT (0, 0, E02_02_01_start-0.002) RELATIVE ISCS
@@ -321,7 +321,7 @@ COMPONENT NFGA_side2 = Absorber(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
     %}
         
 // in-bunker feeder segments	
-COMPONENT E02_02_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_011_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_011_wh/2.0, ymax=E02_02_011_wh/2.0, 
                                       zmin=0.0, zmax=E02_02_01_length/2.0)
     AT (0, 0, E02_02_01_start) RELATIVE ISCS
@@ -330,7 +330,7 @@ COMPONENT E02_02_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_02_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_012_wh/2.0, ymax=E02_02_012_wh/2.0, 
                                       zmin=E02_02_01_length/2.0, zmax=E02_02_01_length)
     AT (0, 0, E02_02_01_start) RELATIVE ISCS
@@ -349,7 +349,7 @@ COMPONENT E02_02_01 = Elliptic_guide_gravity(
     enableGravity=enable_gravity)
     AT (0, 0, E02_02_01_start) RELATIVE ISCS
     
-COMPONENT E02_02_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_021_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_021_wh/2.0, ymax=E02_02_021_wh/2.0, 
                                       zmin=0.0, zmax=E02_02_02_length/2.0)
     AT (0, 0, E02_02_02_start) RELATIVE ISCS
@@ -358,7 +358,7 @@ COMPONENT E02_02_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_02_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_022_wh/2.0, ymax=E02_02_022_wh/2.0, 
                                       zmin=E02_02_02_length/2.0, zmax=E02_02_02_length)
     AT (0, 0, E02_02_02_start) RELATIVE ISCS
@@ -377,7 +377,7 @@ COMPONENT E02_02_02 = Elliptic_guide_gravity(
     enableGravity=enable_gravity)
     AT (0, 0, E02_02_02_start) RELATIVE ISCS
     
-COMPONENT E02_02_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_031_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_031_wh/2.0, ymax=E02_02_031_wh/2.0, 
                                       zmin=0.0, zmax=0.5)
     AT (0, 0, E02_02_03_start) RELATIVE ISCS
@@ -386,7 +386,7 @@ COMPONENT E02_02_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_02_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_032_wh/2.0, ymax=E02_02_032_wh/2.0, 
                                       zmin=0.5, zmax=1.0)
     AT (0, 0, E02_02_03_start) RELATIVE ISCS
@@ -395,7 +395,7 @@ COMPONENT E02_02_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
         ALLOW_BACKPROP;    
         PROP_Z0;
     %}
-COMPONENT E02_02_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
+COMPONENT E02_02_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0, 
                                       ymin=-E02_02_033_wh/2.0, ymax=E02_02_033_wh/2.0, 
                                       zmin=1.0, zmax=E02_02_03_length)
     AT (0, 0, E02_02_03_start) RELATIVE ISCS

simulation/Estia_mf.instr

@@ -1,242 +0,0 @@
-/*******************************************************************************
-*         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,
-                                      rUpData="SNAG_m5p5_Tup.ref", rDownData="SNAG_m5p5_Tdown.ref")
-		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
-		%{
-			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,
-                                      rUpData="SNAG_m5p5_Tup.ref", rDownData="SNAG_m5p5_Tdown.ref")
-		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(nb=21, ndiv=41, filename="mf_hordiv",
-          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0)
-		  AT (0, 0, 0) RELATIVE arm_polarizer
-
-COMPONENT mf_divergence_v = DivPos_monitor(nb=21, ndiv=41, filename="mf_verdiv",
-          xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0, vertical=1)
-		  AT (0, 0, 0) 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

@@ -1,204 +0,0 @@
-/*******************************************************************************
-*         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

@@ -0,0 +1,204 @@
+/*******************************************************************************
+*         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 *
+**************************************/  
+
+
+
+/**************************
+* 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
+
+/* 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

@@ -1,270 +0,0 @@
-/*******************************************************************************
-*         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

@@ -1,210 +0,0 @@
-/*******************************************************************************
-*         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/SNAG_m5p5_Tdown.ref

@@ -1,58 +0,0 @@
-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

@@ -1,58 +0,0 @@
-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/ScanningSlit.comp

@@ -0,0 +1,93 @@
+/*******************************************************************************
+*
+* McStas, neutron ray-tracing package
+*         Copyright (C) 1997-2011, All rights reserved
+*         Risoe National Laboratory, Roskilde, Denmark
+*         Institut Laue Langevin, Grenoble, France
+*
+* Component: ScanningSlit
+*
+* %I
+* Written by: Jochen Stahn, based on the Slip component by
+*             Kim Lefmann and Henrik M. Roennow
+* Date: 15. 10. 2013
+* Version: $Revision: 0.00 $
+* Origin: PSI
+* Release: McStas 1.12c
+*
+* moving rectangular slit with optional insignificance cut
+*
+* %D
+* A simple rectangular slit, where the blades in x-direction move with time.
+* No transmission around the slit is allowed.
+* If cutting option is used, low-weight neutron rays are ABSORBED
+*
+* Example: FastSlit(height=0.1)
+*
+* %P
+* INPUT PARAMETERS
+*
+* height: height of the opening in y, centered.
+*
+* Optional parameters:
+* cut:  Lower limit for allowed weight (1)
+*
+* %E
+*******************************************************************************/
+
+
+DEFINE COMPONENT ScanningSlit
+DEFINITION PARAMETERS ()
+SETTING PARAMETERS (xmin=0, xmax=0, ymin=0, ymax=0, height=0.1, dist=2.4, lmin=4, lmax=10, omeg=0, eps=1.25, reso=0.01, dthet=1.5, posi=33.7)
+OUTPUT PARAMETERS ()
+/* STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) */
+
+DECLARE
+ %{
+  double lact;
+ %}
+
+
+INITIALIZE
+%{
+  if (height > 0) { ymax=height/2; ymin=-ymax; }
+  lmax  *= (1+reso/2) ;
+  lmin  *= (1-reso/2) ;
+  eps   *= PI/180 ;
+  omeg  *= PI/180 ;
+  dthet *= PI/180;
+%}
+
+TRACE
+%{
+    PROP_Z0;
+    lact = 3.956e3*t/posi ;
+    xmin = -dist * (1+reso/2)*( eps + 0.5*dthet - dthet*(lact-lmin)/(lmax-lmin) ) ;
+    xmax = -dist * (1-reso/2)*( eps + 0.5*dthet - dthet*(lact-lmin)/(lmax-lmin) ) ;
+    if ((x<xmin || x>xmax || y<ymin || y>ymax))
+      ABSORB;
+    else
+      SCATTER;
+%}
+
+MCDISPLAY
+%{
+  magnify("xy");
+    double xw, yh;
+    xw = (xmax - xmin)/2.0;
+    yh = (ymax - ymin)/2.0;
+    multiline(3, xmin-xw, (double)ymax, 0.0,
+              (double)xmin, (double)ymax, 0.0,
+              (double)xmin, ymax+yh, 0.0);
+    multiline(3, xmax+xw, (double)ymax, 0.0,
+              (double)xmax, (double)ymax, 0.0,
+              (double)xmax, ymax+yh, 0.0);
+    multiline(3, xmin-xw, (double)ymin, 0.0,
+              (double)xmin, (double)ymin, 0.0,
+              (double)xmin, ymin-yh, 0.0);
+    multiline(3, xmax+xw, (double)ymin, 0.0,
+              (double)xmax, (double)ymin, 0.0,
+              (double)xmax, ymin-yh, 0.0);
+%}
+
+END

simulation/Selene_geometry.instr

@@ -0,0 +1,145 @@
+/*******************************************************************************
+*         McStas instrument definition URL=http://www.mcstas.org
+*
+* Instrument: Selene_geometry
+*
+* %Identification
+* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
+* Date: 01. 11. 2017
+* Origin: PSI
+* Release: McStas 2.4.1
+* Version: 1.0
+* %INSTRUMENT_SITE: Estia (E02)
+*
+* Estia is a vertical sample, focusing reflectometer for small sample
+*
+* %Description
+* This is a test model to analyze the impact of the displacement of the two Selene guide
+* mirrors to each other and the virtual source. With this it is possible to evaluate
+* necessary engineering requirements for the guide carriers.
+*
+* %Parameters
+* 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
+* omegaa: [deg]  sample rotation omega
+*
+* %End
+*******************************************************************************/
+
+DEFINE INSTRUMENT Selene_geometry(double sample_length=0.01, double sample_height=0.01, omegaa=1.2,
+        double tx_1=0.0, double tz_1=0.0, double rz_1=0.0, double ry_1=0.0,
+        double tx_2=0.0, double tz_2=0.0, double rz_2=0.0, double ry_2=0.0
+        )
+
+DECLARE
+%{
+
+%}
+
+INITIALIZE
+%{
+
+%}
+TRACE
+
+COMPONENT origin = Progress_bar()
+    AT (0,0,0) ABSOLUTE
+
+COMPONENT ISCS = Arm()
+    AT (0, 0, 0) RELATIVE origin    
+    
+COMPONENT arm_selene1 = Arm()
+    AT (tx_1, 0, tz_1) RELATIVE ISCS
+    ROTATED (0, ry_1, rz_1) RELATIVE ISCS
+    
+COMPONENT arm_selene2 = Arm()
+    AT (tx_2, 0, 2*selene_c+tz_2) RELATIVE ISCS    
+    ROTATED (0, ry_2, rz_2) RELATIVE ISCS
+
+COMPONENT arm_virtual_source_beam = Arm()
+    AT (0, 0, 0) RELATIVE arm_selene1
+    ROTATED (0, 1.25, 0) RELATIVE ISCS
+// Axes used for the virtual source masks, this depends on the sample angle to achieve the correct footprint
+COMPONENT arm_virtual_source = Arm()
+    AT (0, 0, 0) RELATIVE arm_virtual_source_beam
+    ROTATED (0, omegaa, 0) RELATIVE arm_virtual_source_beam
+    
+    
+COMPONENT source = Moderator(radius = 0.01, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002, 
+                                       target_index=3, Emin=0.75, Emax=6)
+    AT (0,0,-0.1) RELATIVE ISCS
+    ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
+
+/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
+COMPONENT virtual_source_TR = Slit(
+    xmin = 0.0, xmax = 1.0, ymin = -1.0, ymax = sample_height/2)
+    AT (0, 0, -0.5*sample_length) 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, ymax = 1.0)
+    AT (0, 0, 0.5*sample_length) RELATIVE arm_virtual_source
+
+%include "Estia_selene.instr"
+
+COMPONENT M_n3 = Monitor_nD(
+        filename = "M_n3", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c-0.01) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+
+COMPONENT M_n2 = Monitor_nD(
+        filename = "M_n2", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c-0.005) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+    
+COMPONENT M_n1 = Monitor_nD(
+        filename = "M_n1", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c-0.0025) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+    
+COMPONENT M_0 = Monitor_nD(
+        filename = "M_0", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+
+COMPONENT MFP = Monitor_nD(
+        filename = "MFP", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c+0.0001) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa-90, 0) RELATIVE ISCS
+
+COMPONENT M_p1 = Monitor_nD(
+        filename = "M_p1", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c+0.0025) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+
+COMPONENT M_p2 = Monitor_nD(
+        filename = "M_p2", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c+0.005) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+
+COMPONENT M_p3 = Monitor_nD(
+        filename = "M_p3", 
+        options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
+        xwidth = 0.5, yheight = 0.5)
+    AT (0, 0, 4*selene_c+0.01) RELATIVE ISCS
+    ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
+
+FINALLY
+%{
+%}
+
+END

simulation/compile_if_needed.sh

@@ -1,13 +1,9 @@
 #!/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_selene1.instr -nt Estia_baseline.out ] \
-   || [ Estia_mf.instr -nt Estia_baseline.out ] \
-   || [ Estia_selene2.instr -nt Estia_baseline.out ]; then
+   || [ Estia_selene.instr -nt Estia_baseline.out ]; then
 	rm Estia_baseline.c Estia_baseline.out
-	mcstas -o Estia_baseline.c Estia_baseline.instr --trace
-	mpicc -O2 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
+	mcstas -o Estia_baseline.c Estia_baseline.instr
+	mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
 fi

simulation/hostnames

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

simulation/loop_simus.sh

@@ -1,4 +0,0 @@
-foreach z ( `seq -0.5 0.125 0.5` )
-  echo $z
-  ./run_simu.sh $z
-end

simulation/run_simu.sh

@@ -1,121 +1,117 @@
 #!/bin/bash
 
 DEST=../results
-ncount=1e9
-use_cores=6
-sample=4
-omega=0.8
-sample_length=0.005
-sample_height=0.01
-lambda_start=3.0
-lambda_end=32.0
-
-###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
-bash compile_if_needed.sh
-
-DESTi=$DEST/brilliance_nowindow_5x10
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
-
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --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 \
-        enable_windows=0 \
-        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0
-
-DESTi=$DEST/brilliance_5x10
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
-
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --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
-
-
-ncount=1e10
-sample_length=0.001
-sample_height=0.001
-
-DESTi=$DEST/brilliance_1x1
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
-
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --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
-
-
-# ###################### Reference and Ni-layer measurement 10x10mm² sample ####################
-ncount=1e10
-sample_length=0.01
+ncount=4e9
+use_cores=${SLURM_NPROCS:-6}
+sample_length=0.05
 sample_height=0.01
 lambda_start=3.25
-lambda_end=12.75
-sample=1
+lambda_end=25.5
 
-omega=1.0
-DESTi=$DEST/reference_10x10_10
+bash compile_if_needed.sh
+echo "Running with $use_cores threads"
+
+# ###################### Reference and sample measurement pulse skipping ####################
+
+#
+#      Reference in pulse skipping mode
+omega=1.5
+DESTi=$DEST/mooc_reference_ps
+echo "Reference, pulse skipping"
+echo "Results saved to $DESTi"
 if [ -e "$DESTi" ]; then
 	rm -r "$DESTi"
 fi
 
 mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --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=1
-       
-#     Ni Sample
-        
-ncount=6e9
-sample=2
-omega=0.8
-DESTi=$DEST/nickle_10x10_08
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
+        --dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
+        omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
+        sample=1 sample_length=$sample_length sample_height=$sample_height \
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
+#        
+#     Sample models in pulse skipping mode
+for MODEL in 0 1 2 3 4 5
+do
+  omega=1.5
+  DESTi="$DEST/mooc_model_ps_$MODEL"
+  echo "Model $MODEL, pulse skipping"
+  echo "Results saved to $DESTi"
+  if [ -e "$DESTi" ]; then
+          rm -r "$DESTi"
+  fi
 
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
-        omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
-        sample=$sample sample_length=$sample_length sample_height=$sample_height \
-        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
+  mpirun  -np $use_cores Estia_baseline.out \
+          --dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus  \
+          omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
+          sample=2 sample_length=$sample_length sample_height=$sample_height \
+          lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3 \
+          sample_file="alessandra_model_2_nores_00$MODEL.dat"
 
-ncount=2e9
-omega=3.0
-DESTi=$DEST/nickle_10x10_30
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
+done
 
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
-        omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
-        sample=$sample sample_length=$sample_length sample_height=$sample_height \
-        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
-        
-        
-ncount=1e9
-omega=8.0
-DESTi=$DEST/nickle_10x10_80
-if [ -e "$DESTi" ]; then
-	rm -r "$DESTi"
-fi
+#        
+#     Sample models in pulse skipping mode, larger angle
+for MODEL in 0 1 2 3 4 5
+do
+  omega=4.0
+  DESTi="$DEST/mooc_model_ps2_$MODEL"
+  echo "Model $MODEL, pulse skipping"
+  echo "Results saved to $DESTi"
+  if [ -e "$DESTi" ]; then
+          rm -r "$DESTi"
+  fi
 
-mpirun  -np $use_cores Estia_baseline.out \
-        --dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
-        omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
-        sample=$sample sample_length=$sample_length sample_height=$sample_height \
-        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
+  mpirun  -np $use_cores Estia_baseline.out \
+          --dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus  \
+          omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
+          sample=2 sample_length=$sample_length sample_height=$sample_height \
+          lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3 \
+          sample_file="alessandra_model_2_nores_00$MODEL.dat"
 
-# ###################### Reference and Ni-layer measurement 1x1mm² sample ####################
+done
+
+# ###################### Reference and samples in conventional mode ####################
 # 
 # 
+#
+#      Reference in conventional mode
+ncount=4e9
+lambda_end=13.5
+theta_resolution=0.1
+
+omega=1.5
+DESTi=$DEST/mooc_reference_tof
+echo "Reference, omega=$omega, Delta omega=$theta_resolution"
+echo "Results saved to $DESTi"
+if [ -e "$DESTi" ]; then
+	rm -r "$DESTi"
+fi
+
+mpirun  -np $use_cores Estia_baseline.out \
+        --dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus  \
+        omegaa=$omega operationmode=1 theta_resolution=$theta_resolution over_illumination=0.0002 \
+        sample=1 sample_length=$sample_length sample_height=$sample_height \
+        lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
+#        
+#     Sample models in pulse skipping mode
+for MODEL in 0 1 2 3 4 5
+do
+  for omega in 0.3 0.75 1.875 4.6875 11.71875
+  do
+    ncount=$(bc <<< "5000000000/$omega/$omega")
+    theta_resolution=$(bc <<< "scale=10;0.05*$omega")
+    DESTi="$DEST/mooc_model_tof_$MODEL-$omega"
+    echo "Model $MODEL, omega=$omega, Delta omega=$theta_resolution, n=$ncount"
+    echo "Results saved to $DESTi"
+    if [ -e "$DESTi" ]; then
+            rm -r "$DESTi"
+    fi
+
+    mpirun  -np $use_cores Estia_baseline.out \
+            --dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus  \
+            omegaa=$omega operationmode=1 theta_resolution=$theta_resolution over_illumination=0.0002 \
+            sample=2 sample_length=$sample_length sample_height=$sample_height \
+            lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1 \
+            sample_file="alessandra_model_2_nores_00$MODEL.dat"
+  done
+done

simulation/run_simu.tsh

@@ -0,0 +1,23 @@
+#!/bin/tcsh
+#SBATCH -J mooc_McEstia
+#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`"
+
+
+module load mcstas
+
+bash run_simu.sh