Compare commits
2 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 glavic_a
|
118b9f5748 | ||
|
glavic_a
|
bf97ff488e |
@@ -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.
|
||||
|
||||
@@ -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'])
|
||||
|
||||
@@ -114,7 +114,7 @@ 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 %{
|
||||
@@ -127,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);
|
||||
@@ -404,7 +404,7 @@ COMPONENT analyzer2 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, r
|
||||
|
||||
/* 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
|
||||
|
||||
@@ -145,8 +145,7 @@ 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
|
||||
@@ -184,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
|
||||
|
||||
@@ -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
|
||||
+160
-95
@@ -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);
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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')
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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 ####################
|
||||
|
||||
|
||||
@@ -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
|
||||
@@ -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 ####################
|
||||
|
||||
|
||||
@@ -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 ####################
|
||||
|
||||
|
||||
@@ -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`"
|
||||
|
||||
@@ -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."
|
||||
@@ -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`"
|
||||
|
||||
@@ -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`"
|
||||
|
||||
@@ -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`"
|
||||
|
||||
@@ -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`"
|
||||
|
||||
+551
File diff suppressed because one or more lines are too long
Reference in New Issue
Block a user