ESS-Estia/Estia-McStas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

merge into: ESS-Estia/Estia-McStas:polarizer
Branches Tags
ESS-Estia/Estia-McStas:master ESS-Estia/Estia-McStas:paper ESS-Estia/Estia-McStas:polarizer ESS-Estia/Estia-McStas:solid_liquid ESS-Estia/Estia-McStas:mcnp_reference ESS-Estia/Estia-McStas:liquids_addon ESS-Estia/Estia-McStas:shielding ESS-Estia/Estia-McStas:mooc ESS-Estia/Estia-McStas:ersa_wfm
..
pull from: ESS-Estia/Estia-McStas:paper
Branches Tags
ESS-Estia/Estia-McStas:paper ESS-Estia/Estia-McStas:master ESS-Estia/Estia-McStas:polarizer ESS-Estia/Estia-McStas:solid_liquid ESS-Estia/Estia-McStas:mcnp_reference ESS-Estia/Estia-McStas:liquids_addon ESS-Estia/Estia-McStas:shielding ESS-Estia/Estia-McStas:mooc ESS-Estia/Estia-McStas:ersa_wfm

4 Commits
polarizer .. paper

Author SHA1 Message Date
glavic_a 118b9f5748 start building new setup for paper 2025-10-24 16:16:14 +02:00
glavic_a bf97ff488e Also update Polarizer.comp 2025-10-24 15:09:43 +02:00
glavic_a c85792ef22 Update model from DMSC version for McStas 3.x 2025-10-24 14:48:27 +02:00
glavic_a 6e0d611e73 add pycharm to gitignore 2025-10-24 14:37:09 +02:00
27 changed files with 788 additions and 1752 deletions
Expand all files Collapse all files
.gitignore
+1
View File
@@ -8,6 +8,7 @@ old
*.backup
.project
.settings
.idea
*.pyc
*.pyo
hosts
analysis/estia_analyze.py
+1 -1
View File
@@ -8,7 +8,7 @@ import sys, os
from numpy import *
seterr(all='ignore')
import estia_help as eh
import mcstas_reader as mr
from aglib import mcstas_reader as mr
# Scaling factor of source monitor after normalizing by area to get to brilliance.
analysis/mcstas_reader.py
-422
View File
@@ -1,422 +0,0 @@
#-*- 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 list(self.info['data'].keys())
def monitors(self):
output=[]
for val in list(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 list(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 list(self.keys()):
data=self.data_loader.load_item(item)
self._data[item]=data
return data
else:
raise KeyError("Can't find dataset %s"%item)
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 in ['Li', 'List'] and y_col=='p': # Detector_nD
cols=item_info['variables'].split()
data=loadtxt(fname, dtype={'names': cols, 'formats': ['f4']*len(cols)})
return TofData(data, item_info)
else:
raw=loadtxt(fname)
data=raw[:len(raw)//3]
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 list(self.hdf['data'].keys()):
node=self.hdf['data/'+item]
info={}
for key, value in list(node.attrs.items()):
info[key.strip()]=value.strip().decode('ascii')
if info['filename'][-4:]=='.dat' or '_list.' in info['filename']:
self.info['data'][info['component']]=info
else:
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 in ['Li', 'List'] 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=list(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=list(map(float, self.info['xylimits'].split()))
if log:
ax.imshow(self.data, origin='lower', extent=limits, aspect='auto', norm=LogNorm())
else:
ax.imshow(self.data, origin='lower', extent=limits, aspect='auto')
ax.set_xlabel(self.info['xlabel'])
ax.set_ylabel(self.info['ylabel'])
ax.set_title(self.info['component'])
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, str):
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, str):
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
+52 -18
View File
@@ -79,8 +79,6 @@ double selene1_center;
double selene1_shift;
double selene1_rot;
int p_int=0; // a flag that gets incremented if a polarizer mirror scatters
// Selene 2
@@ -107,12 +105,20 @@ double velocity_max ; // m/s neutron velocity of lambda_min
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
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
@@ -121,7 +127,7 @@ INITIALIZE
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;
/* print out some calculated parameter for checking purposes */
printf(" Chopper phase = %.1f deg\n", chopper_phase);
@@ -237,17 +243,15 @@ COMPONENT chopper = DiskChopper(radius=chopper_diameter/2.0, yheight=0.02,
AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
COMPONENT vs_divergence_h = DivPos_monitor(nh=21, ndiv=41, filename="vs_hordiv",
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0,
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(nh=21, ndiv=41, filename="vs_verdiv",
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0,
restore_neutron=1)
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
ROTATED (0,0,90) RELATIVE arm_selene1
/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
COMPONENT virtual_source_TR = Slit(
@@ -345,13 +349,43 @@ COMPONENT si_sample = Mirror(
ROTATED (0, 90, 0) RELATIVE arm_sample
COMPONENT arm_analyzer = Arm()
/* 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 (-selene_theta+(Theta1_analyzer1-Theta2_analyzer1)/2.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_detector
ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector
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,
@@ -365,12 +399,12 @@ COMPONENT analyzer2 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, r
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_analyzer
ROTATED (0,0,0.0) RELATIVE arm_analyzer
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", user1=p_int,
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",
xwidth = 0.5, yheight = 1.0)
WHEN sample!=4
simulation/Estia_feeder.instr
+17 -17
View File
@@ -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 = AbsorberFixed(xmin=NBOA_side_x, xmax=NBOA_side_x+0.012,
COMPONENT NBOA_side = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_011_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_01_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_012_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_021_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_01_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_022_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_031_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_01_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_032_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_01_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_01_033_wedge = Absorber(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 = AbsorberFixed(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
COMPONENT NFGA_side1 = Absorber(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 = AbsorberFixed(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT NFGA_side2 = AbsorberFixed(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
COMPONENT NFGA_side2 = Absorber(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 = AbsorberFixed(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
%}
// in-bunker feeder segments
COMPONENT E02_02_011_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_011_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_02_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_012_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_021_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_02_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_022_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_031_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_02_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_032_wedge = Absorber(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 = AbsorberFixed(xmin=-0.1, xmax=0.0,
ALLOW_BACKPROP;
PROP_Z0;
%}
COMPONENT E02_02_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
COMPONENT E02_02_033_wedge = Absorber(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
+6 -10
View File
@@ -145,15 +145,13 @@ COMPONENT arm_pol2 = Arm()
COMPONENT polarizer2_lin = Pol_mirror(zwidth = pol2_lin_length, yheight = pol_hfar, p_reflect=0,
rUpFunc=TableReflecFunc, rUpPar="SNAG_m5p5_Tup.ref",
rDownFunc=TableReflecFunc, rDownPar="SNAG_m5p5_Tdown.ref", useTables=1)
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
%{
p_int=0; // reset for each trace
if (SCATTERED==2) {
p_int +=2;
}
@@ -185,8 +183,7 @@ COMPONENT replacement_pol2_set = Slit(xwidth=0.2, yheight=0.2)
GROUP polarizer2_set
COMPONENT polarizer1_lin = Pol_mirror(zwidth = pol1_lin_length, yheight = pol_hfar, p_reflect=0,
rUpFunc=TableReflecFunc, rUpPar="SNAG_m5p5_Tup.ref",
rDownFunc=TableReflecFunc, rDownPar="SNAG_m5p5_Tdown.ref", useTables=1)
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
@@ -223,14 +220,13 @@ COMPONENT replacement_pol1_set = Slit(xwidth=0.2, yheight=0.2)
COMPONENT mf_divergence_h = DivPos_monitor(nh=21, ndiv=41, filename="mf_hordiv",
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0)
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(nh=21, ndiv=41, filename="mf_verdiv",
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv_h=2.0)
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
ROTATED (0,0,90) RELATIVE arm_polarizer
simulation/Pol_SF_ideal.comp
-91
View File
@@ -1,91 +0,0 @@
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Written by: Erik B Knudsen
* Date: May 2017
* Version: $Revision$
* Release: McStas 2.4
* Origin: DTU Physics
*
* Component: Pol_RFSF_ideal
*
*
* %I
*
* Ideal model of a spin flipper
*
* %D
* This component simply mirrors the polarization vector of the neutron
* ray in the plane through (0,0,0) with normal nx,ny,nz.
* The flipper is surrounded by a perfectly absorbing box. Neutron rays not hitting
* the box are left untouched.
*
* %P
* Input parameters:
* nx: [ ] x-component of the normal vector of the flipping plane.
* ny: [ ] y-component of the normal vector of the flipping plane.
* nz: [ ] z-component of the normal vector of the flipping plane.
* xwidth: [m] width of the spin flipper.
* yheight: [m] height of the spin flipper.
* zdepth: [m] length of the spin flipper.
*
*
* %E
*******************************************************************************/
DEFINE COMPONENT Pol_SF_ideal
DEFINITION PARAMETERS ()
SETTING PARAMETERS (nx=0, ny=1, nz=0, xwidth=0.1, yheight=0.1, zdepth=0.1)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
TRACE
%{
int hit;
double t0,t1;
hit=box_intersect(&t0,&t1, x,y,z,vx,vy,vz, xwidth,yheight,zdepth);
if(hit){
PROP_DT(t0);
if(fabs(z- -zdepth*0.5)>DBL_EPSILON){
/*neutron must have hit the side walls*/
ABSORB;
}
/*move to center of box and flip*/
PROP_Z0;
SCATTER;
double s=scalar_prod(sx,sy,sz,nx,ny,nz);
if (s!=0){
sx-=2*s*nx;
sy-=2*s*ny;
sz-=2*s*nz;
}
PROP_DT((t1-t0)/2);/*propagate the remaining distance to the box exit*/
if(fabs(z-zdepth*0.5)>DBL_EPSILON){
ABSORB;
}
}
%}
MCDISPLAY
%{
double dx=xwidth/16;
double dy=yheight/8;
box(0,0,0,xwidth,yheight,zdepth);
line(dx,-dy,0,dx,-dy+yheight/2.0,0);
line(-dx,dy,0,-dx,dy-yheight/2.0,0);
line(dx,-dy+yheight/2.0,0, dx+xwidth/16,-dy+yheight-yheight/16,0);
line(dx,-dy+yheight/2.0,0, dx-xwidth/16,-dy+yheight-yheight/16,0);
line(-dx,dy-yheight/2.0,0, -dx+xwidth/16,dy-yheight+yheight/16,0);
line(-dx,dy-yheight/2.0,0, -dx-xwidth/16,dy-yheight+yheight/16,0);
%}
END
simulation/Polariser.comp
+160 -95
View File
@@ -48,9 +48,6 @@
**********************************************************************************************************************/
DEFINE COMPONENT Polariser
DEFINITION PARAMETERS ()
SETTING PARAMETERS (int enable_ref = 1,
int abs_out = 1,
int abs_single = 0,
@@ -72,26 +69,6 @@ SETTING PARAMETERS (int enable_ref = 1,
T_loss=4.0e3,
string reflect_d = "ReflectivitySpinDown.dat",
string reflect_u = "ReflectivitySpinUp.dat")
/*
OUTPUT PARAMETERS (dDThetaHalf,
dZmin,
dZmax,
ab,
bb,
dOffsetY,
qc_substrate,
nRflctUp,
nRflctDn,
dGAx,
dGAy,
dGAz)
DECLARE
%{
%}
*/
SHARE
%{
@@ -100,44 +77,19 @@ SHARE
%include "read_table-lib"
%include "ref-lib"
/**********************************************************************************************************************
Declaration of the variables used in the McStas method Polariser */
t_Table* pTableDn=NULL; /* The table for the spin-down reflectivity */
t_Table* pTableUp=NULL; /* The table for the spin-up reflectivity */
int debug_print=0; /* For debugging purposes */
int abs_reason[]={0,0,0,0,0,0,0,0}; /* Auxiliary memory used for flagging the absorbed neutrons */
/* abs_reason[0]: unable to obtain intersection with the first surface of the polariser */
/* abs_reason[1]: unreliable intersection with the first surface of the polariser */
/* abs_reason[2]: out of bounds (in the x direction) intersection with the first surface of the polariser */
/* abs_reason[3]: reflection at the intersection with the first surface of the polariser */
/* abs_reason[4]: unable to obtain intersection with the second surface of the polariser */
/* abs_reason[5]: unreliable intersection with the second surface of the polariser */
/* abs_reason[6]: out of bounds (in the x direction) intersection with the second surface of the polariser */
/* abs_reason[7]: unphysical (null velocity, erroneous refractive index, negative propagation time, non-positive weight, etc.) */
double pdTmp[]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}; /* Auxiliary memory */
int nRflctDn,nRflctUp;
double dTlrnce,dDThetaHalf,qc_Ni,qc_substrate,ab,bb,dOffsetY,dZmin,dZmax,dGAx,dGAy,dGAz;
double dSpdLght,dMassNeut,dMassNeut,dhbarc;
Coords localG; /* Gravitational-acceleration vector */
Coords eCRD_3DU1; /* Constant unit vector (x direction) */
/*
-----------------------------------------------------------------------------------------------------------------------
Definition of functions used in the main body of the method
*/
double FTN1(double dX) {double dTmp=exp(pdTmp[10]*dX); double dSin=sin(dX);
double FTN1(double dX, double* pdTmp) {double dTmp=exp(pdTmp[10]*dX); double dSin=sin(dX);
return dTmp*((pdTmp[0]*dTmp*dSin+pdTmp[1])*dSin+pdTmp[2]*cos(dX))+pdTmp[3];};
double FTN2(double dX) {double dTmp=exp(pdTmp[10]*dX); double dCos=cos(dX);
double FTN2(double dX, double* pdTmp) {double dTmp=exp(pdTmp[10]*dX); double dCos=cos(dX);
return dTmp*((pdTmp[0]*dTmp*dCos+pdTmp[1])*dCos+pdTmp[2]*sin(dX))+pdTmp[3];};
double FTN3(double dX) {double dT=(exp(pdTmp[10]*dX)*(pdTmp[0]*sin(dX)+pdTmp[1]*cos(dX))+pdTmp[2])/pdTmp[3];
double FTN3(double dX, double* pdTmp) {double dT=(exp(pdTmp[10]*dX)*(pdTmp[0]*sin(dX)+pdTmp[1]*cos(dX))+pdTmp[2])/pdTmp[3];
return pdTmp[4]+(pdTmp[5]+pdTmp[6]*dT)*dT+(pdTmp[7]+(pdTmp[8]+pdTmp[9]*dT)*dT)*tan(dX);};
double FindRoot(double (*FUNCTION)(double ),const double dX1,const double dX2,const double dAccuracy,int* pnSuccess);
double FindRoot(int fct, double *pdTmp, const double dX1,const double dX2,const double dAccuracy,int* pnSuccess);
/*
-----------------------------------------------------------------------------------------------------------------------
@@ -145,7 +97,7 @@ Given a function FUNCTION defined in the interval [dX1,dX2], subdivide the inter
segments, and search for zero crossings of the function. The value of pnNoRoots at input is the maximal number of roots
sought, and is at output reset to the number of bracketing pairs pdX1[1,...,pnNoRoots], pdX2[1,...,pnNoRoots] found.
*/
void BracketRoots(double (*FUNCTION)(double ),
void BracketRoots(int fct, double* pdTmp,
const double dX1,
const double dX2,
const int nNoInts,
@@ -157,11 +109,29 @@ void BracketRoots(double (*FUNCTION)(double ),
double dX,dValue1,dValue2,dStep;
nbb=0;
const int nNoIters=100;
dStep=(dX2-dX1)/(double)nNoInts; /* Determine the spacing appropriate to the mesh */
dValue1=(*FUNCTION)(dX=dX1);
if(*pdX1 == *pdX2)*pdX2 += (double)nNoIters;
if (fct==1) {
dValue1=FTN1(dX=dX1,pdTmp);
} else if (fct==2) {
dValue1=FTN2(dX=dX1,pdTmp);
} else if (fct==3) {
dValue1=FTN3(dX=dX1,pdTmp);
} else {
fprintf(stderr, "Sorry, function with index %d\n is not supported!\n",fct);
exit(-1);
}
for(nCnt=0;nCnt<nNoInts;++nCnt) { /* Loop over all intervals */
dValue2=(*FUNCTION)(dX += dStep);
if (fct==1) {
dValue2=FTN1(dX += dStep,pdTmp);
} else if (fct==2) {
dValue2=FTN2(dX += dStep,pdTmp);
} else if (fct==3) {
dValue2=FTN3(dX += dStep,pdTmp);
}
if(dValue2*dValue1 <= 0.0) {
pdX1[nbb]=dX-dStep;
pdX2[nbb++]=dX;
@@ -180,9 +150,7 @@ Given a function FUNCTION and an initial guess range from dX1 to dX2, the method
a root is bracketed by the returned values dX1 and dX2 (in which case the method returns 1) or until the range becomes
unacceptably large (in which case the method returns 0).
*/
int BracketRoot(double (*FUNCTION)(double ),
double* pdX1,
double* pdX2) {
int BracketRoot(int fct, double *pdTmp, double* pdX1, double* pdX2) {
const int nNoIters=100;
const double dFactor=1.6;
@@ -191,14 +159,39 @@ int BracketRoot(double (*FUNCTION)(double ),
double dValue1,dValue2;
if(*pdX1 == *pdX2)*pdX2 += (double)nNoIters;
dValue1=(*FUNCTION)(*pdX1);
dValue2=(*FUNCTION)(*pdX2);
if (fct==1) {
dValue1=FTN1(*pdX1,pdTmp);
dValue2=FTN1(*pdX2,pdTmp);
} else if (fct==2) {
dValue1=FTN2(*pdX1,pdTmp);
dValue2=FTN2(*pdX2,pdTmp);
} else if (fct==3) {
dValue1=FTN3(*pdX1,pdTmp);
dValue2=FTN3(*pdX2,pdTmp);
} else {
fprintf(stderr, "Sorry, function with index %d\n is not supported!\n",fct);
exit(-1);
}
for(nCnt=0;nCnt<nNoIters;++nCnt) {
if(dValue1*dValue2<0.0)return 1;
if(fabs(dValue1)<fabs(dValue2)) {
dValue1=(*FUNCTION)(*pdX1 += dFactor*(*pdX1-*pdX2));
if (fct==1) {
dValue1=FTN1(*pdX1 += dFactor*(*pdX1-*pdX2),pdTmp);
} else if (fct==2) {
dValue1=FTN2(*pdX1 += dFactor*(*pdX1-*pdX2),pdTmp);
dValue2=FTN2(*pdX2,pdTmp);
} else if (fct==3) {
dValue1=FTN3(*pdX1 += dFactor*(*pdX1-*pdX2),pdTmp);
dValue2=FTN3(*pdX2,pdTmp);
}
} else {
dValue2=(*FUNCTION)(*pdX2 += dFactor*(*pdX2-*pdX1));
if (fct==1) {
dValue2=FTN1(*pdX2 += dFactor*(*pdX2-*pdX1),pdTmp);
} else if (fct==2) {
dValue2=FTN2(*pdX2 += dFactor*(*pdX2-*pdX1),pdTmp);
} else if (fct==3) {
dValue2=FTN3(*pdX2 += dFactor*(*pdX2-*pdX1),pdTmp);
}
}
}
@@ -296,7 +289,7 @@ dAccuracy.
*/
#define SHFT_FindRoot(dA,dB,dC) dA=dB;dB=dC;dC=dA;
#define SIGN_FindRoot(dA,dB) ((dB) >= 0.0 ? fabs(dA) : -fabs(dA))
double FindRoot(double (*FUNCTION)(double ),
double FindRoot(int fct, double* pdTmp,
const double dX1,
const double dX2,
const double dAccuracy,
@@ -314,23 +307,54 @@ double FindRoot(double (*FUNCTION)(double ),
double dC=dX2;
double dD=0.0;
double dE=0.0;
double dValue1=(*FUNCTION)(dA);
double dValue2=(*FUNCTION)(dB);
double dValue1;
double dValue2;
if (fct==1) {
dValue1=FTN1(dA,pdTmp);
dValue2=FTN1(dB,pdTmp);
} else if (fct==2) {
dValue1=FTN1(dA,pdTmp);
dValue2=FTN1(dB,pdTmp);
} else if (fct==3) {
dValue1=FTN1(dA,pdTmp);
dValue2=FTN1(dB,pdTmp);
} else {
fprintf(stderr, "Sorry, function with index %d\n is not supported!\n",fct);
exit(-1);
}
*pnSuccess=0;
if(dValue1*dValue2>0.0) { /* Attempt to bracket the minimum */
if(BracketRoot(FUNCTION,&dA,&dB)) {
if(BracketRoot(fct,pdTmp,&dA,&dB)) {
dC=dB;
dValue1=(*FUNCTION)(dA);
dValue2=(*FUNCTION)(dB);
if (fct==1) {
dValue1=FTN1(dA,pdTmp);
dValue2=FTN1(dB,pdTmp);
} else if (fct==2) {
dValue1=FTN2(dA,pdTmp);
dValue2=FTN2(dB,pdTmp);
} else if (fct==3) {
dValue1=FTN3(dA,pdTmp);
dValue2=FTN3(dB,pdTmp);
}
} else {
BracketRoots(FUNCTION,dX1,dX2,nNoIters,pdX1,pdX2,&nNoRoots);
BracketRoots(fct,pdTmp,dX1,dX2,nNoIters,pdX1,pdX2,&nNoRoots);
if(nNoRoots) { /* First root to be considered */
dA=pdX1[0];
dB=dC=pdX2[0];
dValue1=(*FUNCTION)(dA);
dValue2=(*FUNCTION)(dB);
dC=dB;
if (fct==1) {
dValue1=FTN1(dA,pdTmp);
dValue2=FTN1(dB,pdTmp);
} else if (fct==2) {
dValue1=FTN2(dA,pdTmp);
dValue2=FTN2(dB,pdTmp);
} else if (fct==3) {
dValue1=FTN3(dA,pdTmp);
dValue2=FTN3(dB,pdTmp);
}
} else {
return 0.0;
}
@@ -387,7 +411,13 @@ double FindRoot(double (*FUNCTION)(double ),
} else {
dB += SIGN_FindRoot(dTmp,dX);
}
dValue2=(*FUNCTION)(dB);
if (fct==1) {
dValue2=FTN1(dB,pdTmp);
} else if (fct==2) {
dValue2=FTN2(dB,pdTmp);
} else if (fct==3) {
dValue2=FTN3(dB,pdTmp);
}
}
return 0.0;
@@ -460,7 +490,8 @@ void DetermineIntersectionWithDistalPolariser(double* pdNCt, /* INPUT: The posi
/* 1: irregular neutrons to be absorbed (or propagated, if the method is called for the second surface and abs_single is set to 0) (failed cases) */
/* 2: irregular neutrons to be propagated (failed cases) */
double* pdRt, /* The angle theta corresponding to the intersection (only for regular neutrons) */
double* pdT) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
double* pdT,
double* pdTmp,double ab, double bb, int* abs_reason) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
int nSuccess;
@@ -493,10 +524,10 @@ void DetermineIntersectionWithDistalPolariser(double* pdNCt, /* INPUT: The posi
pdTmp[1]=ab*(dVy*dVz-dGAz*(dY-dOffY));
pdTmp[2]=-ab*dVy*dVy;
pdTmp[3]=dVy*dVy*dZ-dVy*dVz*(dY-dOffY)+0.5*dGAz*(dY-dOffY)*(dY-dOffY);
if(FTN1(-dHalfRange)*FTN1(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN1(-dHalfRange,pdTmp)*FTN1(dHalfRange,pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN1,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
*pdRt=FindRoot(1, pdTmp, -dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -508,10 +539,10 @@ void DetermineIntersectionWithDistalPolariser(double* pdNCt, /* INPUT: The posi
pdTmp[1]=ab*(dVy*dVz-dGAy*dZ);
pdTmp[2]=-ab*dVz*dVz;
pdTmp[3]=dVz*dVz*(dY-dOffY)-dVy*dVz*dZ+0.5*dGAy*dZ*dZ;
if(FTN2(-dHalfRange)*FTN2(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN2(-dHalfRange, pdTmp)*FTN2(dHalfRange, pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN2,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
*pdRt=FindRoot(2, pdTmp,-dHalfRange, dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -529,10 +560,10 @@ void DetermineIntersectionWithDistalPolariser(double* pdNCt, /* INPUT: The posi
pdTmp[7]=-dZ;
pdTmp[8]=-dVz;
pdTmp[9]=-0.5*dGAz;
if(FTN3(-dHalfRange)*FTN3(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN3(-dHalfRange, pdTmp)*FTN3(dHalfRange, pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN3,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
*pdRt=FindRoot(3, pdTmp, -dHalfRange, dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -560,7 +591,8 @@ void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The po
/* 1: irregular neutrons to be absorbed (or propagated, if the method is called for the second surface and abs_single is set to 0) (failed cases) */
/* 2: irregular neutrons to be propagated (failed cases) */
double* pdRt, /* The angle theta corresponding to the intersection (only for regular neutrons) */
double* pdT) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
double* pdT,
double* pdTmp, double ab, double bb, int* abs_reason) { /* The time needed for the neutron to reach the intersection (only for regular neutrons) */
int nSuccess;
@@ -593,10 +625,10 @@ void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The po
pdTmp[1]=ab*(dVy*dVz-dGAz*(dY-dOffY));
pdTmp[2]=ab*dVy*dVy;
pdTmp[3]=dVy*dVy*dZ-dVy*dVz*(dY-dOffY)+0.5*dGAz*(dY-dOffY)*(dY-dOffY);
if(FTN1(-dHalfRange)*FTN1(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN1(-dHalfRange,pdTmp)*FTN1(dHalfRange,pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN1,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
*pdRt=FindRoot(1, pdTmp,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN1: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -608,10 +640,10 @@ void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The po
pdTmp[1]=ab*(dGAy*dZ-dVy*dVz);
pdTmp[2]=-ab*dVz*dVz;
pdTmp[3]=dVz*dVz*(dY-dOffY)-dVy*dVz*dZ+0.5*dGAy*dZ*dZ;
if(FTN2(-dHalfRange)*FTN2(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN2(-dHalfRange,pdTmp)*FTN2(dHalfRange,pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN2,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
*pdRt=FindRoot(2, pdTmp, -dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN2: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -629,10 +661,10 @@ void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The po
pdTmp[7]=dZ;
pdTmp[8]=dVz;
pdTmp[9]=0.5*dGAz;
if(FTN3(-dHalfRange)*FTN3(dHalfRange)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(FTN3(-dHalfRange, pdTmp)*FTN3(dHalfRange, pdTmp)>0.0) { /* No root in [-delta_theta/2,delta_theta/2] */
if(nAbsIrr) { abs_reason[nOff] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
*pdRt=FindRoot(&FTN3,-dHalfRange,dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
*pdRt=FindRoot(3, pdTmp,-dHalfRange, dHalfRange,dTol,&nSuccess); /* Root of the function FTN3: angle theta corresponding to the intersection */
if(!nSuccess) { /* Iteration scheme failed: neutron absorbed or propagated depending on nAbsIrr */
if(nAbsIrr) { abs_reason[nOff+1] += 1; *pnType=1; return; } else { *pnType=2; return; }
}
@@ -648,7 +680,7 @@ void DetermineIntersectionWithProximalPolariser(double* pdNCt, /* INPUT: The po
-----------------------------------------------------------------------------------------------------------------------
The method estimates the time needed for the neutron to reach a specific z
*/
void EstimateExitTime(double* pdT) {
void EstimateExitTime(double* pdT, double* pdTmp) {
double dTmp;
@@ -675,6 +707,39 @@ void EstimateExitTime(double* pdT) {
%}
DECLARE
%{
/**********************************************************************************************************************
Declaration of the variables used in the McStas method Polariser */
t_Table* pTableDn; /* The table for the spin-down reflectivity */
t_Table* pTableUp; /* The table for the spin-up reflectivity */
int debug_print; /* For debugging purposes */
int abs_reason[8]; /* Auxiliary memory used for flagging the absorbed neutrons */
double pdTmp[11]; /* Auxiliary memory */
int nRflctDn;
int nRflctUp;
double dTlrnce;
double dDThetaHalf;
double qc_Ni;
double qc_substrate;
double ab;
double bb;
double dOffsetY;
double dZmin;
double dZmax;
double dGAx;
double dGAy;
double dGAz;
double dSpdLght;
double dMassNeut;
double dhbarc;
Coords localG; /* Gravitational-acceleration vector */
Coords eCRD_3DU1; /* Constant unit vector (x direction) */
%}
/**********************************************************************************************************************
Initialisation of variables and constructs used in the McStas method Polariser */
@@ -826,7 +891,7 @@ First surface of the proximal polariser
pdC[0]=0.0;
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the first surface of the proximal polariser */
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT,pdTmp,ab,bb,abs_reason); /* Intersection with the first surface of the proximal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)ABSORB;
@@ -917,10 +982,10 @@ Second surface of the proximal polariser
pdC[0]=dOffsetY;
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the second surface of the proximal polariser */
DetermineIntersectionWithProximalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT,pdTmp,ab,bb,abs_reason); /* Intersection with the second surface of the proximal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)if(abs_single) { ABSORB; } else { break; }
if(nType == 1) { if(abs_single) { ABSORB; } else { break; } }
if(nType == 2)break;
/* A well-defined intersection with the second surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
@@ -1004,7 +1069,7 @@ First surface of the distal polariser
pdC[0]=dOffsetY;
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the first surface of the distal polariser */
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT,pdTmp,ab,bb,abs_reason); /* Intersection with the first surface of the distal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)ABSORB;
@@ -1096,10 +1161,10 @@ Second surface of the distal polariser
pdC[0]=0.0;
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT); /* Intersection with the second surface of the distal polariser */
DetermineIntersectionWithDistalPolariser(pdNCt,pdGA,pdC,&nType,&dRt,&dT,pdTmp,ab,bb,abs_reason); /* Intersection with the second surface of the distal polariser */
/* Failed cases will be absorbed or propagated */
if(nType == 1)if(abs_single) { ABSORB; } else { break; }
if(nType == 1) { if(abs_single) { ABSORB; } else { break; } }
if(nType == 2)break;
/* A well-defined intersection with the second surface of the polariser has been identified; the question whether this intersection corresponds to an appropriate x value remains */
@@ -1175,7 +1240,7 @@ Second surface of the distal polariser
pdTmp[0]=z-(dZmax-dZmin);
pdTmp[1]=vz;
pdTmp[2]=0.5*dGAz;
EstimateExitTime(&dT);
EstimateExitTime(&dT, pdTmp);
if(dT <= 0.0) { abs_reason[7] += 1; ABSORB; } /* Estimation failed */
PROP_DT(dT);
simulation/PolarizerTest.instr
-108
View File
@@ -1,108 +0,0 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: PolarizerTest
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch)
* Date: 25. 07. 2019
* Origin: PSI
* Release: McStas 2.5
* Version: 1.0
*
* Test for polarization parameters of curved transmission polarizer
*
* %Description
*
*
* %Parameters
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT PolarizerTest(m_up=5.08, m_down=0.6, R0=1.00,
m_residual=0.55, Theta=1.5,
lambda_min=2.0, lambda_max=32.0,
gamma=1.66, focus_height=0.00001,
phi=0.0, polarizer_start=0.2
)
DECLARE
%{
double lambda_mean,delta_lambda,polarizer_length,zero_position;
double source_div, detector_height;
int use_streight=0;
%}
INITIALIZE
%{
lambda_mean=(lambda_min+lambda_max)/2.0;
delta_lambda=(lambda_max-lambda_min)/2.0;
polarizer_length=(exp(Theta/gamma)-1.0)*polarizer_start;
zero_position=exp(0.5*Theta/gamma)*polarizer_start;
source_div=0.1*Theta/180.*PI;
detector_height=source_div/0.1*(polarizer_length+polarizer_start+0.5);
%}
TRACE
COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
COMPONENT focus = Arm()
AT (0,0,0) RELATIVE origin
ROTATED (phi,0,0) RELATIVE origin
COMPONENT Source = Source_simple(xwidth = 0.00001, yheight = focus_height,
dist = 0.1, focus_xw = 0.0001,
focus_yh = source_div,
lambda0 = lambda_mean, dlambda = delta_lambda,
gauss=0)
AT (0, 0, 0) RELATIVE origin
COMPONENT FluxIn = L_monitor(xwidth = 0.5, yheight = 0.5, filename="fluxin",
nL=601, Lmin=lambda_min, Lmax=lambda_max)
AT (0, 0, 0.101) RELATIVE origin
COMPONENT polarizer = Polariser(enable_ref=1, abs_ref=0,
reflect_d=0, reflect_u=0, T_loss=4.0e3,
m_substrate=0.469522, d_substrate = 5e-4,
m_u=m_up, m_d=m_down, R0=R0,
m_residual=m_residual,
lin=polarizer_start, length=polarizer_length,
delta_theta=(Theta)*PI/180.0, h2=0.1, h1=0.05,
both_coated=1, alpha=2.3, W = 0.0014)
WHEN use_streight==0
AT (0, -5e-4, polarizer_start) RELATIVE focus
COMPONENT streight = Mirror(xwidth = 0.2, yheight = 0.2,
m=5.0, center=1, transmit=0)
WHEN use_streight==1
AT (0, 0, zero_position) RELATIVE origin
ROTATED (90.0-gamma,0,0) RELATIVE origin
COMPONENT FluxOut = L_monitor(xwidth = 0.5, yheight = detector_height,
filename="fluxout",
nL=601, Lmin=lambda_min, Lmax=lambda_max)
AT (0, 0, polarizer_length+polarizer_start+0.5) RELATIVE origin
COMPONENT PolMon = MeanPolLambda_monitor(xwidth=0.5, yheight=detector_height,
mx=-1,
nL=601, Lmin=lambda_min, Lmax=lambda_max,
filename="polarization")
AT (0, 0, polarizer_length+polarizer_start+0.5) RELATIVE origin
FINALLY
%{
%}
END
simulation/ScanningSlit.comp
-93
View File
@@ -1,93 +0,0 @@
/*******************************************************************************
*
* 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
-145
View File
@@ -1,145 +0,0 @@
/*******************************************************************************
* 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/pol_sweep_analyze.py
-83
View File
@@ -1,83 +0,0 @@
#!/usr/bin/env python
'''
Run a simulation through various values of polarizer 1/2 distance to focus
and generate a datafile with analyzed values for polarization and transmission
for the short and long wavelength area.
'''
import sys, os
from numpy import *
sys.path.append('../analysis')
from mcstas_reader import McSim
Iref=1.0
L=linspace(1.95, 32.05, 302)
Lc=(L[:-1]+L[1:])/2.
gamma1=1.8
gamma2=1.66
command_s_ref="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=1e10 --dir=../results/polarizer_s_ref --format=NeXus sample=1 omegaa=2.0 over_illumination=0.001 lambda_start=2.0 lambda_end=32.0 sample_length=0.01 enable_polarizer=0"
command_l_ref="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=2e9 --dir=../results/polarizer_l_ref --format=NeXus sample=1 omegaa=6.0 over_illumination=0.001 lambda_start=2.0 lambda_end=32.0 sample_length=0.048 enable_polarizer=0"
command_s="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=5e9 --dir=../results/polarizer_s --format=NeXus sample=1 omegaa=2.0 over_illumination=0.001 pol1_start=%%f pol1_angle=%f pol2_start=%%f pol2_angle=%f lambda_start=2.0 lambda_end=32.0 sample_length=0.01 enable_polarizer=3"%(gamma1, gamma2)
command_l="mpiexec -np 60 --hostfile hostnames Estia_baseline.out --ncount=1e9 --dir=../results/polarizer_l --format=NeXus sample=1 omegaa=6.0 over_illumination=0.001 pol1_start=%%f pol1_angle=%f pol2_start=%%f pol2_angle=%f lambda_start=2.0 lambda_end=32.0 sample_length=0.048 enable_polarizer=3"%(gamma1, gamma2)
lengths=linspace(0.1, 0.6, 25)
def get_polarization(ds):
N,ignore=histogram(ds.data['L'], bins=L)
PN,ignore=histogram(ds.data['L'], bins=L, weights=sqrt(ds.data['sx']**2+ds.data['sy']**2+ds.data['sz']**2))
return PN/maximum(N,1)
def get_intensity(ds):
I, ignore=histogram(ds.data['L'], bins=L, weights=ds.data['p'])
return I
def analyze(ds):
P=get_polarization(ds)
I=get_intensity(ds)
T=I/maximum(Iref, Iref[Iref>0].min())
return P,T
def get_data(fname):
det=McSim(fname)['tof_detector']
P,T=analyze(det)
return [P[17], P[230], T[17], T[230]]
def simulate_and_analyze(pol1_start, pol2_start):
print pol1_start, pol2_start
os.system('rm -rf ../results/polarizer_?')
os.system(command_s%(pol1_start, pol2_start))
os.system(command_l%(pol1_start, pol2_start))
global Iref
Iref=Iref_s
ds=get_data('../results/polarizer_s')
Iref=Iref_l
dl=get_data('../results/polarizer_l')
with open('../results/polarizer_data.dat', 'a') as fh:
fh.write((10*'%.6f '+'\n')%tuple([pol1_start, pol2_start]+ds+dl))
if __name__=='__main__':
if not (os.path.exists('../results/polarizer_s_ref') and
os.path.exists('../results/polarizer_l_ref')):
os.system('rm -rf ../results/polarizer_?_ref')
print 'Simulate reference intensities'
os.system(command_s_ref)
os.system(command_l_ref)
ref=McSim('../results/polarizer_s_ref')['tof_detector']
Iref_s=get_intensity(ref)
ref=McSim('../results/polarizer_l_ref')['tof_detector']
Iref_l=get_intensity(ref)
with open('../results/polarizer_data.dat', 'w') as fh:
fh.write('# pol1 pol2 P3.7_s P25_s T3.7_s T25_s P3.7_l P25_l T3.7_l T25_l\n')
for pol1 in lengths:
for pol2 in lengths:
simulate_and_analyze(pol1, pol2)
with open('../results/polarizer_data.dat', 'a') as fh:
fh.write('\n')
simulation/polarizer_coating.sh
-7
View File
@@ -1,7 +0,0 @@
!#/bin/sh
rm PolarizerTest.c PolarizerTest.out
rm -r ../results/polarizer_coating
mcstas -o PolarizerTest.c PolarizerTest.instr --trace
#mpicc -O2 -o PolarizerTest.out PolarizerTest.c -lm -DUSE_MPI
mpicc -O2 -o PolarizerTest.out PolarizerTest.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus -I/usr/local/include/nexus
LD_LIBRARY_PATH=/usr/local/lib/ mpirun -np 6 ./PolarizerTest.out --ncount=1e8 --format=NeXus --dir=../results/polarizer_coating R0=1.0 gamma=1.66
simulation/run_monitor.sh
-12
View File
@@ -1,12 +0,0 @@
#!/bin/bash
mcstas -o Estia_monitor.c Estia_monitor.instr
mpicc -O3 -o Estia_monitor.out Estia_monitor.c -lm -DUSE_MPI
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_ref --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=2
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_trans --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=1
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor100 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.0001
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor001 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.000001
simulation/run_pol.sh
-77
View File
@@ -1,77 +0,0 @@
#!/bin/bash
DEST=../results
ncount=1e8
use_cores=16
sample=4
sample_length=0.005
sample_height=0.01
lambda_start=3.0
lambda_end=32.0
bash compile_if_needed.sh
###################### Reference measurement 10x5mm² and 10x0.5mm² VS ####################
ncount=1e8
sample_length=0.005
sample_height=0.01
DESTi=$DEST/pol_ref_10x5
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores --hostfile hostnames Estia_baseline.out \
--dir="$DESTi" --format=NeXus --ncount=$ncount \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
enable_polarizer=0 enable_analyzer=0
ncount=5e8
sample_length=0.0005
sample_height=0.01
DESTi=$DEST/pol_ref_10x0p5
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores --hostfile hostnames Estia_baseline.out \
--dir="$DESTi" --format=NeXus --ncount=$ncount \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
enable_polarizer=0 enable_analyzer=0
###################### Polarizers 10x5mm² and 10x0.5mm² VS ####################
ncount=2e8
sample_length=0.005
sample_height=0.01
DESTi=$DEST/pol_10x5
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores --hostfile hostnames Estia_baseline.out \
--dir="$DESTi" --format=NeXus --ncount=$ncount \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
enable_polarizer=1 enable_analyzer=0
ncount=1e9
sample_length=0.0005
sample_height=0.01
DESTi=$DEST/pol_10x0p5
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores --hostfile hostnames Estia_baseline.out \
--dir="$DESTi" --format=NeXus --ncount=$ncount \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=0 \
enable_polarizer=1 enable_analyzer=0
simulation/run_selene.py
-52
View File
@@ -1,52 +0,0 @@
#!/usr/bin/env python
#-*- coding: utf8 -*-
import sys, os
from subprocess import call
from numpy import *
from scipy.optimize import leastsq
CALL='/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out '
CALL=+'--dir=../results/selene_geo --ncount=1e9 '
CALL=+'omegaa=1.0 sample=4 sample_length=0.001 sample_height=0.01 '
CALL=+'lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 '
CALL=+'lambda_min=3.75 selene1_foot1y=%.4f selene1_foot2y=%.4f '
def B(x,w):
# box function with full width w
return float32(abs(x)<=(w/2.))
def G(x,I0,x0,sigma):
# Gaussian with intensity I0, center x0 and standard deviation sigma
return I0*exp(-0.5*(x-x0)**2/sigma**2)
def Intensity(x, p):
I0, x0, sigma, w=p
return convolve(B(xc,w), G(xc, I0, x0, sigma), mode='same')
def Beam(x,I0,x0,sigma,w):
I=I0*where((x-x0)<(-w/2.), exp(-0.5*(x-x0+w/2.)**2/sigma**2),
where((x-x0)>(w/2.), exp(-0.5*(x-x0-w/2.)**2/sigma**2), 1.))
return I
def residuals(p, x, y):
I0, x0, sigma, w=p
return y-Beam(x, I0, x0, sigma, w)
def FWHM(pi):
rng=xc[where(Beam(xc, *pi)>=(pi[0]/2.))[0]]
return rng[-1]-rng[0]
x=linspace(-10.0005, 10.0005, 20001)
xc=(x[1:]+x[:-1])/2.
def analyze(fi):
sim=mr.McSim(fi)
det=sim['tof_sample']
ignore,y=det.project1d('x', bins=x/1000.)
p,res=leastsq(residuals, (y.max(), (xc*y).sum()/y.sum(), 0.01, 0.1), (xc, y))
return p,det
if __name__=='__main__':
pass
simulation/run_simu_2.sh
-86
View File
@@ -1,86 +0,0 @@
#!/bin/bash
DEST=../results
use_cores=6
###################### Compile model if necessary ####################
bash compile_if_neede.sh
######## Reference and Ni-layer conventional measurement 10x10mm² sample #############
ncount=1e10
sample_length=0.01
sample_height=0.01
lambda_start=3.25
lambda_min=3.75
lambda_end=12.75
sample=1
omega=1.0
DESTi=$DEST/tof_reference_10x10_10
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=1 theta_resolution=0.01 over_illumination=0.0001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
# Ni Sample
ncount=2e10
sample=2
omega=0.5
DESTi=$DEST/tof_nickle_10x10_05
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=1 theta_resolution=0.03 over_illumination=0.0001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
ncount=1e10
omega=1.2
DESTi=$DEST/tof_nickle_10x10_12
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=1 theta_resolution=0.03 over_illumination=0.0001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
ncount=5e9
omega=3.0
DESTi=$DEST/tof_nickle_10x10_30
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=1 theta_resolution=0.03 over_illumination=0.0001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
ncount=2e9
omega=7.5
DESTi=$DEST/tof_nickle_10x10_75
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=1 theta_resolution=0.03 over_illumination=0.0001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
simulation/run_simu_3.sh
-46
View File
@@ -1,46 +0,0 @@
#!/bin/bash
DEST=../results
ncount=1e9
use_cores=6
sample=4
omega=0.8
sample_length=0.01
sample_height=0.01
lambda_start=3.0
lambda_end=32.0
bash compile_if_neede.sh
################# Reference and Ni-layer measurement 2 pulse skipping ####################
ncount=4e9
sample_length=0.01
sample_height=0.01
lambda_start=3.9
lambda_end=24.0
sample=1
omega=2.0
DESTi=$DEST/single_skip_reference
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.0002 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
sample=2
DESTi=$DEST/single_skip_nickle
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.0002 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
simulation/run_simu_chopper.sh
-119
View File
@@ -1,119 +0,0 @@
#!/bin/bash
DEST=../results
ncount=1e9
use_cores=6
sample=1
omega=15.0
sample_length=0.02
sample_height=0.01
lambda_start=3.5
lambda_min=5.0
lambda_end=15.5
frame_usage=0.92
###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
bash compile_if_neede.sh
###################### Reference and Ni-layer measurement 10x10mm² sample ####################
frame_usage=0.92
DESTi=$DEST/chopper_f1_092
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.94
DESTi=$DEST/chopper_f1_094
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.96
DESTi=$DEST/chopper_f1_096
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.974
DESTi=$DEST/chopper_f1_097s
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.98
DESTi=$DEST/chopper_f1_098
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.99
DESTi=$DEST/chopper_f1_099
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
frame_usage=0.98
lambda_start=2.75
lambda_min=3.75
lambda_end=12.75
DESTi=$DEST/chopper_f1_375A
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
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=0 enable_chopper=1 \
frame_usage=$frame_usage lambda_min=$lambda_min
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
simulation/run_vs.sh
-32
View File
@@ -1,32 +0,0 @@
#!/bin/bash
DEST=../results
ncount=1e8
use_cores=6
sample=1
omega=0.8
sample_length=0.01
sample_height=0.015
if [ Estia_vs.instr -nt Estia_vs.out ]; then
rm Estia_vs.c Estia_vs.out
mcstas -t -o Estia_vs.c Estia_vs.instr
mpicc -O3 -o Estia_vs.out Estia_vs.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
fi
omega=2.0
DESTi=$DEST/vs_data
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_vs.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_gravity=1 enable_chopper=0 source_power=5
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
simulation/sbatch.sh
-40
View File
@@ -1,40 +0,0 @@
#!/bin/bash
#SBATCH -J mcEstia_1
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout.log
#SBATCH -e stderr.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.009 0.001 0.009)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
simulation/sbatch.tsh
-38
View File
@@ -1,38 +0,0 @@
#!/bin/tcsh
#SBATCH -J mcEstia_S
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=2:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout.log
#SBATCH -e stderr.log
#SBATCH --partition=ll_short
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
#module unload intel
#module load intel
module load mcstas
bash compile_if_needed.sh
mpirun -np $SLURM_NPROCS Estia_baseline.out --format=NeXus \
--dir=$*
# --dir=../results/$1 --ncount=1e9 \
# sample=1 enable_polarizer=1 \
# omegaa=6.0 sample_length=0.048 sample_height=0.01
echo "Program finished with exit code $? at: `date`"
echo "Compressing hdf5 file..."
h5repack -i ../results/$1/mccode.h5 -o ../results/$1/mccode.h5.c -f /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:GZIP=5 -l /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:CHUNK=3072x8
mv ../results/$1/mccode.h5.c ../results/$1/mccode.h5
echo "Scipt finished."
simulation/sbatch_2.sh
-40
View File
@@ -1,40 +0,0 @@
#!/bin/bash
#SBATCH -J mcEstia_2
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout2.log
#SBATCH -e stderr2.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.09 0.01 -0.01)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
simulation/sbatch_3.sh
-40
View File
@@ -1,40 +0,0 @@
#!/bin/bash
#SBATCH -J mcEstia_3
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout3.log
#SBATCH -e stderr3.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq 0.01 0.01 0.09)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
simulation/sbatch_4.sh
-40
View File
@@ -1,40 +0,0 @@
#!/bin/bash
#SBATCH -J mcEstia_4
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout4.log
#SBATCH -e stderr4.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.5 0.1 -0.1)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
simulation/sbatch_5.sh
-40
View File
@@ -1,40 +0,0 @@
#!/bin/bash
#SBATCH -J mcEstia_5
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout5.log
#SBATCH -e stderr5.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq 0.1 0.1 0.5)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
testing.ipynb
+551
View File
File diff suppressed because one or more lines are too long