Compare commits
2 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 glavic_a
|
8434fecdde | ||
|
glavic_a
|
9697c9a584 |
@@ -8,7 +8,6 @@ old
|
||||
*.backup
|
||||
.project
|
||||
.settings
|
||||
.idea
|
||||
*.pyc
|
||||
*.pyo
|
||||
hosts
|
||||
|
||||
@@ -1,9 +0,0 @@
|
||||
#!/bin/bash
|
||||
# compress McStas hdf5 files to increase performance and save disk space
|
||||
|
||||
for fi in $*
|
||||
do
|
||||
echo $fi
|
||||
h5repack -i $fi -o $fi.c -f /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:GZIP=5 -l /entry1/data/detector_list_p_x_y_t_L_sx_sy_sz/events:CHUNK=3072x8
|
||||
mv $fi.c $fi
|
||||
done
|
||||
@@ -8,7 +8,7 @@ import sys, os
|
||||
from numpy import *
|
||||
seterr(all='ignore')
|
||||
import estia_help as eh
|
||||
from aglib import mcstas_reader as mr
|
||||
import mcstas_reader as mr
|
||||
|
||||
|
||||
# Scaling factor of source monitor after normalizing by area to get to brilliance.
|
||||
|
||||
@@ -0,0 +1,422 @@
|
||||
#-*- coding: utf-8 -*-
|
||||
'''
|
||||
Simple support library for reading, analyzing and plotting of McStas results
|
||||
from the Estia instrument simulations.
|
||||
Meant to be used from IPython Notebook or QtConsole, but can also be run stand alone.
|
||||
'''
|
||||
|
||||
import os, sys
|
||||
from numpy import *
|
||||
|
||||
try:
|
||||
import h5py
|
||||
except ImportError:
|
||||
print "h5py not found, modern NeXuS format will not be readable."
|
||||
|
||||
try:
|
||||
from IPython import display #@UnusedImport
|
||||
from IPython.core.pylabtools import print_figure
|
||||
from matplotlib.figure import Figure
|
||||
from matplotlib.backends.backend_agg import FigureCanvasAgg
|
||||
from matplotlib.colors import LogNorm
|
||||
except ImportError:
|
||||
# IPython and/or matplotlib not available, no interactive functionalities
|
||||
display=None
|
||||
|
||||
MAX_EVTS_BATCH=50000
|
||||
|
||||
class McSim(object):
|
||||
'''
|
||||
Object representing complete simulation from McStas.
|
||||
Supports either the old style McStas output with separate ASCII files
|
||||
or HDF5 single file output.
|
||||
|
||||
Different monitors can be accessed as keys like in a dictionary. The data is only loaded when the monitor is first accessed.
|
||||
'''
|
||||
|
||||
def __init__(self, path):
|
||||
'''
|
||||
Initialize the simulation. path should be the file name to either the NeXuS or the mccode.sim file.
|
||||
'''
|
||||
self._data={}
|
||||
if path.endswith('.h5'):
|
||||
self._init_hdf(path)
|
||||
elif path.endswith('.sim'):
|
||||
self._init_old(path)
|
||||
else:
|
||||
if os.path.exists(os.path.join(path,'mccode.h5')):
|
||||
self._init_hdf(os.path.join(path,'mccode.h5'))
|
||||
elif os.path.exists(os.path.join(path, 'mccode.sim')):
|
||||
self._init_old(os.path.join(path, 'mccode.sim'))
|
||||
else:
|
||||
raise IOError, "Can't locate mccode.h5 of mccode.sim file in %s"%path
|
||||
|
||||
def _init_hdf(self, path):
|
||||
self.hdf=h5py.File(path, 'r')
|
||||
self.data_loader=DataLoaderHDF(self.hdf)
|
||||
self.info=self.data_loader.info
|
||||
|
||||
def _init_old(self, path):
|
||||
self.info=HeaderFile(path)
|
||||
self.data_loader=DataLoaderOld(self.info, os.path.dirname(path))
|
||||
|
||||
def keys(self):
|
||||
return self.info['data'].keys()
|
||||
|
||||
def monitors(self):
|
||||
output=[]
|
||||
for val in self.info['data'].values():
|
||||
xy=(val['xvar'], val['yvar'])
|
||||
if not xy in output:
|
||||
output.append(xy)
|
||||
output.sort()
|
||||
return output
|
||||
|
||||
def plot(self, monitors=None):
|
||||
graphs={}
|
||||
for key, val in self.info['data'].items():
|
||||
xy=(val['xvar'], val['yvar'])
|
||||
if monitors is not None and xy!=monitors:
|
||||
continue
|
||||
data=self[key]
|
||||
if xy in graphs:
|
||||
graphs[xy].append(data)
|
||||
else:
|
||||
graphs[xy]=[data]
|
||||
for xy, datasets in sorted(graphs.items()):
|
||||
cols=min(len(datasets), 3)
|
||||
rows=len(datasets)/3+1
|
||||
|
||||
fig=Figure(figsize=(12, 5*rows), dpi=300, facecolor='#FFFFFF')
|
||||
FigureCanvasAgg(fig)
|
||||
|
||||
for i, data in enumerate(datasets):
|
||||
ax=fig.add_subplot(rows, cols, i+1)
|
||||
data.plot(ax=ax)
|
||||
graphs[xy]=fig
|
||||
if monitors is None:
|
||||
return graphs
|
||||
else:
|
||||
return fig
|
||||
|
||||
def __getitem__(self, item):
|
||||
if item in self._data:
|
||||
return self._data[item]
|
||||
elif item in self.keys():
|
||||
data=self.data_loader.load_item(item)
|
||||
self._data[item]=data
|
||||
return data
|
||||
else:
|
||||
raise KeyError, "Can't find dataset %s"%item
|
||||
|
||||
class HeaderFile(object):
|
||||
'''
|
||||
Analyze McStas mccode.sim header files.
|
||||
'''
|
||||
_data=None
|
||||
|
||||
@property
|
||||
def data(self):
|
||||
return self._data['data']
|
||||
|
||||
def __init__(self, path):
|
||||
self._data={}
|
||||
data=open(path, 'r').read()
|
||||
data_lines=data.splitlines()
|
||||
if not data.startswith('McStas simulation description file'):
|
||||
raise IOError, 'Not a valid McStas description file.'
|
||||
self._data['start_time']=data_lines[1].split(':', 1)[1].strip()
|
||||
self._data['program_name']=data_lines[2].split(':', 1)[1].strip()
|
||||
self._data['data']=self.get_data(data)
|
||||
|
||||
def get_data(self, data):
|
||||
output={}
|
||||
start_idx=0
|
||||
while start_idx<len(data):
|
||||
start_idx=data.find('begin data', start_idx)
|
||||
if start_idx==-1:
|
||||
break
|
||||
end_idx=data.find('end data', start_idx)
|
||||
block=data[start_idx:end_idx].splitlines()[1:]
|
||||
block_info={}
|
||||
for bline in block:
|
||||
item, value=bline.strip().split(':', 1)
|
||||
block_info[item]=value.strip()
|
||||
if item=='component':
|
||||
output[value.strip()]=block_info
|
||||
|
||||
start_idx=end_idx
|
||||
return output
|
||||
|
||||
def __getitem__(self, item):
|
||||
return self._data[item]
|
||||
|
||||
class DataLoaderOld(object):
|
||||
'''
|
||||
Load and analyze a old style McStas format with a simulation and a set of data text files.
|
||||
'''
|
||||
def __init__(self, info, root):
|
||||
self.info=info
|
||||
self.root=root
|
||||
|
||||
def load_item(self, item):
|
||||
item_info=self.info['data'][item]
|
||||
if not item_info['type'].startswith('array_2d'):
|
||||
return self.load_item_1d(item)
|
||||
|
||||
fname=os.path.join(self.root, item_info['filename'])
|
||||
x_col=item_info['xvar']
|
||||
y_col=item_info['yvar']
|
||||
if x_col=='Li' and y_col=='p': # Detector_nD
|
||||
cols=item_info['variables'].split()
|
||||
data=loadtxt(fname, dtype={'names': cols, 'formats': ['f4']*len(cols)})
|
||||
return TofData(data, item_info)
|
||||
else:
|
||||
raw=loadtxt(fname)
|
||||
data=raw[:len(raw)/3]
|
||||
return Dataset(data, item_info)
|
||||
|
||||
def load_item_1d(self, item):
|
||||
item_info=self.info['data'][item]
|
||||
fname=os.path.join(self.root, item_info['filename'])
|
||||
raw=loadtxt(fname).T
|
||||
data=raw[1]
|
||||
errors=raw[2]
|
||||
return Dataset1D(data, errors, item_info)
|
||||
|
||||
class DataLoaderHDF(object):
|
||||
'''
|
||||
Load and analyze a new style NeXuS file format.
|
||||
'''
|
||||
def __init__(self, hdf):
|
||||
self.hdf=hdf['entry1']
|
||||
self.info={}
|
||||
self.info['data']={}
|
||||
for item in self.hdf['data'].keys():
|
||||
node=self.hdf['data/'+item]
|
||||
info={}
|
||||
for key, value in node.attrs.items():
|
||||
info[key.strip()]=value.strip()
|
||||
if info['filename'][-4:]=='.dat' or '_list.' in info['filename']:
|
||||
self.info['data'][info['component']]=info
|
||||
else:
|
||||
self.info['data'][info['filename']]=info
|
||||
info['datapath']='data/'+item
|
||||
|
||||
def load_item(self, item):
|
||||
item_info=self.info['data'][item]
|
||||
if not item_info['type'].startswith('array_2d'):
|
||||
return self.load_item_1d(item)
|
||||
|
||||
node=self.hdf[item_info['datapath']]
|
||||
x_col=item_info['xvar']
|
||||
y_col=item_info['yvar']
|
||||
if x_col=='Li' and y_col=='p': # Detector_nD
|
||||
cols=item_info['variables'].split()
|
||||
evds=node['events']
|
||||
if len(evds)<=MAX_EVTS_BATCH:
|
||||
data=evds.value.astype(float32).view(
|
||||
dtype={'names': cols, 'formats': ['f4']*len(cols)}).flatten()
|
||||
else:
|
||||
ds=[]
|
||||
sys.stdout.write('Reading large dataset:\n')
|
||||
for i in range(len(evds)//MAX_EVTS_BATCH+1):
|
||||
sys.stdout.write('\r%i/%i'%(i*MAX_EVTS_BATCH, len(evds)))
|
||||
sys.stdout.flush()
|
||||
ds.append(evds[i*MAX_EVTS_BATCH:(i+1)*MAX_EVTS_BATCH])
|
||||
sys.stdout.write('\r%i/%i\n'%(len(evds), len(evds)))
|
||||
data=vstack(ds).astype(float32).view(
|
||||
dtype={'names': cols, 'formats': ['f4']*len(cols)}).flatten()
|
||||
return TofData(data, item_info)
|
||||
else:
|
||||
data=node['data'].value.T
|
||||
return Dataset(data, item_info)
|
||||
|
||||
def load_item_1d(self, item):
|
||||
item_info=self.info['data'][item]
|
||||
node=self.hdf[item_info['datapath']]
|
||||
data=node['data'].value
|
||||
errors=node['errors'].value
|
||||
return Dataset1D(data, errors, item_info)
|
||||
|
||||
class Dataset1D(object):
|
||||
'''
|
||||
Representation of a standard McStas dataset of one variable.
|
||||
'''
|
||||
def __init__(self, data, errors, info):
|
||||
self.data=data
|
||||
self.errors=errors
|
||||
self.info=info
|
||||
|
||||
def plot(self, log=False, ax=None):
|
||||
if ax is None:
|
||||
import pylab
|
||||
ax=pylab.gca()
|
||||
|
||||
limits=map(float, self.info['xlimits'].split())
|
||||
x=linspace(limits[0], limits[1], len(self.data))
|
||||
|
||||
ax.errorbar(x, self.data, yerr=self.errors)
|
||||
if log:
|
||||
ax.set_yscale('log')
|
||||
else:
|
||||
ax.set_yscale('linear')
|
||||
ax.set_xlabel(self.info['xlabel'])
|
||||
ax.set_ylabel(self.info['ylabel'])
|
||||
ax.set_title(self.info['component'])
|
||||
|
||||
def _repr_png_(self):
|
||||
'''
|
||||
Image representation form ipython console/notebook.
|
||||
'''
|
||||
fig=Figure(figsize=(8, 5), dpi=300, facecolor='#FFFFFF')
|
||||
FigureCanvasAgg(fig)
|
||||
ax=fig.add_subplot(111)
|
||||
|
||||
self.plot(ax=ax)
|
||||
|
||||
return print_figure(fig, dpi=72)
|
||||
|
||||
|
||||
class Dataset(object):
|
||||
'''
|
||||
Representation of a standard McStas dataset with 2D view.
|
||||
'''
|
||||
def __init__(self, data, info):
|
||||
self.data=data
|
||||
self.info=info
|
||||
|
||||
def plot(self, log=False, ax=None):
|
||||
if ax is None:
|
||||
import pylab
|
||||
ax=pylab.gca()
|
||||
|
||||
limits=map(float, self.info['xylimits'].split())
|
||||
|
||||
if log:
|
||||
ax.imshow(self.data, 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, basestring):
|
||||
fltr=eval(fltr, globals(), columns)
|
||||
I, x=histogram(columns[col][fltr], bins=bins, weights=w[fltr])
|
||||
return x, I
|
||||
|
||||
def project2d(self, xcol, ycol, bins=50, fltr=None, newcols=None):
|
||||
'''
|
||||
Generate 2D binned data for arbitrary binning and columns.
|
||||
|
||||
The fltr argument can be a string with a filtering condition on the available columns.
|
||||
Syntax for this filtering follows numpy array convetions like (x>5)&(abs(L)<0.3) would
|
||||
be a valid statement.
|
||||
|
||||
newcols can be a list of (name, code) tuples, that calculate new columns from
|
||||
existing ones, like newcols=[('q', '4*pi/L*sin(theta)')].
|
||||
'''
|
||||
columns=dict([(coli, self.data[coli]) for coli in self.data.dtype.names])
|
||||
if newcols is not None:
|
||||
for name, code in newcols:
|
||||
columns[name]=eval(code, globals(), columns)
|
||||
|
||||
if fltr is None:
|
||||
I, y, x=histogram2d(columns[ycol], columns[xcol],
|
||||
bins=bins, weights=self.data['p'])
|
||||
else:
|
||||
if isinstance(fltr, basestring):
|
||||
fltr=eval(fltr, globals(), columns)
|
||||
I, y, x=histogram2d(columns[ycol][fltr], columns[xcol][fltr],
|
||||
bins=bins, weights=columns['p'][fltr])
|
||||
return x, y, I
|
||||
|
||||
def plot(self, xcol='x', ycol='y', log=False, ax=None, bins=50, fltr=None, newcols=None):
|
||||
if ax is None:
|
||||
import pylab
|
||||
ax=pylab.gca()
|
||||
|
||||
x, y, I=self.project2d(xcol, ycol, bins=bins, fltr=fltr, newcols=newcols)
|
||||
|
||||
if log:
|
||||
ax.pcolormesh(x, y, I, norm=LogNorm())
|
||||
else:
|
||||
ax.pcolormesh(x, y, I)
|
||||
|
||||
col_names=self.info['title'].split()
|
||||
cols=self.info['variables'].split()
|
||||
|
||||
try:
|
||||
ax.set_xlabel(xcol+'-'+col_names[cols.index(xcol)])
|
||||
except ValueError:
|
||||
ax.set_xlabel(xcol)
|
||||
try:
|
||||
ax.set_ylabel(ycol+'-'+col_names[cols.index(ycol)])
|
||||
except ValueError:
|
||||
ax.set_ylabel(ycol)
|
||||
ax.set_title(self.info['component'])
|
||||
|
||||
|
||||
def plot1d(self, col='x', log=False, ax=None, bins=50, fltr=None, newcols=None):
|
||||
if ax is None:
|
||||
import pylab
|
||||
ax=pylab.gca()
|
||||
|
||||
x, I=self.project1d(col, bins=bins, fltr=fltr, newcols=newcols)
|
||||
|
||||
if log:
|
||||
ax.semilogy((x[:-1]+x[1:])/2., I)
|
||||
else:
|
||||
ax.plot((x[:-1]+x[1:])/2., I)
|
||||
col_names=self.info['title'].split()
|
||||
cols=self.info['variables'].split()
|
||||
|
||||
try:
|
||||
ax.set_xlabel(col+'-'+col_names[cols.index(col)])
|
||||
except ValueError:
|
||||
ax.set_xlabel(col)
|
||||
ax.set_ylabel('Intensity')
|
||||
ax.set_title(self.info['component'])
|
||||
|
||||
@@ -0,0 +1,753 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Dose_calculator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Rodion Kolevatov
|
||||
* Date: August 2018
|
||||
* Version: $Revision: 1.0 $
|
||||
* Release: ---
|
||||
* Origin: IFE
|
||||
*
|
||||
* Calculating dose rate at the outer surface of lateral shielding of uniform thickness
|
||||
* made of material specified in the component input (concrete, iron or lead).
|
||||
* Shielding starts at a distance 0.5*Innerspace from the center of the guide.
|
||||
*
|
||||
* %D
|
||||
* A Shielding_logger together with a set of Shielding_iterators is required.
|
||||
* Iterators should contain Monitor_nD components for recording capture rates.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
* Shielding type: Concrete, Iron or Lead
|
||||
* Names of text files where the coating capture rates are recorded
|
||||
* Inner space in the shielding given by parameter Innerspace (in meters).
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Dose_calculator
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS (string Material = "Fe", double Innerspace=0.3, double Thickness=0.2, string SteelTubing="", double TubingThickness=0.0, string NiCaptureFile="NiCapture.dat",
|
||||
string TiCaptureFile = "TiCapture.dat", string TotalCaptureFile = "TotalCapture.dat", string OutputFile="Dose.dat")
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef LIN_INT_ROUTINE
|
||||
#define LIN_INT_ROUTINE 1
|
||||
//Multi-d linear interpolation routine. Array of args, number of args, sizes of args in interpolation grid, argument grid in single line, data in single line
|
||||
double lint (double * args, int Narg, int * sizearg, double * argtable, double * datatable)
|
||||
{
|
||||
if (Narg==1) {
|
||||
// printf ("Narg=1, arg = %g\n", *args);
|
||||
int point=0;
|
||||
if ((*args)<(*argtable) ) return *datatable;
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return * (datatable+ *sizearg-1); // if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) point++;
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup = (*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*(*(datatable+point))+weightlow*(*(datatable+point-1));
|
||||
}//arg within range of values
|
||||
}// if Narg==1
|
||||
else if (Narg >1){//if more than one argument
|
||||
//lookup how large is the data (Narg-1)D slice for fixed value of the first argument
|
||||
int slicesize=1;
|
||||
int i,point=0;
|
||||
for ( i=1;i<Narg;i++) slicesize*=sizearg[i];
|
||||
// printf("Narg > 1, slicesize = %d, arggridstart = %g, argument = %g, arggridend = %g\n",slicesize,*argtable,*args, *(argtable+ *sizearg-1) );
|
||||
//search weights for the first argument and get results with one argument less
|
||||
if ((*args)<(*argtable) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable);
|
||||
// if the value is too low -- return what is on the lower bound
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), (datatable+ slicesize*(*sizearg)));
|
||||
// if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) {
|
||||
// printf ("*(argtable+point) = %g\n", *(argtable+point));
|
||||
point++;}
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup =(*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+point*slicesize )
|
||||
+weightlow*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+(point-1)*slicesize);
|
||||
} // argument within the ragne
|
||||
} //Narg>1
|
||||
};
|
||||
#endif
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
#include <dirent.h>
|
||||
#include <unistd.h>
|
||||
|
||||
//*** NEUTRON CONVERSION DATA *****//
|
||||
#ifndef NEUTRON_GAMMA_DATA
|
||||
#define NEUTRON_GAMMA_DATA 1
|
||||
//Photons per capture in NiTi coating
|
||||
double fraction[]={0.003565789, 0.009694611, 0.028838269, 0.149117282, 0.117348519, 0.042269932, 0.386434136, 0.092896131, 0.046767054, 0.032512648, 0.043547581, 0.040537585, 0.510133557, 0.320364465};
|
||||
double energy [] = {0.150,
|
||||
0.200,
|
||||
0.300,
|
||||
0.400,
|
||||
0.500,
|
||||
0.600,
|
||||
0.800,
|
||||
1.,
|
||||
1.50,
|
||||
2.000,
|
||||
3.0,
|
||||
4.0,
|
||||
5.000,
|
||||
6.000,
|
||||
7.0,
|
||||
8.000,
|
||||
9.0,
|
||||
10.000,
|
||||
11.000};
|
||||
int nEgroup= 19;
|
||||
//Spectrum of capture photons from Ni and Ti
|
||||
//Ni gamma capture specrum
|
||||
double fractionNi[]={0.003722,//150
|
||||
0.015513,//200
|
||||
0.048432,//300
|
||||
0.047485,//400
|
||||
0.198059,//500
|
||||
0.002809,//600
|
||||
0.003142,//800
|
||||
0.06579,//1000
|
||||
0.054534,//1500
|
||||
0.028872,//2000
|
||||
0.063882,//3000
|
||||
0.050286,//4000
|
||||
0.035567,//5000
|
||||
0.080326,//6000
|
||||
0.156571,//7000
|
||||
0.130355,//8000
|
||||
0.534153,//9000
|
||||
0.000558,//10000
|
||||
0.000260};//11000
|
||||
|
||||
|
||||
//Ti capture gamma spectrum
|
||||
double fractionTi[]={0.009326, //0.15
|
||||
0.0, //0.2
|
||||
0.001938,//0.3
|
||||
0.302632,//0.4
|
||||
0.000882,//0.5
|
||||
0.000668, //0.6
|
||||
0.000914, //0.8
|
||||
0.023046,//1.0
|
||||
0.947104,//1.5
|
||||
0.232257,//2.0
|
||||
0.086418,//3.0
|
||||
0.162878,//4.0
|
||||
0.11209,//5.0
|
||||
0.013419,//6.0
|
||||
0.0875051, //7.0
|
||||
0.009083,//8.0
|
||||
0.001923, //9.0
|
||||
0.002189,//10.0
|
||||
0.000332};//11.0
|
||||
//double fractionB[]={0.0, 0.0, 0.0, 0.0, 0.93, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0};
|
||||
double fractionB[]={0,//0.15
|
||||
0,//0.2
|
||||
0,//0.3
|
||||
2.32983E-06,//0.4
|
||||
0.935769511,//0.5
|
||||
1.94926E-06,//0.6
|
||||
2.22356E-05,//0.8
|
||||
3.29279E-05,//1.0
|
||||
0.000248357,//1.5
|
||||
0.000398671,//2.0
|
||||
0.000410038,//3.0
|
||||
0.001042575,//4.0
|
||||
0.000961964,//5.0
|
||||
8.714E-05,//6.0
|
||||
0.000248803,//7.0
|
||||
0.00017111,//8.0
|
||||
3.9074E-05,//9.0
|
||||
7.28E-07,//10.0
|
||||
6.27e-6};//11.0
|
||||
|
||||
/*
|
||||
//conversion test input
|
||||
double fraction [] = {0.5,0.5};
|
||||
double energy [] = {2.0,2.0};
|
||||
int nEgroup= 2;
|
||||
*/
|
||||
//*** SHIELDING DATA TABLES ****//
|
||||
// Mashkovich set + NIST data for 15 MeV
|
||||
// Linear attenuation for concrete
|
||||
double AttenuationArgsConc[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataConc[]={31.7, 28.5, 24.6, 21.9, 20.0, 18.5, 16.3, 14.6, 13.1, 11.9, 10.3, 8.74, 8.37, 7.34, 6.65, 6.19, 5.61, 5.29, 4.8208}; // in m^-1
|
||||
int AttenuationSizeConc[]={18};
|
||||
// Iron
|
||||
double AttenuationArgsFe[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataFe[]={139., 106., 83.3, 71.7, 64.6, 59.5, 52.0, 46.7, 42.2, 38.1, 33.3, 29.1, 28.4, 26.0, 24.8, 24.0, 23.4, 23.4, 24.3}; //in m^-1
|
||||
int AttenuationSizeFe[]={19};
|
||||
double AttenuationArgsPb[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataPb[]={2180., 1070., 425., 244., 170., 133., 95.2, 77.1, 65.8, 57.7, 50.8, 47.6, 46.8, 47.2, 48.1, 49.4, 52.0, 55.0, 64.2}; //in m^-1
|
||||
int AttenuationSizePb[]={19};
|
||||
/*
|
||||
//NIST set
|
||||
//Concrete
|
||||
double AttenuationArgsConc[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataConc[]={
|
||||
33.028,
|
||||
29.486,
|
||||
25.231,
|
||||
22.5009,
|
||||
20.5045,
|
||||
18.9428,
|
||||
16.6221,
|
||||
14.9385,
|
||||
13.3561,
|
||||
12.1624,
|
||||
10.4811,
|
||||
8.5123,
|
||||
7.3991,
|
||||
6.6884,
|
||||
6.2031,
|
||||
5.5936,
|
||||
5.2394,
|
||||
4.8208};
|
||||
int AttenuationSizeConc[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsFe[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataFe[]={
|
||||
1.55E+02,
|
||||
1.15E+02,
|
||||
8.65E+01,
|
||||
7.40E+01,
|
||||
6.63E+01,
|
||||
6.07E+01,
|
||||
5.27E+01,
|
||||
4.72E+01,
|
||||
4.21E+01,
|
||||
3.84E+01,
|
||||
3.36E+01,
|
||||
2.85E+01,
|
||||
2.61E+01,
|
||||
2.48E+01,
|
||||
2.41E+01,
|
||||
2.36E+01,
|
||||
2.36E+01,
|
||||
2.43E+01};
|
||||
int AttenuationSizeFe[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsPb[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in meV
|
||||
double AttenuationDataPb[]={
|
||||
2.28E+03,
|
||||
1.13E+03,
|
||||
4.57E+02,
|
||||
2.63E+02,
|
||||
1.83E+02,
|
||||
1.42E+02,
|
||||
1.01E+02,
|
||||
8.05E+01,
|
||||
6.66E+01,
|
||||
5.92E+01,
|
||||
5.22E+01,
|
||||
4.80E+01,
|
||||
4.76E+01,
|
||||
4.84E+01,
|
||||
4.98E+01,
|
||||
5.30E+01,
|
||||
5.64E+01,
|
||||
6.42E+01
|
||||
}; //in m^-1
|
||||
int AttenuationSizePb[]={18};
|
||||
*/
|
||||
|
||||
/*
|
||||
//Linear attenuation test input
|
||||
double AttenuationArgs[]={0.15, 10.0}; // in meV
|
||||
double AttenuationData[]={7.0, 7.0 }; // in m^-1
|
||||
int AttenuationSize[]={2};
|
||||
*/
|
||||
|
||||
// Dose buildup factors concrete
|
||||
double BuildupDataConc[]={1., 1.74, 2.26, 2.95, 3.79, 4.51, 5.57, 6.51, 3.18,
|
||||
1., 2.82, 5.13, 11.2, 24.2, 42.7, 87.6, 153., 353.,
|
||||
1., 2.52, 4.66, 10.8, 25.6, 48.2, 107., 198., 497.,
|
||||
1., 2.27, 4.03, 8.97, 20.2, 30.4, 75.6, 131, 292,
|
||||
1., 1.98, 3.24, 6.42, 12.7, 20.7, 37.2, 57.1, 106,
|
||||
1., 1.77, 2.65, 4.61, 7.97, 11.7, 18.6, 26.0, 42.2,
|
||||
1., 1.67, 2.38, 3.84, 6.20, 8.71, 13.1, 17.7, 27.4,
|
||||
1., 1.61, 2.18, 3.37, 5.23, 7.15, 10.5, 13.9, 20.9,
|
||||
1., 1.49, 1.93, 2.80, 4.14, 5.52, 7.86, 10.2, 15.5,
|
||||
1., 1.41, 1.76, 2.45, 3.51, 4.59, 6.43, 8.31, 12.2,
|
||||
1., 1.35, 1.64, 2.22, 3.10, 4.01, 5.57, 7.19, 10.6,
|
||||
1., 1.26, 1.46, 1.86, 2.50, 3.16, 4.34, 5.59, 8.27};
|
||||
|
||||
int BuildupSizeConc[] ={12, //energy
|
||||
9}; // mu d
|
||||
double BuildupArgsConc[]={
|
||||
0.05, 0.15, 0.3, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 // energy in MeV
|
||||
,
|
||||
0., 1., 2., 4., 7., 10., 15., 20., 30. //mu d
|
||||
};
|
||||
//Dose buildup factors Fe
|
||||
double BuildupDataFe[]={
|
||||
1., 1.5, 2.2, 3.1, 4.1, 4.6, 5.4, 5.9, //0.1
|
||||
1., 2.0, 3.1, 5.3, 8.9, 14., 22., 31., //0.2
|
||||
1., 2.1, 3.3, 6.0, 12., 23., 49., 84., //0.4
|
||||
1., 1.98, 3.09, 5.98, 11.7, 19.2, 35.4, 55.6, //0.5
|
||||
1., 1.87, 2.89, 5.39, 10.2, 16.2, 28.3, 42.7, //1.0
|
||||
1., 1.76, 2.43, 4.13, 7.25, 10.9, 17.6, 25.1, //2.0
|
||||
1., 1.55, 2.15, 3.51, 5.85, 8.51, 13.5, 19.1, //3.0
|
||||
1., 1.45, 1.94, 3.03, 4.91, 7.11, 11.2, 16.0, //4.0
|
||||
1., 1.34, 1.72, 2.58, 4.14, 6.02, 9.89, 14.7, //6.0
|
||||
1., 1.27, 1.56, 2.23, 3.49, 5.07, 8.50, 13.0, //8.0
|
||||
1., 1.20, 1.42, 1.95, 2.99, 4.35, 7.54, 12.4, //9.0
|
||||
1., 1.48, 1.86, 2.72, 4.30 , 6.37, 11.4, 19.1};//15.0
|
||||
int BuildupSizeFe[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsFe[]={0.1, 0.2, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
//Dose buildup factors Pb
|
||||
double BuildupDataPb[]={1., 1.01, 1.03, 1.06, 1.15, 1.16, 1.18, 1.19,
|
||||
1., 1.11, 1.17, 1.25, 1.34, 1.41, 1.5, 1.56,
|
||||
1., 1.17, 1.29, 1.46, 1.58, 1.72, 1.89, 2.02,
|
||||
1., 1.24, 1.42, 1.69, 2.00, 2.27, 2.65, 2.73,
|
||||
1., 1.37, 1.69, 2.26, 3.02, 3.74, 4.81, 5.86,
|
||||
1., 1.39, 1.76, 2.51, 3.66, 4.84, 6.87, 9.00,
|
||||
1., 1.34, 1.68, 2.43, 3.75, 5.30, 8.44, 12.3,
|
||||
1., 1.27, 1.56, 2.25, 3.61, 5.44, 9.80, 16.3,
|
||||
1., 1.21, 1.46, 2.08, 3.44, 5.55, 11.7, 23.6,
|
||||
1., 1.18, 1.40, 1.97, 3.34, 5.69, 13.8, 32.7,
|
||||
1., 1.14, 1.30, 1.74, 2.89, 5.07, 14.1, 44.6,
|
||||
1., 1.11, 1.23, 1.58, 2.52, 4.34, 12.5, 39.2};
|
||||
|
||||
int BuildupSizePb[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsPb[]={0.15, 0.30, 0.40, 0.5, 1.0, 2.0, 3.0, 4.0, 5.1, 6.0, 8.0, 10.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
|
||||
|
||||
// Flux-to-dose conversion
|
||||
/* //NRB-99?
|
||||
double FtoDArgs[]= {0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.80, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*0.752E-10, 3600*1.00E-10, 3600*1.51E-10, 3600*2.00E-10, 3600*2.47E-10, 3600*2.91E-10, 3600*3.73E-10, 3600*4.48E-10, 3600*7.49E-10, 3600*12.0E-10, 3600*16.0E-10, 3600*19.9E-10, 3600*23.8E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 13;
|
||||
int FtoDSize[] = {13};
|
||||
*/
|
||||
//ESS
|
||||
double FtoDArgs[]= {
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.511,
|
||||
0.6,
|
||||
0.662,
|
||||
0.8,
|
||||
1.,
|
||||
1.12,
|
||||
1.33,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
6.13,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in MeV
|
||||
double FtoDData[] = {
|
||||
2.69E-07,
|
||||
3.60E-07,
|
||||
5.44E-07,
|
||||
7.20E-07,
|
||||
8.89E-07,
|
||||
9.07E-07,
|
||||
1.05E-06,
|
||||
1.14E-06,
|
||||
1.34E-06,
|
||||
1.62E-06,
|
||||
1.76E-06,
|
||||
2.01E-06,
|
||||
2.20E-06,
|
||||
2.69E-06,
|
||||
3.51E-06,
|
||||
4.21E-06,
|
||||
4.82E-06,
|
||||
5.40E-06,
|
||||
5.47E-06,
|
||||
6.70E-06,
|
||||
7.92E-06,
|
||||
1.09E-05
|
||||
}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 22;
|
||||
int FtoDSize[] = {22};
|
||||
|
||||
|
||||
|
||||
// test input flux to dose
|
||||
/*
|
||||
double FtoDArgs[]= {0.15, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*7.0E-10,3600*7.0E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 2;
|
||||
int FtoDSize[] = {2};
|
||||
*/
|
||||
#endif //end of insertion of datatables.
|
||||
/*** end of insertion to DECLARE section **/
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
|
||||
/** Define material for dose calculation **/
|
||||
|
||||
int imat=0, itubing=0;
|
||||
if(strcmp(SteelTubing,"yes")==0||strcmp(SteelTubing,"Yes")==0||strcmp(SteelTubing,"YES")==0||strcmp(SteelTubing,"y")==0||strcmp(SteelTubing,"Y")==0||strcmp(SteelTubing,"1")==0){
|
||||
itubing=1;
|
||||
imat=0;
|
||||
printf("Calculating dose for STEEL TUBING\n Outer material will be CONCRETE\n");}
|
||||
else if((strcmp(Material,"Fe")==0)||(strcmp(Material,"fe")==0)||(strcmp(Material,"Iron")==0)||(strcmp(Material,"iron")==0)||(strcmp(Material,"IRON")==0)){
|
||||
itubing=0;
|
||||
imat=1;
|
||||
printf("Calculating dose for IRON\n");}
|
||||
else if ((strcmp(Material,"Pb")==0)||(strcmp(Material,"pb")==0)||(strcmp(Material,"Lead")==0)||(strcmp(Material,"lead")==0)||(strcmp(Material,"LEAD")==0))
|
||||
{itubing=0;
|
||||
imat=2;
|
||||
printf("Calculating dose for LEAD\n");}
|
||||
else if((strcmp(Material,"Conc")==0)||(strcmp(Material,"Concrete")==0)||(strcmp(Material,"conc")==0)||(strcmp(Material,"concrete")==0))
|
||||
{itubing=0;
|
||||
imat=0;
|
||||
printf("Calculating dose for CONCRETE\n");}
|
||||
else
|
||||
{
|
||||
printf ("No/wrong material specified for dose calculation, possible options: Fe, Pb, concrete\n");
|
||||
printf ("Using default: concrete\n");
|
||||
}
|
||||
/**Reading datafiles with capture per bin**/
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
double d1, d2, d3, d4;
|
||||
int ibin=0,NBINS;
|
||||
char line[1000],line1[1000];
|
||||
char dirname[1000];
|
||||
memset(dirname,0,sizeof(dirname));
|
||||
strcat(strcat(dirname,mcdirname),separator);
|
||||
|
||||
char filename[1000];
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FILE *dataIn;
|
||||
|
||||
memset(filename,0,sizeof(filename));//resetting filename to NULL string
|
||||
strcat(strcat(filename,dirname),NiCaptureFile);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",NiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4) ibin++;
|
||||
};
|
||||
|
||||
|
||||
NBINS=ibin;
|
||||
|
||||
//printf("Shielding calculator:\n Input file contains %d bins.\n", NBINS);
|
||||
|
||||
double* zhat = malloc(NBINS*sizeof(double));
|
||||
double* Ihat= malloc(NBINS*sizeof(double));
|
||||
double* IhatNi= malloc(NBINS*sizeof(double));
|
||||
double* IhatTi= malloc(NBINS*sizeof(double));
|
||||
rewind(dataIn);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatNi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
memset(filename,0,sizeof(filename));// setting filename to NULL string
|
||||
strcat(strcat(filename,dirname),TiCaptureFile); //setting filename to TiCapture file
|
||||
//printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatTi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
fclose(dataIn);
|
||||
|
||||
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),TotalCaptureFile);
|
||||
//printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TotalCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
Ihat[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
int i,iz,iEgroup; // Auxiliary variables
|
||||
|
||||
double zbinlength = zhat[2]-zhat[1];
|
||||
|
||||
|
||||
|
||||
double RRR = (0.5*Innerspace)+Thickness+itubing*TubingThickness;
|
||||
double* doseUnshielded= malloc(NBINS*sizeof(double));
|
||||
double* doseShielded= malloc(NBINS*sizeof(double));
|
||||
|
||||
|
||||
|
||||
|
||||
/*Moving along the bins, calculating the dose */
|
||||
|
||||
|
||||
|
||||
for ( iz=0; iz< NBINS; iz++)
|
||||
{double zc=zhat[iz]; // the z-value for the calculation
|
||||
double z;
|
||||
doseShielded[iz]=0.0;
|
||||
doseUnshielded[iz]=0.0;
|
||||
|
||||
/*computing the integral*/
|
||||
double dmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
doseUnshielded[iz]+=(fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926/(RRR*RRR+z*z)*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData);
|
||||
if (imat==2){//LEAD
|
||||
dmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*Thickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);}
|
||||
else if (imat==1){//IRON
|
||||
dmfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*Thickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);}
|
||||
else {//CONCRETE, imat=0, default
|
||||
dmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*Thickness*sqrt(RRR*RRR+z*z)/RRR + itubing*lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TubingThickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);}
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
} /* Position scan iz */
|
||||
|
||||
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
FILE *dose_out;
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),OutputFile);
|
||||
printf("Dose calculator:\nWriting to file %s\n",filename);
|
||||
dose_out=fopen(filename,"w");
|
||||
|
||||
|
||||
fprintf(dose_out,"#Dose at the outer surface of lateral %s shielding of thickness %g m\n#Guide is in the center of the innner space of width %g m\n#Presense of Steel Tubing = %s; of %g m\n#Position z,m\tOverall loss, n/s/m\tNi capture, n/s/m,\tTi capture, n/s/m\tDose unshielded, uSv/h\tDose shielded, uSv/h\n",Material, Thickness, Innerspace, SteelTubing,TubingThickness);
|
||||
|
||||
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
fprintf( dose_out,"%g\t%.3g\t%.3g\t%.3g\t%.3g\t%.3g\n",zhat[i], Ihat[i]/zbinlength, IhatNi[i]/zbinlength, IhatTi[i]/zbinlength,doseUnshielded[i],doseShielded[i]);
|
||||
}
|
||||
fclose(dose_out);
|
||||
free(zhat); free(Ihat); free(IhatNi); free(IhatTi);
|
||||
free(doseShielded); free(doseUnshielded);
|
||||
|
||||
printf("Done with writing dose rates\n");
|
||||
%}
|
||||
|
||||
END
|
||||
File diff suppressed because it is too large
Load Diff
+189
-74
@@ -49,13 +49,12 @@
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE INSTRUMENT ESS_reflectometer_Estia
|
||||
(double omegaa = 1.2, int sample = 0, double sample_length = 0.01, double sample_height = 0.01,
|
||||
(double omegaa = 25, int sample = 0, double sample_length = 0.02, double sample_height = 0.015,
|
||||
int operationmode = 0, double over_illumination = 0.0, double theta_resolution = 0.04,
|
||||
double lambda_min = 3.75, double lambda_start = 3.0, double lambda_end = 12.0,
|
||||
double lambda_min = 3.75, double lambda_start = 0.1, double lambda_end = 12.0,
|
||||
int enable_chopper = 0, int enable_gravity=0, int enable_windows=1,
|
||||
int enable_polarizer = 0, int enable_analyzer = 0,
|
||||
double pol1_start=0.3, double pol1_angle=1.66, double pol2_start=0.3, double pol2_angle=1.66,
|
||||
double source_power = 2,
|
||||
int enable_polarizer = 0, int enable_analyzer = 0,
|
||||
double source_power = 5,
|
||||
double selene1_foot1y =0.0, double selene1_foot2y = 0.0
|
||||
)
|
||||
|
||||
@@ -105,20 +104,12 @@ 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
|
||||
int p_int;
|
||||
%}
|
||||
|
||||
USERVARS %{
|
||||
int p_int; // a flag that gets incremented if a polarizer mirror scatters
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
@@ -127,11 +118,11 @@ 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);
|
||||
//printf(" Chopper open = %.1f deg\n", chopper_open);
|
||||
printf(" Chopper open = %.1f deg\n", chopper_open);
|
||||
|
||||
|
||||
dist_ana_vfocus=analyzer_max_height/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
|
||||
@@ -144,8 +135,8 @@ Theta2_analyzer2=atan((dist_ana_vfocus+analyzer2_start+analyzer2_length)/dist_an
|
||||
selene1_center=(selene1_foot1+selene1_foot2)/2.;
|
||||
selene1_shift=(selene1_foot1y+selene1_foot2y)/1000.0;
|
||||
selene1_rot=atan((selene1_foot2y-selene1_foot1y)/(selene1_foot2-selene1_foot1)/1000.0)*180.0/PI;
|
||||
//printf(" Selene 1 rotation = %.4f deg\n", selene1_rot);
|
||||
//printf(" Selene 1 shift = %.1f mm\n", selene1_shift*1e3);
|
||||
printf(" Selene 1 rotation = %.4f deg\n", selene1_rot);
|
||||
printf(" Selene 1 shift = %.1f mm\n", selene1_shift*1e3);
|
||||
|
||||
%}
|
||||
TRACE
|
||||
@@ -211,7 +202,7 @@ COMPONENT moderator = ESS_butterfly(
|
||||
sector = "E", beamline = 2, yheight = 0.03, cold_frac = 0.9,
|
||||
focus_xw = E02_01_01_Cu_in_xmax-E02_01_01_Cu_in_xmin+0.002,
|
||||
focus_yh = E02_01_01_Cu_in_yheight+0.002,
|
||||
target_index=3, Lmin = lambda_start, Lmax = lambda_end,
|
||||
target_index=4, Lmin = lambda_start, Lmax = lambda_end,
|
||||
n_pulses = 1+enable_chopper, acc_power=source_power)
|
||||
AT (0, 0, 0) RELATIVE origin
|
||||
|
||||
@@ -243,15 +234,6 @@ 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(nb=21, ndiv=41, filename="vs_hordiv",
|
||||
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0,
|
||||
restore_neutron=1)
|
||||
AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT vs_divergence_v = DivPos_monitor(nb=21, ndiv=41, filename="vs_verdiv",
|
||||
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0,
|
||||
restore_neutron=1, vertical=1)
|
||||
AT (0, 0, -0.5*sample_length-0.001) RELATIVE arm_selene1
|
||||
|
||||
/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
|
||||
COMPONENT virtual_source_TR = Slit(
|
||||
@@ -271,6 +253,20 @@ COMPONENT virtual_source_BL = Slit(
|
||||
AT (over_illumination, 0, 0.5*sample_length) RELATIVE arm_virtual_source
|
||||
|
||||
|
||||
/******* Start of scatter logger. Logging scatterings and absorption in the guide system. *********/
|
||||
COMPONENT log_P_start=Shielding_logger()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
#define scatter_logger_stop log_P_start
|
||||
%}
|
||||
|
||||
|
||||
|
||||
|
||||
/*************************************
|
||||
* Geometry of Selene guide separate *
|
||||
*************************************/
|
||||
@@ -280,6 +276,156 @@ COMPONENT virtual_source_BL = Slit(
|
||||
|
||||
%include "Estia_selene2.instr"
|
||||
|
||||
|
||||
|
||||
COMPONENT Lambda_EndGuide = L_monitor(
|
||||
nL = 1000, filename = "Lambda_EndGuide", Lmin=0.5, Lmax=100,
|
||||
restore_neutron = 1)
|
||||
AT (0, 0, -slit_distance-0.01) RELATIVE arm_sample_beam
|
||||
|
||||
|
||||
|
||||
/*Stop of scatter logger. Record scatterings in what is between start and stop, that is, parabolic feeder. */
|
||||
COMPONENT log_P_stop=Shielding_logger_stop(logger=log_P_start)
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P1_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P1_start = Shielding_log_iterator_Ni_new()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P1_start
|
||||
%}
|
||||
JUMP arm_iter_P1_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P1_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP01=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="NiCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P1_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P1_stop = Shielding_log_iterator_stop(iterator=iter_P1_start)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a11i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P1_start WHEN(MC_GETPAR(iter_P1_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
|
||||
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION2. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P2_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P2_start = Shielding_log_iterator_Ti_new()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P2_start
|
||||
%}
|
||||
JUMP arm_iter_P2_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P2_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP02=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="TiCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P2_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P2_stop = Shielding_log_iterator_stop(iterator=iter_P2_start)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a12i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P2_start WHEN(MC_GETPAR(iter_P2_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR2 SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION3. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P3_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P3_start = Shielding_log_iterator_total()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P3_start
|
||||
%}
|
||||
JUMP arm_iter_P3_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P3_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP03=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="TotalCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P3_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P3_stop = Shielding_log_iterator_stop(iterator=iter_P3_start, last=1)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a13i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P3_start WHEN(MC_GETPAR(iter_P3_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR3 SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**********************************
|
||||
* Optical components within cave *
|
||||
**********************************/
|
||||
@@ -300,8 +446,8 @@ COMPONENT ac_slit = Slit(
|
||||
* Sample area *
|
||||
***************/
|
||||
COMPONENT tof_sample = Monitor_nD(
|
||||
filename = "tof_sample", user1=p_int,
|
||||
options = "x limits=[-0.025 0.025] bins=1000 y limits=[-0.025 0.025] bins=1000 xdiv limits=[-0.75 0.75] bins=150 ydiv limits=[-2.0 2.0] bins=400 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sy limits=[-1 1] bins=2000 user1 list all",
|
||||
filename = "tof_sample",
|
||||
options = "x limits=[-0.025 0.025] bins=1000 y limits=[-0.025 0.025] bins=1000 xdiv limits=[-0.75 0.75] bins=150 ydiv limits=[-2.0 2.0] bins=400 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 list all",
|
||||
xwidth=0.05, yheight = 0.05)
|
||||
WHEN sample==4
|
||||
AT (0, 0, 0) RELATIVE arm_sample_beam
|
||||
@@ -349,67 +495,36 @@ COMPONENT si_sample = Mirror(
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/* Rotate spin 90 degrees*/
|
||||
COMPONENT arm_hack = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
EXTEND %{
|
||||
double tmp = sy;
|
||||
sy = sx;
|
||||
sx = -tmp;
|
||||
%}
|
||||
|
||||
|
||||
|
||||
|
||||
/* flipped arm detector to position analyzer correctly for beam path 1*/
|
||||
COMPONENT arm_detector2 = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_detector
|
||||
ROTATED (0,0,180) RELATIVE arm_detector
|
||||
|
||||
|
||||
|
||||
COMPONENT arm_analyzer = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_detector2
|
||||
ROTATED (-selene_theta+(Theta1_analyzer1-Theta2_analyzer1)/2.0, 0, 0) RELATIVE arm_detector2
|
||||
|
||||
COMPONENT virtual_analyzer_flipper = Arm() // Gone -> Pol_SF_ideal(ny=1, xwidth=1, yheight=1, zdepth=0.001)
|
||||
WHEN (enable_analyzer>2)
|
||||
AT (0, 0, analyzer_flipper_start) RELATIVE arm_analyzer
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer
|
||||
EXTEND %{
|
||||
if(INSTRUMENT_GETPAR(enable_analyzer)>2) sx *=-1;
|
||||
%}
|
||||
AT (0, 0, 0) RELATIVE arm_detector
|
||||
ROTATED (-selene_theta+(Theta1_analyzer1-Theta2_analyzer1)/2.0, 0, 0) RELATIVE arm_detector
|
||||
|
||||
COMPONENT arm_analyzer2 = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_detector2
|
||||
ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector2
|
||||
AT (0, 0, 0) RELATIVE arm_detector
|
||||
ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector
|
||||
|
||||
/* polarization analyser */
|
||||
COMPONENT analyzer1 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer1_start, length=analyzer1_length,
|
||||
COMPONENT analyzer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer1_start, length=analyzer1_length,
|
||||
delta_theta=(Theta1_analyzer1+Theta2_analyzer1+0.05)*PI/180.0, h2=0.14, h1=0.05, abs_ref=0,
|
||||
m_u=5.5, m_d=0.45, both_coated=0, alpha=2.3, W = 0.0014)
|
||||
WHEN enable_analyzer
|
||||
AT (0, 0.0, analyzer1_start) RELATIVE arm_analyzer
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer
|
||||
|
||||
COMPONENT analyzer2 = Polariser(enable_ref=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer2_start, length=analyzer2_length,
|
||||
COMPONENT analyzer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer2_start, length=analyzer2_length,
|
||||
delta_theta=(Theta1_analyzer2+Theta2_analyzer2+0.05)*PI/180.0, h2=0.2, h1=0.05, abs_out=0,
|
||||
m_u=5.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
|
||||
WHEN enable_analyzer==2
|
||||
AT (0, 0.0, analyzer2_start) RELATIVE arm_analyzer2
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer2
|
||||
AT (0, 0.0, analyzer2_start) RELATIVE arm_analyzer
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer
|
||||
|
||||
/* detector */
|
||||
COMPONENT tof_detector = Monitor_nD(
|
||||
filename = "tof_detector",
|
||||
options = "x limits=[-0.25 0.25] bins=1000 y limits=[-0.5 0.5] bins=1000 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sx sy sz user1, list all",
|
||||
xwidth = 0.5, yheight = 1.0)
|
||||
WHEN sample!=4
|
||||
AT (0, 0, detector_arm+0.00001) RELATIVE arm_detector
|
||||
ROTATED (0, 0, 0) RELATIVE arm_detector
|
||||
// COMPONENT tof_detector = Monitor_nD(
|
||||
// filename = "tof_detector",
|
||||
// options = "x limits=[-0.25 0.25] bins=1000 y limits=[-0.5 0.5] bins=1000 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sx limits=[-1 1] bins=1000 sy limits=[-1 1] bins=1000, list all",
|
||||
// xwidth = 0.5, yheight = 1.0)
|
||||
// WHEN sample!=4
|
||||
// AT (0, 0, detector_arm+0.00001) RELATIVE arm_detector
|
||||
|
||||
/***********************************************************************/
|
||||
|
||||
|
||||
@@ -144,7 +144,7 @@ COMPONENT NBOA_window=Al_window(thickness=NBOA_Al_entrance_length)
|
||||
AT (0,0,NBOA_Al_entrance_start) RELATIVE ISCS
|
||||
|
||||
|
||||
COMPONENT NBOA_side = Absorber(xmin=NBOA_side_x, xmax=NBOA_side_x+0.012,
|
||||
COMPONENT NBOA_side = AbsorberFixed(xmin=NBOA_side_x, xmax=NBOA_side_x+0.012,
|
||||
ymin=-0.2, ymax=0.2,
|
||||
zmin=0.0, zmax=E02_01_01_Cu_length+E02_01_01_length+E02_01_02_length+E02_01_03_length)
|
||||
AT (0, 0, E02_01_01_Cu_start) RELATIVE ISCS
|
||||
@@ -181,7 +181,7 @@ COMPONENT NBOA_Cu_collimator=Guide_gravity(
|
||||
* feeder neutron guide *
|
||||
************************/
|
||||
|
||||
COMPONENT E02_01_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_011_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_011_wh/2.0, ymax=E02_01_011_wh/2.0,
|
||||
zmin=0.0, zmax=E02_01_01_length-0.5)
|
||||
AT (0, 0, E02_01_01_start) RELATIVE ISCS
|
||||
@@ -190,7 +190,7 @@ COMPONENT E02_01_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_01_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_012_wh/2.0, ymax=E02_01_012_wh/2.0,
|
||||
zmin=0.5-E02_01_01_Cu_length, zmax=E02_01_01_length)
|
||||
AT (0, 0, E02_01_01_start) RELATIVE ISCS
|
||||
@@ -201,7 +201,7 @@ COMPONENT E02_01_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
%}
|
||||
|
||||
/* Insert part of the feeder before the monolith exit window*/
|
||||
COMPONENT E02_01_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_01 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_01_length-NBOA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_01_01_start, loutyh= NBOA_c*2-E02_01_01_start-E02_01_01_length,
|
||||
linxw = E02_01_01_start, loutxw= NBOA_c*2-E02_01_01_start-E02_01_01_length,
|
||||
@@ -210,7 +210,7 @@ COMPONENT E02_01_01 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, E02_01_01_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT E02_01_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_021_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_021_wh/2.0, ymax=E02_01_021_wh/2.0,
|
||||
zmin=0.0, zmax=E02_01_02_length/2.0)
|
||||
AT (0, 0, E02_01_02_start) RELATIVE ISCS
|
||||
@@ -219,7 +219,7 @@ COMPONENT E02_01_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_01_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_022_wh/2.0, ymax=E02_01_022_wh/2.0,
|
||||
zmin=E02_01_02_length/2.0, zmax=E02_01_02_length)
|
||||
AT (0, 0, E02_01_02_start) RELATIVE ISCS
|
||||
@@ -229,7 +229,7 @@ COMPONENT E02_01_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_01_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_02 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_02_length-NBOA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_01_02_start, loutyh= NBOA_c*2-E02_01_02_start-E02_01_02_length,
|
||||
linxw = E02_01_02_start, loutxw= NBOA_c*2-E02_01_02_start-E02_01_02_length,
|
||||
@@ -238,7 +238,7 @@ COMPONENT E02_01_02 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, E02_01_02_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT E02_01_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_031_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_031_wh/2.0, ymax=E02_01_031_wh/2.0,
|
||||
zmin=0.0, zmax=0.5)
|
||||
AT (0, 0, E02_01_03_start) RELATIVE ISCS
|
||||
@@ -247,7 +247,7 @@ COMPONENT E02_01_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_01_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_032_wh/2.0, ymax=E02_01_032_wh/2.0,
|
||||
zmin=0.5, zmax=1.0)
|
||||
AT (0, 0, E02_01_03_start) RELATIVE ISCS
|
||||
@@ -256,7 +256,7 @@ COMPONENT E02_01_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_01_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_01_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_01_033_wh/2.0, ymax=E02_01_033_wh/2.0,
|
||||
zmin=1.0, zmax=E02_01_03_length)
|
||||
AT (0, 0, E02_01_03_start) RELATIVE ISCS
|
||||
@@ -266,7 +266,7 @@ COMPONENT E02_01_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_01_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_03 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_03_length, dimensionsAt = "mid",
|
||||
linyh = E02_01_03_start, loutyh= NBOA_c*2-E02_01_03_start-E02_01_03_length,
|
||||
linxw = E02_01_03_start, loutxw= NBOA_c*2-E02_01_03_start-E02_01_03_length,
|
||||
@@ -299,7 +299,7 @@ COMPONENT NFGA_Al_window_in=Al_window(
|
||||
WHEN enable_windows
|
||||
AT (0,0,NFGA_Al_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT NFGA_side1 = Absorber(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
|
||||
COMPONENT NFGA_side1 = AbsorberFixed(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
|
||||
ymin=-0.2, ymax=0.2,
|
||||
zmin=0.002, zmax=0.502)
|
||||
AT (0, 0, E02_02_01_start-0.002) RELATIVE ISCS
|
||||
@@ -309,7 +309,7 @@ COMPONENT NFGA_side1 = Absorber(xmin=NFGA_side_x1, xmax=NFGA_side_x1+0.008,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT NFGA_side2 = Absorber(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
|
||||
COMPONENT NFGA_side2 = AbsorberFixed(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
|
||||
ymin=-0.2, ymax=0.2,
|
||||
zmin=0.502, zmax=E02_02_01_length+E02_02_02_length+E02_02_03_length+0.002)
|
||||
AT (0, 0, E02_02_01_start-0.002) RELATIVE ISCS
|
||||
@@ -321,7 +321,7 @@ COMPONENT NFGA_side2 = Absorber(xmin=NFGA_side_x2, xmax=NFGA_side_x2+0.008,
|
||||
%}
|
||||
|
||||
// in-bunker feeder segments
|
||||
COMPONENT E02_02_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_011_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_011_wh/2.0, ymax=E02_02_011_wh/2.0,
|
||||
zmin=0.0, zmax=E02_02_01_length/2.0)
|
||||
AT (0, 0, E02_02_01_start) RELATIVE ISCS
|
||||
@@ -330,7 +330,7 @@ COMPONENT E02_02_011_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_02_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_012_wh/2.0, ymax=E02_02_012_wh/2.0,
|
||||
zmin=E02_02_01_length/2.0, zmax=E02_02_01_length)
|
||||
AT (0, 0, E02_02_01_start) RELATIVE ISCS
|
||||
@@ -340,7 +340,7 @@ COMPONENT E02_02_012_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_01 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_01_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_01_start, loutyh= NFGA_c*2-E02_02_01_start-E02_02_01_length,
|
||||
linxw = E02_02_01_start, loutxw= NFGA_c*2-E02_02_01_start-E02_02_01_length,
|
||||
@@ -349,7 +349,7 @@ COMPONENT E02_02_01 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, E02_02_01_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT E02_02_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_021_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_021_wh/2.0, ymax=E02_02_021_wh/2.0,
|
||||
zmin=0.0, zmax=E02_02_02_length/2.0)
|
||||
AT (0, 0, E02_02_02_start) RELATIVE ISCS
|
||||
@@ -358,7 +358,7 @@ COMPONENT E02_02_021_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_02_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_022_wh/2.0, ymax=E02_02_022_wh/2.0,
|
||||
zmin=E02_02_02_length/2.0, zmax=E02_02_02_length)
|
||||
AT (0, 0, E02_02_02_start) RELATIVE ISCS
|
||||
@@ -368,7 +368,7 @@ COMPONENT E02_02_022_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_02 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_02_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_02_start, loutyh= NFGA_c*2-E02_02_02_start-E02_02_02_length,
|
||||
linxw = E02_02_02_start, loutxw= NFGA_c*2-E02_02_02_start-E02_02_02_length,
|
||||
@@ -377,7 +377,7 @@ COMPONENT E02_02_02 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, E02_02_02_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT E02_02_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_031_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_031_wh/2.0, ymax=E02_02_031_wh/2.0,
|
||||
zmin=0.0, zmax=0.5)
|
||||
AT (0, 0, E02_02_03_start) RELATIVE ISCS
|
||||
@@ -386,7 +386,7 @@ COMPONENT E02_02_031_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_02_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_032_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_032_wh/2.0, ymax=E02_02_032_wh/2.0,
|
||||
zmin=0.5, zmax=1.0)
|
||||
AT (0, 0, E02_02_03_start) RELATIVE ISCS
|
||||
@@ -395,7 +395,7 @@ COMPONENT E02_02_032_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
%}
|
||||
COMPONENT E02_02_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
COMPONENT E02_02_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
ymin=-E02_02_033_wh/2.0, ymax=E02_02_033_wh/2.0,
|
||||
zmin=1.0, zmax=E02_02_03_length)
|
||||
AT (0, 0, E02_02_03_start) RELATIVE ISCS
|
||||
@@ -405,7 +405,7 @@ COMPONENT E02_02_033_wedge = Absorber(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_03 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_03_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_03_start, loutyh= NFGA_c*2-E02_02_03_start-E02_02_03_length,
|
||||
linxw = E02_02_03_start, loutxw= NFGA_c*2-E02_02_03_start-E02_02_03_length,
|
||||
|
||||
+21
-167
@@ -28,80 +28,20 @@ DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
|
||||
// Polarizer parameters used for the detailed design
|
||||
double pol1_xstart = 1.007; // [m] up-stream of focus point
|
||||
double pol1_ystart = 0.0064; // [m] (relative to c-axis)
|
||||
double pol1_xend = 0.365; // [m] up-stream of focus point
|
||||
double pol1_yend = 0.013; // [m]
|
||||
double pol1_gamma = 1.65; // [deg] optimal incident angle for beams passing through focus point
|
||||
double pol1_m = 5.08;
|
||||
double pol1_lin_xstart = 1.007; //[m]
|
||||
double pol1_lin_ystart = 0.0064; //[m]
|
||||
double pol1_lin_xend = 0.812; //[m]
|
||||
double pol1_lin_yend = -0.0014; //[m]
|
||||
double pol1_lin_m = 5.5;
|
||||
|
||||
|
||||
double pol2_xstart = 1.286; // [m] up-stream of focus point
|
||||
double pol2_ystart = 0.009; // [m] (relative to c-axis)
|
||||
double pol2_xend = 0.500; // [m] up-stream of focus point
|
||||
double pol2_yend = 0.0175; // [m]
|
||||
double pol2_gamma = 1.70; // [deg] optimal incident angle for beams passing through focus point
|
||||
double pol2_m = 5.08;
|
||||
double pol2_lin_xstart = 1.286; //[m]
|
||||
double pol2_lin_ystart = 0.009; //[m]
|
||||
double pol2_lin_xend = 1.093; //[m]
|
||||
double pol2_lin_yend = 0.0014; //[m]
|
||||
double pol2_lin_m = 5.08;
|
||||
|
||||
double pol_hclose = 0.12; // [m] height of entrance
|
||||
double pol_hfar = 0.12; // [m] height of exit
|
||||
|
||||
// parameters for polarizer component calculated from the above
|
||||
double pol1_length, pol2_length, pol1_lin_length, pol2_lin_length;
|
||||
double Theta_pol1, Theta_pol2, Theta_rot1, Theta_rot2;
|
||||
double Theta1_pol1, Theta1_pol2, Theta2_pol1, Theta2_pol2;
|
||||
double pol1_lin_rot, pol2_lin_rot;
|
||||
|
||||
double polarizer_max_width = 0.0015; // Maximum sample height to be covered by te polarization analyzers
|
||||
double polarizer_start = 0.3; // Distance to sample to start the first analyzer
|
||||
double polarizer_length = 0.6; // length of first analyzer
|
||||
double Theta1_polarizer, Theta2_polarizer, dist_pol_vfocus; // quantities calculated out of values above and 1.5 degree covered divergence
|
||||
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
dist_pol_vfocus=polarizer_max_width/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
|
||||
|
||||
pol1_length = pol1_xstart-pol1_xend;
|
||||
pol2_length = pol2_xstart-pol2_xend;
|
||||
pol1_lin_length = pol1_lin_xstart-pol1_lin_xend;
|
||||
pol2_lin_length = pol2_lin_xstart-pol2_lin_xend;
|
||||
|
||||
Theta1_pol1=atan2(pol1_ystart, pol1_xstart); // [rad] angle between c-axis and entrance point
|
||||
Theta2_pol1=atan2(pol1_yend, pol1_xend); // [rad] angle between c-axis and exit point
|
||||
Theta1_pol2=atan2(pol2_ystart, pol2_xstart); // [rad] angle between c-axis and entrance point
|
||||
Theta2_pol2=atan2(pol2_yend, pol2_xend); // [rad] angle between c-axis and exit point
|
||||
|
||||
Theta_pol1=(Theta2_pol1-Theta1_pol1)*180.0/PI; // [deg] full covered divergence angle
|
||||
Theta_pol2=(Theta2_pol2-Theta1_pol2)*180.0/PI; // [deg]
|
||||
Theta_rot1=(Theta1_pol1+Theta2_pol1)/2.0*180.0/PI; // [deg] rotation of center of polarizer
|
||||
Theta_rot2=(Theta1_pol2+Theta2_pol2)/2.0*180.0/PI; // [deg]
|
||||
|
||||
pol1_lin_rot=atan2(pol1_lin_ystart-pol1_lin_yend,
|
||||
pol1_lin_xstart-pol1_lin_xend)*180.0/PI; // [deg] rotation angle of polarizer
|
||||
pol2_lin_rot=atan2(pol2_lin_ystart-pol2_lin_yend,
|
||||
pol2_lin_xstart-pol2_lin_xend)*180.0/PI; // [deg] rotation angle of polarizer
|
||||
|
||||
|
||||
|
||||
printf(" Polarizer 1 angle = %.2f deg\n", Theta_rot1);
|
||||
printf(" Polarizer 1 divergence = %.2f deg\n", Theta_pol1);
|
||||
printf(" Polarizer 1 distance = %.1f mm\n", pol1_xstart*1000.0);
|
||||
printf(" Polarizer 1 length = %.1f mm\n", pol1_length*1000.0);
|
||||
printf(" PolStraight 1 angle = %.2f deg\n", pol1_lin_rot);
|
||||
|
||||
printf(" Polarizer 2 angle = %.2f deg\n", Theta_rot1);
|
||||
printf(" Polarizer 2 divergence = %.2f deg\n", Theta_pol2);
|
||||
printf(" Polarizer 2 distance = %.1f mm\n", pol2_xstart*1000.0);
|
||||
printf(" Polarizer 2 length = %.1f mm\n", pol2_length*1000.0);
|
||||
printf(" PolStraight 2 angle = %.2f deg\n", pol2_lin_rot);
|
||||
Theta1_polarizer=atan((dist_pol_vfocus+polarizer_start)/dist_pol_vfocus*polarizer_max_width/2.0/polarizer_start)*180.0/PI;
|
||||
Theta2_polarizer=atan((dist_pol_vfocus+polarizer_start+polarizer_length)/dist_pol_vfocus*polarizer_max_width/2.0/(polarizer_start+polarizer_length))*180.0/PI;
|
||||
|
||||
%}
|
||||
TRACE
|
||||
@@ -127,112 +67,26 @@ REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0
|
||||
/**************************************
|
||||
* Middle focus between Selene guides *
|
||||
**************************************/
|
||||
COMPONENT arm_polref = Arm()
|
||||
COMPONENT arm_polarizer = Arm()
|
||||
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
|
||||
ROTATED (selene_theta, -selene_theta, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT arm_polarizer = Arm()
|
||||
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
|
||||
ROTATED (0, 0, 0) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT arm_pol1 = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_polarizer
|
||||
ROTATED (0, -Theta_rot1, 0) RELATIVE arm_polarizer
|
||||
|
||||
COMPONENT arm_pol2 = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_polarizer
|
||||
ROTATED (0, -Theta_rot2, 0) RELATIVE arm_polarizer
|
||||
COMPONENT polarizer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=-(polarizer_start+polarizer_length), length=polarizer_length,
|
||||
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
|
||||
WHEN enable_polarizer
|
||||
AT (0, 0, -(polarizer_start+polarizer_length)) RELATIVE arm_polarizer
|
||||
ROTATED (0,0,90.0) RELATIVE arm_polarizer
|
||||
|
||||
|
||||
COMPONENT polarizer2_lin = Pol_mirror(zwidth = pol2_lin_length, yheight = pol_hfar, p_reflect=0,
|
||||
rUpData="SNAG_m5p5_Tup.ref", rDownData="SNAG_m5p5_Tdown.ref")
|
||||
WHEN (enable_polarizer>0)
|
||||
AT ((pol2_lin_ystart+pol2_lin_yend)/2.0, 0, -(pol2_lin_xstart+pol2_lin_xend)/2.0) RELATIVE arm_polarizer
|
||||
ROTATED (0,-pol2_lin_rot,0.0) RELATIVE arm_polarizer
|
||||
GROUP polarizer2_set
|
||||
EXTEND
|
||||
%{
|
||||
if (SCATTERED==2) {
|
||||
p_int +=2;
|
||||
}
|
||||
%}
|
||||
|
||||
COMPONENT polarizer2 = Polariser(lin=-pol2_xstart, length=pol2_length,
|
||||
enable_ref=1, abs_ref=1, abs_out=1, both_coated=1,
|
||||
d_substrate = 5e-4, T_loss=4.0e3,
|
||||
m_u=pol2_m, m_d=0.6, m_residual=0.55,
|
||||
alpha=2.3, W = 0.0014, reflect_d=0, reflect_u=0,
|
||||
delta_theta=Theta_pol2*PI/180.0,
|
||||
h2=pol_hfar, h1=pol_hclose)
|
||||
WHEN (enable_polarizer>0)
|
||||
AT (0, 0, -pol2_xstart) RELATIVE arm_pol2
|
||||
ROTATED (0,0,90.0) RELATIVE arm_pol2
|
||||
GROUP polarizer2_set
|
||||
EXTEND
|
||||
%{
|
||||
if (SCATTERED) {
|
||||
ALLOW_BACKPROP;
|
||||
PROP_Z0;
|
||||
p_int +=1;
|
||||
} else ABSORB;
|
||||
%}
|
||||
|
||||
COMPONENT replacement_pol2_set = Slit(xwidth=0.2, yheight=0.2)
|
||||
WHEN (enable_polarizer==0)
|
||||
AT (0, 0, -0.2) RELATIVE arm_polarizer
|
||||
GROUP polarizer2_set
|
||||
|
||||
COMPONENT polarizer1_lin = Pol_mirror(zwidth = pol1_lin_length, yheight = pol_hfar, p_reflect=0,
|
||||
rUpData="SNAG_m5p5_Tup.ref", rDownData="SNAG_m5p5_Tdown.ref")
|
||||
WHEN (enable_polarizer>0)
|
||||
AT ((pol1_lin_ystart+pol1_lin_yend)/2.0, 0, -(pol1_lin_xstart+pol1_lin_xend)/2.0) RELATIVE arm_polarizer
|
||||
ROTATED (0,-pol1_lin_rot,0.0) RELATIVE arm_polarizer
|
||||
GROUP polarizer1_set
|
||||
EXTEND
|
||||
%{
|
||||
if (SCATTERED==2) {
|
||||
p_int +=8;
|
||||
}
|
||||
%}
|
||||
|
||||
COMPONENT polarizer1 = Polariser(lin=-pol1_xstart, length=pol1_length,
|
||||
enable_ref=1, abs_ref=1, abs_out=1, both_coated=1,
|
||||
d_substrate = 5e-4, T_loss=4.0e3,
|
||||
m_u=pol1_m, m_d=0.6, m_residual=0.55,
|
||||
alpha=2.3, W = 0.0014, reflect_d=0, reflect_u=0,
|
||||
delta_theta=Theta_pol1*PI/180.0,
|
||||
h2=pol_hfar, h1=pol_hclose)
|
||||
WHEN (enable_polarizer>0)
|
||||
AT (0, 0, -pol1_xstart) RELATIVE arm_pol1
|
||||
ROTATED (0,0,90.0) RELATIVE arm_pol1
|
||||
GROUP polarizer1_set
|
||||
EXTEND
|
||||
%{
|
||||
if (SCATTERED) {
|
||||
p_int +=4;
|
||||
}
|
||||
%}
|
||||
|
||||
COMPONENT replacement_pol1_set = Slit(xwidth=0.2, yheight=0.2)
|
||||
WHEN (enable_polarizer==0)
|
||||
AT (0, 0, -0.1) RELATIVE arm_polarizer
|
||||
GROUP polarizer1_set
|
||||
/* The proximal polariser comes next */
|
||||
COMPONENT polarizer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=polarizer_start, length=polarizer_length,
|
||||
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
|
||||
WHEN (enable_polarizer>1)
|
||||
AT (0, 0, polarizer_start) RELATIVE arm_polarizer
|
||||
ROTATED (0.0,0,90.0) RELATIVE arm_polarizer
|
||||
|
||||
|
||||
|
||||
|
||||
COMPONENT mf_divergence_h = DivPos_monitor(nb=21, ndiv=41, filename="mf_hordiv",
|
||||
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0)
|
||||
AT (0, 0, 0) RELATIVE arm_polarizer
|
||||
|
||||
COMPONENT mf_divergence_v = DivPos_monitor(nb=21, ndiv=41, filename="mf_verdiv",
|
||||
xmin=-0.02, xmax=0.02, ymin=-0.02, ymax=0.02, maxdiv=2.0, vertical=1)
|
||||
AT (0, 0, 0) RELATIVE arm_polarizer
|
||||
|
||||
|
||||
|
||||
COMPONENT virtual_flipper = Pol_SF_ideal(ny=1, xwidth=0.1, yheight=0.1, zdepth=0.001)
|
||||
WHEN (enable_polarizer>2)
|
||||
AT (0, 0, 0) RELATIVE arm_polarizer
|
||||
|
||||
|
||||
FINALLY
|
||||
|
||||
@@ -102,7 +102,7 @@ REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0
|
||||
/* Absorber to cut direct view beam (Copper in CAD model) */
|
||||
COMPONENT slit_before_selene_guide_1 = Slit(
|
||||
xmin = 0, xmax = selene_entry+0.001,
|
||||
ymin = -selene_entry-0.005, ymax = 0)
|
||||
ymin = 0, ymax = selene_entry+0.005)
|
||||
AT (0, 0, selene_distance-0.01) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT block_before_selene_guide_1 = Absorber(
|
||||
@@ -113,145 +113,145 @@ COMPONENT block_before_selene_guide_1 = Absorber(
|
||||
AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
|
||||
|
||||
/* Selene 1 elliptic guide */
|
||||
COMPONENT E02_03_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_01 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
|
||||
linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_01, mbottom=scoating_01, mtop=0,
|
||||
mright=0, mleft=scoating_01, mtop=scoating_01, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_02 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
|
||||
linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_02, mbottom=scoating_02, mtop=0,
|
||||
mright=0, mleft=scoating_02, mtop=scoating_02, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_03 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
|
||||
linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_03, mbottom=scoating_03, mtop=0,
|
||||
mright=0, mleft=scoating_03, mtop=scoating_03, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_04 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_04 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
|
||||
linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_04, mbottom=scoating_04, mtop=0,
|
||||
mright=0, mleft=scoating_04, mtop=scoating_04, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_05 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_05 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
|
||||
linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_05, mbottom=scoating_05, mtop=0,
|
||||
mright=0, mleft=scoating_05, mtop=scoating_05, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_06 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_06 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
|
||||
linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_06, mbottom=scoating_06, mtop=0,
|
||||
mright=0, mleft=scoating_06, mtop=scoating_06, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_07 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_07 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
|
||||
linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_07, mbottom=scoating_07, mtop=0,
|
||||
mright=0, mleft=scoating_07, mtop=scoating_07, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_08 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_08 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
|
||||
linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_08, mbottom=scoating_08, mtop=0,
|
||||
mright=0, mleft=scoating_08, mtop=scoating_08, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_09 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_09 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
|
||||
linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_09, mbottom=scoating_09, mtop=0,
|
||||
mright=0, mleft=scoating_09, mtop=scoating_09, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_10 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_10 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
|
||||
linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_10, mbottom=scoating_10, mtop=0,
|
||||
mright=0, mleft=scoating_10, mtop=scoating_10, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_11 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_11 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
|
||||
linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_11, mbottom=scoating_11, mtop=0,
|
||||
mright=0, mleft=scoating_11, mtop=scoating_11, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_12 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_12 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
|
||||
linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_12, mbottom=scoating_12, mtop=0,
|
||||
mright=0, mleft=scoating_12, mtop=scoating_12, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_13 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_13 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
|
||||
linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_13, mbottom=scoating_13, mtop=0,
|
||||
mright=0, mleft=scoating_13, mtop=scoating_13, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_14 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_14 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
|
||||
linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_14, mbottom=scoating_14, mtop=0,
|
||||
mright=0, mleft=scoating_14, mtop=scoating_14, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_15 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_03_15 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
|
||||
linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=0, mleft=scoating_15, mbottom=scoating_15, mtop=0,
|
||||
mright=0, mleft=scoating_15, mtop=scoating_15, mbottom=0,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
/* Absorber to cut direct view beam (Copper in CAD model) */
|
||||
COMPONENT slit_after_selene_guide_1 = Slit(
|
||||
xmin = 0, xmax = selene_entry+0.005,
|
||||
ymin = -selene_entry-0.005, ymax = 0)
|
||||
ymin = 0, ymax = selene_entry+0.01)
|
||||
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT block_after_selene_guide_1 = Absorber(
|
||||
|
||||
@@ -42,7 +42,7 @@ TRACE
|
||||
/* Absorber to cut direct view beam (Bor-Al in CAD model) */
|
||||
COMPONENT slit_before_selene_guide_2 = Slit(
|
||||
xmin = -selene_entry-0.005, xmax=0,
|
||||
ymin = 0, ymax = selene_entry+0.01)
|
||||
ymin = -selene_entry-0.01, ymax = 0)
|
||||
AT (0, 0, selene_distance-0.002) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT block_before_selene_guide_2 = Absorber(
|
||||
@@ -54,145 +54,145 @@ COMPONENT block_before_selene_guide_2 = Absorber(
|
||||
|
||||
|
||||
/* Selene 2 elliptic guide first half */
|
||||
COMPONENT E02_04_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_01 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
|
||||
linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_15, mleft=0, mbottom=0, mtop=scoating_15,
|
||||
mright=scoating_15, mleft=0, mtop=0, mbottom=scoating_15,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_02 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
|
||||
linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_14, mleft=0, mbottom=0, mtop=scoating_14,
|
||||
mright=scoating_14, mleft=0, mtop=0, mbottom=scoating_14,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_03 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
|
||||
linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_13, mleft=0, mbottom=0, mtop=scoating_13,
|
||||
mright=scoating_13, mleft=0, mtop=0, mbottom=scoating_13,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_04 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_04 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
|
||||
linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_12, mleft=0, mbottom=0, mtop=scoating_12,
|
||||
mright=scoating_12, mleft=0, mtop=0, mbottom=scoating_12,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_05 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_05 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
|
||||
linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_11, mleft=0, mbottom=0, mtop=scoating_11,
|
||||
mright=scoating_11, mleft=0, mtop=0, mbottom=scoating_11,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_06 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_06 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
|
||||
linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_10, mleft=0, mbottom=0, mtop=scoating_10,
|
||||
mright=scoating_10, mleft=0, mtop=0, mbottom=scoating_10,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_07 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_07 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
|
||||
linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_09, mleft=0, mbottom=0, mtop=scoating_09,
|
||||
mright=scoating_09, mleft=0, mtop=0, mbottom=scoating_09,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_08 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_08 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
|
||||
linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_08, mleft=0, mbottom=0, mtop=scoating_08,
|
||||
mright=scoating_08, mleft=0, mtop=0, mbottom=scoating_08,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_09 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_09 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
|
||||
linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_07, mleft=0, mbottom=0, mtop=scoating_07,
|
||||
mright=scoating_07, mleft=0, mtop=0, mbottom=scoating_07,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_10 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_10 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
|
||||
linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_06, mleft=0, mbottom=0, mtop=scoating_06,
|
||||
mright=scoating_06, mleft=0, mtop=0, mbottom=scoating_06,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_11 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_11 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
|
||||
linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_05, mleft=0, mbottom=0, mtop=scoating_05,
|
||||
mright=scoating_05, mleft=0, mtop=0, mbottom=scoating_05,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_12 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_12 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
|
||||
linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_04, mleft=0, mbottom=0, mtop=scoating_04,
|
||||
mright=scoating_04, mleft=0, mtop=0, mbottom=scoating_04,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_13 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_13 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
|
||||
linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_03, mleft=0, mbottom=0, mtop=scoating_03,
|
||||
mright=scoating_03, mleft=0, mtop=0, mbottom=scoating_03,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_14 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_14 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
|
||||
linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_02, mleft=0, mbottom=0, mtop=scoating_02,
|
||||
mright=scoating_02, mleft=0, mtop=0, mbottom=scoating_02,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_15 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_04_15 = Elliptic_guide_gravity_custom(
|
||||
l=selene_segment, dimensionsAt = "mid",
|
||||
linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
|
||||
linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
|
||||
xwidth=selene_b*2, yheight=selene_b*2,
|
||||
mright=scoating_01, mleft=0, mbottom=0, mtop=scoating_01,
|
||||
mright=scoating_01, mleft=0, mtop=0, mbottom=scoating_01,
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
/* Absorber to cut direct view beam (Bor-Al in CAD model) */
|
||||
COMPONENT slit_after_selene_guide_2 = Slit(
|
||||
xmin = -selene_entry-0.001, xmax=0, xmax = 0.0,
|
||||
ymin = 0, ymax = selene_entry+0.01)
|
||||
ymin = -selene_entry-0.005, ymax = 0)
|
||||
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT block_after_selene_guide_2 = Absorber(
|
||||
|
||||
+515
-697
File diff suppressed because it is too large
Load Diff
@@ -1,58 +0,0 @@
|
||||
4.37173E-05 0.253838498
|
||||
0.001967279 0.548126562
|
||||
0.004240578 0.781912929
|
||||
0.006513878 0.831781609
|
||||
0.008656026 0.935568791
|
||||
0.010841891 0.996046699
|
||||
0.012984039 0.994254859
|
||||
0.015038752 0.780745054
|
||||
0.017443204 0.571166135
|
||||
0.019497917 0.421992388
|
||||
0.021771217 0.347878292
|
||||
0.023913365 0.314863976
|
||||
0.026317816 0.273648917
|
||||
0.028241378 0.248910367
|
||||
0.030514678 0.229613667
|
||||
0.032700543 0.210315136
|
||||
0.034711539 0.201079179
|
||||
0.03711599 0.179153501
|
||||
0.039301855 0.174572119
|
||||
0.041444003 0.161570874
|
||||
0.043629868 0.159330731
|
||||
0.045772016 0.158627627
|
||||
0.047957881 0.133760075
|
||||
0.050231181 0.128230852
|
||||
0.052373329 0.135460425
|
||||
0.054559194 0.125901964
|
||||
0.056701342 0.120064418
|
||||
0.058974642 0.122922579
|
||||
0.061160507 0.127084935
|
||||
0.063171503 0.12641248
|
||||
0.06548852 0.117637177
|
||||
0.067630668 0.105706827
|
||||
0.069903968 0.112415362
|
||||
0.072089833 0.119429932
|
||||
0.074231981 0.105676178
|
||||
0.076417846 0.107716872
|
||||
0.078428842 0.107129046
|
||||
0.080833293 0.10262909
|
||||
0.082888007 0.094164854
|
||||
0.085161306 0.092239896
|
||||
0.087303454 0.085655788
|
||||
0.089620471 0.081615114
|
||||
0.091631467 0.076984785
|
||||
0.093817332 0.076157254
|
||||
0.096090632 0.070972035
|
||||
0.09823278 0.060551277
|
||||
0.100287493 0.071366358
|
||||
0.102691945 0.065426801
|
||||
0.104834093 0.059045344
|
||||
0.107019958 0.056229269
|
||||
0.109162106 0.061290519
|
||||
0.111347971 0.051160701
|
||||
0.113490119 0.053836338
|
||||
0.115763418 0.05839439
|
||||
0.117949284 0.057514252
|
||||
0.120091431 0.047536765
|
||||
0.122364731 0.04877829
|
||||
0.300000000 0.000000000
|
||||
@@ -1,58 +0,0 @@
|
||||
4.37173E-05 0.224586133
|
||||
0.001967279 0.518057252
|
||||
0.004240578 0.775672976
|
||||
0.006513878 0.825879049
|
||||
0.008656026 0.899985882
|
||||
0.010841891 0.967789786
|
||||
0.012984039 0.99196296
|
||||
0.015038752 0.992909754
|
||||
0.017443204 0.994177124
|
||||
0.019497917 0.991424846
|
||||
0.021771217 0.993887344
|
||||
0.023913365 0.998648926
|
||||
0.026317816 0.999361915
|
||||
0.028241378 0.998678679
|
||||
0.030514678 0.998109885
|
||||
0.032700543 0.997772638
|
||||
0.034711539 0.997392554
|
||||
0.03711599 0.999600524
|
||||
0.039301855 0.994779551
|
||||
0.041444003 0.996274647
|
||||
0.043629868 0.996149207
|
||||
0.045772016 0.998305774
|
||||
0.047957881 0.99658353
|
||||
0.050231181 0.995950077
|
||||
0.052373329 0.995563046
|
||||
0.054559194 0.995479284
|
||||
0.056701342 0.993347667
|
||||
0.058974642 0.990508095
|
||||
0.061160507 0.990476431
|
||||
0.063171503 0.988754593
|
||||
0.06548852 0.987959462
|
||||
0.067630668 0.987058572
|
||||
0.069903968 0.988939946
|
||||
0.072089833 0.986821236
|
||||
0.074231981 0.985149296
|
||||
0.076417846 0.987362592
|
||||
0.078428842 0.986507086
|
||||
0.080833293 0.98295765
|
||||
0.082888007 0.983313364
|
||||
0.085161306 0.983817788
|
||||
0.087303454 0.983554288
|
||||
0.089620471 0.979763465
|
||||
0.091631467 0.975348682
|
||||
0.093817332 0.975887375
|
||||
0.096090632 0.972236927
|
||||
0.09823278 0.965515882
|
||||
0.100287493 0.962852402
|
||||
0.102691945 0.956266464
|
||||
0.104834093 0.953509902
|
||||
0.107019958 0.945236222
|
||||
0.109162106 0.808306574
|
||||
0.111347971 0.186001557
|
||||
0.113490119 0.078929784
|
||||
0.115763418 0.051976893
|
||||
0.117949284 0.05338487
|
||||
0.120091431 0.04844847
|
||||
0.122364731 0.053399168
|
||||
0.300000000 0.000000000
|
||||
Binary file not shown.
@@ -0,0 +1,93 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright (C) 1997-2011, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: ScanningSlit
|
||||
*
|
||||
* %I
|
||||
* Written by: Jochen Stahn, based on the Slip component by
|
||||
* Kim Lefmann and Henrik M. Roennow
|
||||
* Date: 15. 10. 2013
|
||||
* Version: $Revision: 0.00 $
|
||||
* Origin: PSI
|
||||
* Release: McStas 1.12c
|
||||
*
|
||||
* moving rectangular slit with optional insignificance cut
|
||||
*
|
||||
* %D
|
||||
* A simple rectangular slit, where the blades in x-direction move with time.
|
||||
* No transmission around the slit is allowed.
|
||||
* If cutting option is used, low-weight neutron rays are ABSORBED
|
||||
*
|
||||
* Example: FastSlit(height=0.1)
|
||||
*
|
||||
* %P
|
||||
* INPUT PARAMETERS
|
||||
*
|
||||
* height: height of the opening in y, centered.
|
||||
*
|
||||
* Optional parameters:
|
||||
* cut: Lower limit for allowed weight (1)
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
|
||||
DEFINE COMPONENT ScanningSlit
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS (xmin=0, xmax=0, ymin=0, ymax=0, height=0.1, dist=2.4, lmin=4, lmax=10, omeg=0, eps=1.25, reso=0.01, dthet=1.5, posi=33.7)
|
||||
OUTPUT PARAMETERS ()
|
||||
/* STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) */
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
double lact;
|
||||
%}
|
||||
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (height > 0) { ymax=height/2; ymin=-ymax; }
|
||||
lmax *= (1+reso/2) ;
|
||||
lmin *= (1-reso/2) ;
|
||||
eps *= PI/180 ;
|
||||
omeg *= PI/180 ;
|
||||
dthet *= PI/180;
|
||||
%}
|
||||
|
||||
TRACE
|
||||
%{
|
||||
PROP_Z0;
|
||||
lact = 3.956e3*t/posi ;
|
||||
xmin = -dist * (1+reso/2)*( eps + 0.5*dthet - dthet*(lact-lmin)/(lmax-lmin) ) ;
|
||||
xmax = -dist * (1-reso/2)*( eps + 0.5*dthet - dthet*(lact-lmin)/(lmax-lmin) ) ;
|
||||
if ((x<xmin || x>xmax || y<ymin || y>ymax))
|
||||
ABSORB;
|
||||
else
|
||||
SCATTER;
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
magnify("xy");
|
||||
double xw, yh;
|
||||
xw = (xmax - xmin)/2.0;
|
||||
yh = (ymax - ymin)/2.0;
|
||||
multiline(3, xmin-xw, (double)ymax, 0.0,
|
||||
(double)xmin, (double)ymax, 0.0,
|
||||
(double)xmin, ymax+yh, 0.0);
|
||||
multiline(3, xmax+xw, (double)ymax, 0.0,
|
||||
(double)xmax, (double)ymax, 0.0,
|
||||
(double)xmax, ymax+yh, 0.0);
|
||||
multiline(3, xmin-xw, (double)ymin, 0.0,
|
||||
(double)xmin, (double)ymin, 0.0,
|
||||
(double)xmin, ymin-yh, 0.0);
|
||||
multiline(3, xmax+xw, (double)ymin, 0.0,
|
||||
(double)xmax, (double)ymin, 0.0,
|
||||
(double)xmax, ymin-yh, 0.0);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,145 @@
|
||||
/*******************************************************************************
|
||||
* McStas instrument definition URL=http://www.mcstas.org
|
||||
*
|
||||
* Instrument: Selene_geometry
|
||||
*
|
||||
* %Identification
|
||||
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
|
||||
* Date: 01. 11. 2017
|
||||
* Origin: PSI
|
||||
* Release: McStas 2.4.1
|
||||
* Version: 1.0
|
||||
* %INSTRUMENT_SITE: Estia (E02)
|
||||
*
|
||||
* Estia is a vertical sample, focusing reflectometer for small sample
|
||||
*
|
||||
* %Description
|
||||
* This is a test model to analyze the impact of the displacement of the two Selene guide
|
||||
* mirrors to each other and the virtual source. With this it is possible to evaluate
|
||||
* necessary engineering requirements for the guide carriers.
|
||||
*
|
||||
* %Parameters
|
||||
* sample_length: [m] Size of sample in beam direction, also controls virtual source opening
|
||||
* sample_height: [m] Size of sample in vertical direction, also controls virtual source opening
|
||||
* omegaa: [deg] sample rotation omega
|
||||
*
|
||||
* %End
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE INSTRUMENT Selene_geometry(double sample_length=0.01, double sample_height=0.01, omegaa=1.2,
|
||||
double tx_1=0.0, double tz_1=0.0, double rz_1=0.0, double ry_1=0.0,
|
||||
double tx_2=0.0, double tz_2=0.0, double rz_2=0.0, double ry_2=0.0
|
||||
)
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
|
||||
%}
|
||||
TRACE
|
||||
|
||||
COMPONENT origin = Progress_bar()
|
||||
AT (0,0,0) ABSOLUTE
|
||||
|
||||
COMPONENT ISCS = Arm()
|
||||
AT (0, 0, 0) RELATIVE origin
|
||||
|
||||
COMPONENT arm_selene1 = Arm()
|
||||
AT (tx_1, 0, tz_1) RELATIVE ISCS
|
||||
ROTATED (0, ry_1, rz_1) RELATIVE ISCS
|
||||
|
||||
COMPONENT arm_selene2 = Arm()
|
||||
AT (tx_2, 0, 2*selene_c+tz_2) RELATIVE ISCS
|
||||
ROTATED (0, ry_2, rz_2) RELATIVE ISCS
|
||||
|
||||
COMPONENT arm_virtual_source_beam = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_selene1
|
||||
ROTATED (0, 1.25, 0) RELATIVE ISCS
|
||||
// Axes used for the virtual source masks, this depends on the sample angle to achieve the correct footprint
|
||||
COMPONENT arm_virtual_source = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_virtual_source_beam
|
||||
ROTATED (0, omegaa, 0) RELATIVE arm_virtual_source_beam
|
||||
|
||||
|
||||
COMPONENT source = Moderator(radius = 0.01, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002,
|
||||
target_index=3, Emin=0.75, Emax=6)
|
||||
AT (0,0,-0.1) RELATIVE ISCS
|
||||
ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
|
||||
|
||||
/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
|
||||
COMPONENT virtual_source_TR = Slit(
|
||||
xmin = 0.0, xmax = 1.0, ymin = -1.0, ymax = sample_height/2)
|
||||
AT (0, 0, -0.5*sample_length) RELATIVE arm_virtual_source
|
||||
|
||||
/* The actual virtual source mask, two L-shaped absorbers (second bottom-left) */
|
||||
COMPONENT virtual_source_BL = Slit(
|
||||
xmin = -1.0, xmax = 0.0, ymin = -sample_height/2, ymax = 1.0)
|
||||
AT (0, 0, 0.5*sample_length) RELATIVE arm_virtual_source
|
||||
|
||||
%include "Estia_selene.instr"
|
||||
|
||||
COMPONENT M_n3 = Monitor_nD(
|
||||
filename = "M_n3",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c-0.01) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_n2 = Monitor_nD(
|
||||
filename = "M_n2",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c-0.005) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_n1 = Monitor_nD(
|
||||
filename = "M_n1",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c-0.0025) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_0 = Monitor_nD(
|
||||
filename = "M_0",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT MFP = Monitor_nD(
|
||||
filename = "MFP",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c+0.0001) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa-90, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_p1 = Monitor_nD(
|
||||
filename = "M_p1",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c+0.0025) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_p2 = Monitor_nD(
|
||||
filename = "M_p2",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c+0.005) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
COMPONENT M_p3 = Monitor_nD(
|
||||
filename = "M_p3",
|
||||
options = "x limits=[-0.05 0.05] bins=1000 y limits=[-0.05 0.05] bins=1000 lambda limits=[0 35] bins=350",
|
||||
xwidth = 0.5, yheight = 0.5)
|
||||
AT (0, 0, 4*selene_c+0.01) RELATIVE ISCS
|
||||
ROTATED (0, 1.25+omegaa, 0) RELATIVE ISCS
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,949 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Shielding_calculator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Rodion Kolevatov
|
||||
* Date: August 2018
|
||||
* Version: $Revision: 1.0 $
|
||||
* Release: ---
|
||||
* Origin: IFE
|
||||
*
|
||||
* Calculating lateral shielding thickness required to attenuate dose rate from gammas
|
||||
* arising from coating capture to a level specified by MaxRate.
|
||||
*
|
||||
* %D
|
||||
* Start of the trace-region in which SCATTER events should be logged.
|
||||
* Whenever a SCATTER occurs in components between this one and its counterpart Scatter_logger_stop the neutron state is logged.
|
||||
* The log is kept in memory - but only for one single
|
||||
* numerical neutron at a time, so there should be no real danger of memory overrun.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
* Allowed dose rate beyond shielding, MaxRate
|
||||
* Names of text files where the coating capture rates are recorded
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_calculator
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS (double MaxRate = 3.0, double Innerspace = 0.3, string NiCaptureFile = "NiCapture.dat",
|
||||
string TiCaptureFile = "TiCapture.dat", string TotalCaptureFile = "TotalCapture.dat", string OutputFile="Shielding.dat")
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef LIN_INT_ROUTINE
|
||||
#define LIN_INT_ROUTINE 1
|
||||
//Multi-d linear interpolation routine. Array of args, number of args, sizes of args in interpolation grid, argument grid in single line, data in single line
|
||||
double lint (double * args, int Narg, int * sizearg, double * argtable, double * datatable)
|
||||
{
|
||||
if (Narg==1) {
|
||||
// printf ("Narg=1, arg = %g\n", *args);
|
||||
int point=0;
|
||||
if ((*args)<(*argtable) ) return *datatable;
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return * (datatable+ *sizearg-1); // if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) point++;
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup = (*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*(*(datatable+point))+weightlow*(*(datatable+point-1));
|
||||
}//arg within range of values
|
||||
}// if Narg==1
|
||||
else if (Narg >1){//if more than one argument
|
||||
//lookup how large is the data (Narg-1)D slice for fixed value of the first argument
|
||||
int slicesize=1;
|
||||
int i,point=0;
|
||||
for ( i=1;i<Narg;i++) slicesize*=sizearg[i];
|
||||
// printf("Narg > 1, slicesize = %d, arggridstart = %g, argument = %g, arggridend = %g\n",slicesize,*argtable,*args, *(argtable+ *sizearg-1) );
|
||||
//search weights for the first argument and get results with one argument less
|
||||
if ((*args)<(*argtable) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable);
|
||||
// if the value is too low -- return what is on the lower bound
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), (datatable+ slicesize*(*sizearg)));
|
||||
// if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) {
|
||||
// printf ("*(argtable+point) = %g\n", *(argtable+point));
|
||||
point++;}
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup =(*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+point*slicesize )
|
||||
+weightlow*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+(point-1)*slicesize);
|
||||
} // argument within the ragne
|
||||
} //Narg>1
|
||||
};
|
||||
#endif
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
#include <dirent.h>
|
||||
#include <unistd.h>
|
||||
/******Parameters relevant for calculations of the shielding*******/
|
||||
#define stepPb 0.001 // Step for calculations with steel shielding (2mm)
|
||||
#define stepFe 0.001 // Step for calculations with steel shielding (5mm)
|
||||
#define stepbt 0.005 // Step for calculations with concrete shielding. (1cm)
|
||||
#define Rtubing 0.05 //in meters. Vacuum tubing inner radius, used for calculation of combined shielding in Febt case.
|
||||
#define TFe 0.1 // in meters. Thickness of steel vacuum tubing. (Assuming guide model as guide in evacuated steel pipe with wall thickness TFe).
|
||||
#define stepcomb 0.005 // Step for calculation of combined shielding, TFe layer of steel plus concrete. (2.5 cm)
|
||||
|
||||
//*** NEUTRON CONVERSION DATA *****//
|
||||
#ifndef NEUTRON_GAMMA_DATA
|
||||
#define NEUTRON_GAMMA_DATA 1
|
||||
//Photons per capture in NiTi coating
|
||||
double fraction[]={0.003565789, 0.009694611, 0.028838269, 0.149117282, 0.117348519, 0.042269932, 0.386434136, 0.092896131, 0.046767054, 0.032512648, 0.043547581, 0.040537585, 0.510133557, 0.320364465};
|
||||
double energy [] = {0.150,
|
||||
0.200,
|
||||
0.300,
|
||||
0.400,
|
||||
0.500,
|
||||
0.600,
|
||||
0.800,
|
||||
1.,
|
||||
1.50,
|
||||
2.000,
|
||||
3.0,
|
||||
4.0,
|
||||
5.000,
|
||||
6.000,
|
||||
7.0,
|
||||
8.000,
|
||||
9.0,
|
||||
10.000,
|
||||
11.000};
|
||||
int nEgroup= 19;
|
||||
//Spectrum of capture photons from Ni and Ti
|
||||
//Ni gamma capture specrum
|
||||
double fractionNi[]={0.003722,//150
|
||||
0.015513,//200
|
||||
0.048432,//300
|
||||
0.047485,//400
|
||||
0.198059,//500
|
||||
0.002809,//600
|
||||
0.003142,//800
|
||||
0.06579,//1000
|
||||
0.054534,//1500
|
||||
0.028872,//2000
|
||||
0.063882,//3000
|
||||
0.050286,//4000
|
||||
0.035567,//5000
|
||||
0.080326,//6000
|
||||
0.156571,//7000
|
||||
0.130355,//8000
|
||||
0.534153,//9000
|
||||
0.000558,//10000
|
||||
0.000260};//11000
|
||||
|
||||
|
||||
//Ti capture gamma spectrum
|
||||
double fractionTi[]={0.009326, //0.15
|
||||
0.0, //0.2
|
||||
0.001938,//0.3
|
||||
0.302632,//0.4
|
||||
0.000882,//0.5
|
||||
0.000668, //0.6
|
||||
0.000914, //0.8
|
||||
0.023046,//1.0
|
||||
0.947104,//1.5
|
||||
0.232257,//2.0
|
||||
0.086418,//3.0
|
||||
0.162878,//4.0
|
||||
0.11209,//5.0
|
||||
0.013419,//6.0
|
||||
0.0875051, //7.0
|
||||
0.009083,//8.0
|
||||
0.001923, //9.0
|
||||
0.002189,//10.0
|
||||
0.000332};//11.0
|
||||
//double fractionB[]={0.0, 0.0, 0.0, 0.0, 0.93, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0};
|
||||
double fractionB[]={0,//0.15
|
||||
0,//0.2
|
||||
0,//0.3
|
||||
2.32983E-06,//0.4
|
||||
0.935769511,//0.5
|
||||
1.94926E-06,//0.6
|
||||
2.22356E-05,//0.8
|
||||
3.29279E-05,//1.0
|
||||
0.000248357,//1.5
|
||||
0.000398671,//2.0
|
||||
0.000410038,//3.0
|
||||
0.001042575,//4.0
|
||||
0.000961964,//5.0
|
||||
8.714E-05,//6.0
|
||||
0.000248803,//7.0
|
||||
0.00017111,//8.0
|
||||
3.9074E-05,//9.0
|
||||
7.28E-07,//10.0
|
||||
6.27e-6};//11.0
|
||||
|
||||
/*
|
||||
//conversion test input
|
||||
double fraction [] = {0.5,0.5};
|
||||
double energy [] = {2.0,2.0};
|
||||
int nEgroup= 2;
|
||||
*/
|
||||
//*** SHIELDING DATA TABLES ****//
|
||||
// Mashkovich set + NIST data for 15 MeV
|
||||
// Linear attenuation for concrete
|
||||
double AttenuationArgsConc[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataConc[]={31.7, 28.5, 24.6, 21.9, 20.0, 18.5, 16.3, 14.6, 13.1, 11.9, 10.3, 8.74, 8.37, 7.34, 6.65, 6.19, 5.61, 5.29, 4.8208}; // in m^-1
|
||||
int AttenuationSizeConc[]={18};
|
||||
// Iron
|
||||
double AttenuationArgsFe[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataFe[]={139., 106., 83.3, 71.7, 64.6, 59.5, 52.0, 46.7, 42.2, 38.1, 33.3, 29.1, 28.4, 26.0, 24.8, 24.0, 23.4, 23.4, 24.3}; //in m^-1
|
||||
int AttenuationSizeFe[]={19};
|
||||
double AttenuationArgsPb[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataPb[]={2180., 1070., 425., 244., 170., 133., 95.2, 77.1, 65.8, 57.7, 50.8, 47.6, 46.8, 47.2, 48.1, 49.4, 52.0, 55.0, 64.2}; //in m^-1
|
||||
int AttenuationSizePb[]={19};
|
||||
/*
|
||||
//NIST set
|
||||
//Concrete
|
||||
double AttenuationArgsConc[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataConc[]={
|
||||
33.028,
|
||||
29.486,
|
||||
25.231,
|
||||
22.5009,
|
||||
20.5045,
|
||||
18.9428,
|
||||
16.6221,
|
||||
14.9385,
|
||||
13.3561,
|
||||
12.1624,
|
||||
10.4811,
|
||||
8.5123,
|
||||
7.3991,
|
||||
6.6884,
|
||||
6.2031,
|
||||
5.5936,
|
||||
5.2394,
|
||||
4.8208};
|
||||
int AttenuationSizeConc[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsFe[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataFe[]={
|
||||
1.55E+02,
|
||||
1.15E+02,
|
||||
8.65E+01,
|
||||
7.40E+01,
|
||||
6.63E+01,
|
||||
6.07E+01,
|
||||
5.27E+01,
|
||||
4.72E+01,
|
||||
4.21E+01,
|
||||
3.84E+01,
|
||||
3.36E+01,
|
||||
2.85E+01,
|
||||
2.61E+01,
|
||||
2.48E+01,
|
||||
2.41E+01,
|
||||
2.36E+01,
|
||||
2.36E+01,
|
||||
2.43E+01};
|
||||
int AttenuationSizeFe[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsPb[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in meV
|
||||
double AttenuationDataPb[]={
|
||||
2.28E+03,
|
||||
1.13E+03,
|
||||
4.57E+02,
|
||||
2.63E+02,
|
||||
1.83E+02,
|
||||
1.42E+02,
|
||||
1.01E+02,
|
||||
8.05E+01,
|
||||
6.66E+01,
|
||||
5.92E+01,
|
||||
5.22E+01,
|
||||
4.80E+01,
|
||||
4.76E+01,
|
||||
4.84E+01,
|
||||
4.98E+01,
|
||||
5.30E+01,
|
||||
5.64E+01,
|
||||
6.42E+01
|
||||
}; //in m^-1
|
||||
int AttenuationSizePb[]={18};
|
||||
*/
|
||||
|
||||
/*
|
||||
//Linear attenuation test input
|
||||
double AttenuationArgs[]={0.15, 10.0}; // in meV
|
||||
double AttenuationData[]={7.0, 7.0 }; // in m^-1
|
||||
int AttenuationSize[]={2};
|
||||
*/
|
||||
|
||||
// Dose buildup factors concrete
|
||||
double BuildupDataConc[]={1., 1.74, 2.26, 2.95, 3.79, 4.51, 5.57, 6.51, 3.18,
|
||||
1., 2.82, 5.13, 11.2, 24.2, 42.7, 87.6, 153., 353.,
|
||||
1., 2.52, 4.66, 10.8, 25.6, 48.2, 107., 198., 497.,
|
||||
1., 2.27, 4.03, 8.97, 20.2, 30.4, 75.6, 131, 292,
|
||||
1., 1.98, 3.24, 6.42, 12.7, 20.7, 37.2, 57.1, 106,
|
||||
1., 1.77, 2.65, 4.61, 7.97, 11.7, 18.6, 26.0, 42.2,
|
||||
1., 1.67, 2.38, 3.84, 6.20, 8.71, 13.1, 17.7, 27.4,
|
||||
1., 1.61, 2.18, 3.37, 5.23, 7.15, 10.5, 13.9, 20.9,
|
||||
1., 1.49, 1.93, 2.80, 4.14, 5.52, 7.86, 10.2, 15.5,
|
||||
1., 1.41, 1.76, 2.45, 3.51, 4.59, 6.43, 8.31, 12.2,
|
||||
1., 1.35, 1.64, 2.22, 3.10, 4.01, 5.57, 7.19, 10.6,
|
||||
1., 1.26, 1.46, 1.86, 2.50, 3.16, 4.34, 5.59, 8.27};
|
||||
|
||||
int BuildupSizeConc[] ={12, //energy
|
||||
9}; // mu d
|
||||
double BuildupArgsConc[]={
|
||||
0.05, 0.15, 0.3, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 // energy in MeV
|
||||
,
|
||||
0., 1., 2., 4., 7., 10., 15., 20., 30. //mu d
|
||||
};
|
||||
//Dose buildup factors Fe
|
||||
double BuildupDataFe[]={
|
||||
1., 1.5, 2.2, 3.1, 4.1, 4.6, 5.4, 5.9, //0.1
|
||||
1., 2.0, 3.1, 5.3, 8.9, 14., 22., 31., //0.2
|
||||
1., 2.1, 3.3, 6.0, 12., 23., 49., 84., //0.4
|
||||
1., 1.98, 3.09, 5.98, 11.7, 19.2, 35.4, 55.6, //0.5
|
||||
1., 1.87, 2.89, 5.39, 10.2, 16.2, 28.3, 42.7, //1.0
|
||||
1., 1.76, 2.43, 4.13, 7.25, 10.9, 17.6, 25.1, //2.0
|
||||
1., 1.55, 2.15, 3.51, 5.85, 8.51, 13.5, 19.1, //3.0
|
||||
1., 1.45, 1.94, 3.03, 4.91, 7.11, 11.2, 16.0, //4.0
|
||||
1., 1.34, 1.72, 2.58, 4.14, 6.02, 9.89, 14.7, //6.0
|
||||
1., 1.27, 1.56, 2.23, 3.49, 5.07, 8.50, 13.0, //8.0
|
||||
1., 1.20, 1.42, 1.95, 2.99, 4.35, 7.54, 12.4, //9.0
|
||||
1., 1.48, 1.86, 2.72, 4.30 , 6.37, 11.4, 19.1};//15.0
|
||||
int BuildupSizeFe[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsFe[]={0.1, 0.2, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
//Dose buildup factors Pb
|
||||
double BuildupDataPb[]={1., 1.01, 1.03, 1.06, 1.15, 1.16, 1.18, 1.19,
|
||||
1., 1.11, 1.17, 1.25, 1.34, 1.41, 1.5, 1.56,
|
||||
1., 1.17, 1.29, 1.46, 1.58, 1.72, 1.89, 2.02,
|
||||
1., 1.24, 1.42, 1.69, 2.00, 2.27, 2.65, 2.73,
|
||||
1., 1.37, 1.69, 2.26, 3.02, 3.74, 4.81, 5.86,
|
||||
1., 1.39, 1.76, 2.51, 3.66, 4.84, 6.87, 9.00,
|
||||
1., 1.34, 1.68, 2.43, 3.75, 5.30, 8.44, 12.3,
|
||||
1., 1.27, 1.56, 2.25, 3.61, 5.44, 9.80, 16.3,
|
||||
1., 1.21, 1.46, 2.08, 3.44, 5.55, 11.7, 23.6,
|
||||
1., 1.18, 1.40, 1.97, 3.34, 5.69, 13.8, 32.7,
|
||||
1., 1.14, 1.30, 1.74, 2.89, 5.07, 14.1, 44.6,
|
||||
1., 1.11, 1.23, 1.58, 2.52, 4.34, 12.5, 39.2};
|
||||
|
||||
int BuildupSizePb[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsPb[]={0.15, 0.30, 0.40, 0.5, 1.0, 2.0, 3.0, 4.0, 5.1, 6.0, 8.0, 10.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
|
||||
|
||||
// Flux-to-dose conversion
|
||||
/* //NRB-99?
|
||||
double FtoDArgs[]= {0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.80, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*0.752E-10, 3600*1.00E-10, 3600*1.51E-10, 3600*2.00E-10, 3600*2.47E-10, 3600*2.91E-10, 3600*3.73E-10, 3600*4.48E-10, 3600*7.49E-10, 3600*12.0E-10, 3600*16.0E-10, 3600*19.9E-10, 3600*23.8E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 13;
|
||||
int FtoDSize[] = {13};
|
||||
*/
|
||||
//ESS
|
||||
double FtoDArgs[]= {
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.511,
|
||||
0.6,
|
||||
0.662,
|
||||
0.8,
|
||||
1.,
|
||||
1.12,
|
||||
1.33,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
6.13,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in MeV
|
||||
double FtoDData[] = {
|
||||
2.69E-07,
|
||||
3.60E-07,
|
||||
5.44E-07,
|
||||
7.20E-07,
|
||||
8.89E-07,
|
||||
9.07E-07,
|
||||
1.05E-06,
|
||||
1.14E-06,
|
||||
1.34E-06,
|
||||
1.62E-06,
|
||||
1.76E-06,
|
||||
2.01E-06,
|
||||
2.20E-06,
|
||||
2.69E-06,
|
||||
3.51E-06,
|
||||
4.21E-06,
|
||||
4.82E-06,
|
||||
5.40E-06,
|
||||
5.47E-06,
|
||||
6.70E-06,
|
||||
7.92E-06,
|
||||
1.09E-05
|
||||
}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 22;
|
||||
int FtoDSize[] = {22};
|
||||
|
||||
|
||||
|
||||
// test input flux to dose
|
||||
/*
|
||||
double FtoDArgs[]= {0.15, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*7.0E-10,3600*7.0E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 2;
|
||||
int FtoDSize[] = {2};
|
||||
*/
|
||||
#endif //end of insertion of datatables.
|
||||
/*** end of insertion to DECLARE section **/
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
/**Reading datafiles with capture per bin**/
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
double d1, d2, d3, d4;
|
||||
int ibin=0,NBINS;
|
||||
char line[1000],line1[1000];
|
||||
char dirname[1000];
|
||||
memset(dirname,0,sizeof(dirname));
|
||||
strcat(strcat(dirname,mcdirname),separator);
|
||||
|
||||
char filename[1000];
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FILE *dataIn;
|
||||
|
||||
memset(filename,0,sizeof(filename));//resetting filename to NULL string
|
||||
strcat(strcat(filename,dirname),NiCaptureFile);
|
||||
|
||||
printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",NiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4) ibin++;
|
||||
};
|
||||
|
||||
|
||||
NBINS=ibin;
|
||||
|
||||
//printf("Shielding calculator:\n Input file contains %d bins.\n", NBINS);
|
||||
|
||||
double* zhat = malloc(NBINS*sizeof(double));
|
||||
double* Ihat= malloc(NBINS*sizeof(double));
|
||||
double* IhatNi= malloc(NBINS*sizeof(double));
|
||||
double* IhatTi= malloc(NBINS*sizeof(double));
|
||||
rewind(dataIn);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatNi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
memset(filename,0,sizeof(filename));// setting filename to NULL string
|
||||
strcat(strcat(filename,dirname),TiCaptureFile); //setting filename to TiCapture file
|
||||
printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatTi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
fclose(dataIn);
|
||||
|
||||
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),TotalCaptureFile);
|
||||
printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TotalCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
Ihat[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
int i,ilayerConc=1,ilayerFe=1,ilayerFeConc=1,ilayerPb=1,iz,iEgroup; // Auxiliary variables
|
||||
|
||||
double zbinlength = zhat[2]-zhat[1];
|
||||
|
||||
|
||||
|
||||
double* RvaluesPb = malloc(NBINS*sizeof(double));
|
||||
double* RvaluesFe= malloc(NBINS*sizeof(double));
|
||||
double* Rvaluesbt= malloc(NBINS*sizeof(double));
|
||||
double* thicknessPb= malloc(NBINS*sizeof(double));
|
||||
double* thicknessFe= malloc(NBINS*sizeof(double));
|
||||
double* thicknessbt= malloc(NBINS*sizeof(double));
|
||||
double* RvaluesFebt= malloc(NBINS*sizeof(double));
|
||||
double* thicknessFebt= malloc(NBINS*sizeof(double));
|
||||
|
||||
double* doseFe= malloc(NBINS*sizeof(double));
|
||||
double* dosebt= malloc(NBINS*sizeof(double));
|
||||
double* doseFebt= malloc(NBINS*sizeof(double));
|
||||
double* dosePb= malloc(NBINS*sizeof(double));
|
||||
|
||||
|
||||
|
||||
|
||||
/* for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
doseFe[i]=0.0; dosebt[i]=0.0;doseFebt[i]=0.0;dosePb[i]=0.0;
|
||||
}
|
||||
*/
|
||||
|
||||
for ( iz=0; iz< NBINS; iz++)
|
||||
{
|
||||
double zc=zhat[iz]; // the z-value for the calculation
|
||||
double z;
|
||||
|
||||
//CONCRETE
|
||||
// ilayerConc may have already some value calculated for an upstream piece.
|
||||
// We start from it and if shielding is too thick, reduce this number.
|
||||
do {
|
||||
double thickbt = (ilayerConc)*stepbt;
|
||||
double Rbt = (0.5*Innerspace) + thickbt;
|
||||
Rvaluesbt[iz]=Rbt;
|
||||
dosebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
dBtmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*(Rbt-(0.5*Innerspace))*sqrt(Rbt*Rbt+z*z)/Rbt;
|
||||
double args[]={En,dBtmfp};
|
||||
dosebt[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dBtmfp)/(Rbt*Rbt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerConc--;
|
||||
} while ((dosebt[iz] < MaxRate)&&(ilayerConc>0));
|
||||
|
||||
//increase ilayerConc until the dose rate requirements are satisfied.
|
||||
do {
|
||||
double thickbt = (ilayerConc)*stepbt;
|
||||
double Rbt = (0.5*Innerspace) + thickbt;
|
||||
Rvaluesbt[iz]=Rbt;
|
||||
dosebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
dBtmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*(Rbt-(0.5*Innerspace))*sqrt(Rbt*Rbt+z*z)/Rbt;
|
||||
double args[]={En,dBtmfp};
|
||||
dosebt[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dBtmfp)/(Rbt*Rbt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerConc++;
|
||||
} while ((dosebt[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
//CONCRETE+TUBING
|
||||
do{
|
||||
double thickFebt= (ilayerFeConc)*stepcomb;
|
||||
double RFebt=(0.5*Innerspace)+TFe+thickFebt;
|
||||
RvaluesFebt[iz]=RFebt;
|
||||
doseFebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
|
||||
dFebtmfp=(lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*thickFebt+
|
||||
lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TFe)
|
||||
*sqrt(RFebt*RFebt+z*z)/RFebt;
|
||||
double args[]={En,dFebtmfp};
|
||||
|
||||
doseFebt[iz]+= (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFebtmfp)/(RFebt*RFebt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFeConc--;
|
||||
} while ((doseFebt[iz] < MaxRate)&&(ilayerFeConc>0));
|
||||
|
||||
|
||||
do {
|
||||
double thickFebt= (ilayerFeConc)*stepcomb;
|
||||
double RFebt=(0.5*Innerspace)+TFe+thickFebt;
|
||||
RvaluesFebt[iz]=RFebt;
|
||||
doseFebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
|
||||
dFebtmfp=(lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*thickFebt+
|
||||
lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TFe)
|
||||
*sqrt(RFebt*RFebt+z*z)/RFebt;
|
||||
double args[]={En,dFebtmfp};
|
||||
|
||||
doseFebt[iz]+= (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFebtmfp)/(RFebt*RFebt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFeConc++;
|
||||
} while ((doseFebt[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//IRON
|
||||
do {
|
||||
double thickFe = (ilayerFe)*stepFe;
|
||||
double RFe = (0.5*Innerspace) + thickFe;
|
||||
RvaluesFe[iz]=RFe;
|
||||
doseFe[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
|
||||
dFemfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*thickFe*sqrt(RFe*RFe+z*z)/RFe;
|
||||
double args[]={En,dFemfp};
|
||||
doseFe[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFemfp)/(RFe*RFe+z*z)*
|
||||
lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFe--;
|
||||
} while ((doseFe[iz] <MaxRate)&&(ilayerFe>0));
|
||||
|
||||
|
||||
|
||||
do {
|
||||
double thickFe = (ilayerFe)*stepFe;
|
||||
double RFe = (0.5*Innerspace) + thickFe;
|
||||
RvaluesFe[iz]=RFe;
|
||||
doseFe[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
dFemfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*(RFe-(0.5*Innerspace))*sqrt(RFe*RFe+z*z)/RFe;
|
||||
double args[]={En,dFemfp};
|
||||
doseFe[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFemfp)/(RFe*RFe+z*z)*
|
||||
lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFe++;
|
||||
} while ((doseFe[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//LEAD
|
||||
do {
|
||||
double thickPb = (ilayerPb)*stepPb;
|
||||
double RPb = (0.5*Innerspace) + thickPb;
|
||||
RvaluesPb[iz]=RPb;
|
||||
dosePb[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dPbmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
dPbmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*(RPb-(0.5*Innerspace))*sqrt(RPb*RPb+z*z)/RPb;
|
||||
double args[]={En,dPbmfp};
|
||||
dosePb[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dPbmfp)/(RPb*RPb+z*z)*
|
||||
lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerPb--;
|
||||
} while ((dosePb[iz] <MaxRate)&&(ilayerPb>0));
|
||||
|
||||
|
||||
do {
|
||||
double thickPb = (ilayerPb)*stepPb;
|
||||
double RPb = (0.5*Innerspace) + thickPb;
|
||||
RvaluesPb[iz]=RPb;
|
||||
dosePb[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dPbmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
dPbmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*(RPb-(0.5*Innerspace))*sqrt(RPb*RPb+z*z)/RPb;
|
||||
double args[]={En,dPbmfp};
|
||||
dosePb[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dPbmfp)/(RPb*RPb+z*z)*
|
||||
lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerPb++;
|
||||
} while ((dosePb[iz] >=MaxRate));
|
||||
|
||||
|
||||
} /* Position scan iz */
|
||||
|
||||
/* Selection of appropriate radius for the gamma shielding */
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
int index=0;
|
||||
thicknessFe[i]=RvaluesFe[i]-(0.5*Innerspace);
|
||||
thicknessbt[i]=Rvaluesbt[i]-(0.5*Innerspace);
|
||||
thicknessPb[i]=RvaluesPb[i]-(0.5*Innerspace);
|
||||
thicknessFebt[i]=RvaluesFebt[i]-(0.5*Innerspace)-TFe;
|
||||
}
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
FILE *cost_out;
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),OutputFile);
|
||||
printf("Shielding calculator:\nWriting to file %s\n",filename);
|
||||
cost_out=fopen(filename,"w");
|
||||
|
||||
|
||||
fprintf(cost_out,"#Required shielding dimensions for the dose rate not to exceed %g uSv/hr at the surface with shielding inner dimensions %g meters\n#Position z,m\tNeutrons Ihat(z),n/s/bin\tNeutrons I(z),n/s/m\tCaptured by Ni /s/bin \tCapturedby Ni /s/m\tCaptured by Ti /s/bin\tCaptured by Ti /s/m\tThickness Pb, m\tThickness Fe, m\tThickness concrete, m (no tubing)\tThickness concrete, m (steel tubing %g cm)\tOuter radius Pb,m\tOuter radius Fe,m\tOuter radius concrete, m\tOuter radius concrete,m (with tubing)\n",MaxRate,Innerspace,100*TFe);
|
||||
|
||||
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
fprintf( cost_out,"%g\t%.4g\t%.4g\t%.4g\t%.4g\t%.4g\t%.4g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",zhat[i],Ihat[i], Ihat[i]/zbinlength, IhatNi[i], IhatNi[i]/zbinlength, IhatTi[i],IhatTi[i]/zbinlength,thicknessPb[i], thicknessFe[i],thicknessbt[i],thicknessFebt[i], RvaluesPb[i],RvaluesFe[i],Rvaluesbt[i],RvaluesFebt[i]);
|
||||
}
|
||||
fclose(cost_out);
|
||||
free(zhat); free(Ihat); free(IhatNi); free(IhatTi);
|
||||
free(RvaluesPb);
|
||||
free(RvaluesFe);
|
||||
free(Rvaluesbt);
|
||||
free(thicknessPb);
|
||||
free(thicknessFe);
|
||||
free(thicknessbt);
|
||||
free(RvaluesFebt);
|
||||
free(thicknessFebt);
|
||||
|
||||
free(doseFe);
|
||||
free(dosebt);
|
||||
free(doseFebt);
|
||||
free(dosePb);
|
||||
|
||||
printf("Done with writing shielding\n");
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,210 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Updated by: Rodon Kolevatov
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_Ni_new
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef M1THICKNESS
|
||||
#define M1THICKNESS 1500.0 //thickness of the m=1 coating
|
||||
#endif
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states
|
||||
with weight corresponding to capture in the Nickel layers*/
|
||||
int exit_neutron_Ni(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*momentum transfer at reflection. ASSUMES NO GRAVITY???*/
|
||||
double velocity=sqrt((S1->_vx)*(S1->_vx)+(S1->_vy)*(S1->_vy)+(S1->_vz)*(S1->_vz));
|
||||
double qtransf = V2Q*sqrt((S1->_vx-S0->_vx)*(S1->_vx-S0->_vx)+(S1->_vy-S0->_vy)*(S1->_vy-S0->_vy)+(S1->_vz-S0->_vz)*(S1->_vz-S0->_vz));
|
||||
double qinm = qtransf/0.0218;
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight of capture in Ni is determined analytically*/
|
||||
double captureprob;
|
||||
if ((S1->_mvalue<0)||(qinm<0.0001)) {ns_tilde[10]=0.0; captureprob=0.0;}
|
||||
else if (qinm<=1.01){
|
||||
/*q<qc(Ni), loss due to diffusion beyond the coating cutoff*/
|
||||
/* captureprob = 4.0582E-09*2200.0/velocity/((0.085787L*18.5 + 0.008321L*5.7)*1.0E-08+4.0582E-09*2200.0/velocity)*(S0->_p-S1->_p)/S0->_p; */
|
||||
double sigmadif ;//diffuse scattering cross section
|
||||
if (velocity > 1950.0){ sigmadif = 18.5; //lambda < 2A, totally incoherent scattering
|
||||
} else if (velocity < 1300.0){sigmadif = 5.2; // lambda >3A, totally coherent regime, only incoherent part contributes
|
||||
} else sigmadif=5.2+13.3*(1950-velocity)/(1950.0-1300.0);
|
||||
/*Now checking if the coating is m=1 or has higher m*/
|
||||
|
||||
if(S1->_mvalue<1.05) /*m=1 coating=> "triple thickness approximation" with Im k determined for reflectivity dip*/
|
||||
{
|
||||
double imk; /*imaginary part of momentum in the layer*/
|
||||
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
|
||||
captureprob = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
|
||||
} else { // reflection loss below qcNi in case of multilayer.
|
||||
/*Taking the minimum of what is for m=1 coating and what is when assume that neutron is physically lost beyond the cutoff.*/
|
||||
double imk; /*imaginary part of momentum in the layer*/
|
||||
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
|
||||
double c1, c2;
|
||||
c1 = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
|
||||
c2 = 0.005*(S1->_mvalue+0.1)*(S0->_p-S1->_p)/S0->_p; // for the matters of conservative estimate might be divided by 1-R at mvalue+0.1.
|
||||
captureprob = (c1>c2) ? c1 : c2 ;
|
||||
}
|
||||
} //end of "if (qinm<=1.02)"
|
||||
else if (qinm<=S1->_mvalue+0.1) { /*q>q_c(Ni) and reflection, absorption per hit */
|
||||
captureprob=0.005*qinm;
|
||||
} else { /* transmission beyond smirrorcutoff, some reserve for m=1 coating*/
|
||||
if(S1->_mvalue<1.05) captureprob=0.0025*1.3*1.3/qinm;
|
||||
else captureprob=0.0025*(S1->_mvalue+0.1)*(S1->_mvalue+0.1)/qinm;}
|
||||
/*check change in weight*/
|
||||
if((S0->_p-S1->_p)>1.e-5*S0->_p) ns_tilde[10]=S0->_p*captureprob; else ns_tilde[10]=0.0;
|
||||
/* if "mvalue" is -1, then absorption happened not on the coating but somewhere else. The weight of capture in nickel should be set to zero.*/
|
||||
|
||||
if(ns_tilde[10]>(S0->_p-S1->_p)) { printf("%f\t%f\t%f\t%E\t%E\t%f\t%E\n",velocity, qinm, S1->_mvalue, S0->_p, S1->_p, captureprob, ns_tilde[10]); exit(0);}
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function_Ni) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function_Ni=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function_Ni=exit_neutron_Ni;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
//printf("Entering Iterator_Ni\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function_Ni(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Ni: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,186 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_Ti_new
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states
|
||||
with weights corresponding to capture in Ti layers*/
|
||||
int exit_neutron_Ti(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*momentum transfer at reflection*/
|
||||
double velocity=sqrt((S1->_vx)*(S1->_vx)+(S1->_vy)*(S1->_vy)+(S1->_vz)*(S1->_vz));
|
||||
double qtransf = V2Q*sqrt((S1->_vx-S0->_vx)*(S1->_vx-S0->_vx)+(S1->_vy-S0->_vy)*(S1->_vy-S0->_vy)+(S1->_vz-S0->_vz)*(S1->_vz-S0->_vz));
|
||||
double qinm = qtransf/0.0218;
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight of capture in Ni is determined analytically*/
|
||||
double captureprob;
|
||||
/* if "mvalue" is -1, than absorption happened not on the coating but somewhere else. The weight of capture in nickel should be set to zero.*/
|
||||
if ((S1->_mvalue<0)||(qinm<0.0001)) {ns_tilde[10]=0.0; captureprob=0.0;}
|
||||
else if (S1->_mvalue<1.02){captureprob=0.0; //negligible probability for capture in Ti for m=1 coatings
|
||||
}
|
||||
else // now for coatings with m>1
|
||||
{if (qinm<=1.02){
|
||||
/*q<qc(Ni) means diffusion to higher divergences and interaction with coating there*/
|
||||
captureprob=0.0045*(S1->_mvalue-0.9)*(S0->_p-S1->_p)/S0->_p;
|
||||
} else if (qinm<=S1->_mvalue+0.1) { /*q>q_c(Ni) and reflection, absorption per hit */
|
||||
captureprob=0.0045*(qinm-1.0);
|
||||
} else { /* transmission beyond smirrorcutoff*/
|
||||
captureprob=0.00225*(S1->_mvalue-0.9)*(S1->_mvalue+0.1)/qinm;}}
|
||||
/*check change in weight*/
|
||||
if((S0->_p-S1->_p)>1.e-5*S0->_p) ns_tilde[10]=S0->_p*captureprob; else ns_tilde[10]=0.0;
|
||||
//printf("%f\t%f\t%f\t%E\t%E\t%f\t%E\n",velocity, qinm, S1->_mvalue, S0->_p, S1->_p, captureprob, ns_tilde[10]);
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function_Ti) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function_Ti=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function_Ti=exit_neutron_Ti;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
//printf("Entering Iterator_Ti\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function_Ti(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Ti: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,116 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator_stop.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration stop element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the end of the trace-region in which pseudo-neutrons are handled. Please see the Scatter_log_iterator-component for more details.
|
||||
* N.B. This component should be immediately followed by a construction like:
|
||||
* COMPONENT a1 = Arm()
|
||||
* AT (0,0,0) ABSOLUTE
|
||||
* JUMP a0 WHEN(MC_GETPAR(iterator_stop_comp_name,loop))
|
||||
*
|
||||
* This is to extract the value of the loop variable from the innards of this component
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* iterator [ ] Instance name of the Scatter_log_iterator log component preceeding this one.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_stop
|
||||
DEFINITION PARAMETERS (iterator)
|
||||
SETTING PARAMETERS (int last=0)
|
||||
OUTPUT PARAMETERS (loop)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int loop;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
loop=1;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iterator
|
||||
scatter_iterator_stop:
|
||||
loop=1;
|
||||
struct Generalized_State_t *s1=MC_GETPAR(iterator,s1);
|
||||
|
||||
if (s1->_p==-1){
|
||||
/*we have reached the end - unset loop and reset neutron state to whatever it was before we entered the pseudo neutron iterator*/
|
||||
loop=0;
|
||||
double *ns=MC_GETPAR(iterator,nstate_initial);
|
||||
x=ns[0];y=ns[1];z=ns[2];
|
||||
vx=ns[3];vy=ns[4];vz=ns[5];
|
||||
t=ns[6];
|
||||
sx=ns[7];sy=ns[8];sz=ns[9];
|
||||
p=ns[10];
|
||||
|
||||
free(ns);
|
||||
MC_GETPAR(iterator,nstate_initial)=NULL;
|
||||
/* Restore std ABSORB */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
// printf("No ABSORB after the iterator\n");
|
||||
} else if (s1->_p==-2) {
|
||||
/*if the weight equals -2 it means that the neutron was absorbed in components within the scatter_logger (see scatter logger definition)
|
||||
and should not be propagated further after this iterator*/
|
||||
// printf("Performing ABSORB after the iterator\n");
|
||||
/*we have reached the end - unset loop and reset neutron state to whatever it was before we entered the pseudo neutron iterator*/
|
||||
loop=0;
|
||||
double *ns=MC_GETPAR(iterator,nstate_initial);
|
||||
x=ns[0];y=ns[1];z=ns[2];
|
||||
vx=ns[3];vy=ns[4];vz=ns[5];
|
||||
t=ns[6];
|
||||
sx=ns[7];sy=ns[8];sz=ns[9];
|
||||
p=ns[10];
|
||||
|
||||
free(ns);
|
||||
MC_GETPAR(iterator,nstate_initial)=NULL;
|
||||
/* Restore std ABSORB */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
if (last>0) ABSORB; //Will sample new neutron.
|
||||
}
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,165 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_total
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states*/
|
||||
int exit_neutron(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight is difference old-new to mean the neutrons "deposited" in the guide wall*/
|
||||
ns_tilde[10]=S0->_p-S1->_p;
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function=exit_neutron;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
//printf("Entering Iterator_total\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Tot: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,217 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_logger.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen, Peter K Willendrup & Esben Klinkby
|
||||
* Date: January 2013
|
||||
* Version: $Revision: 1.12 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics / DTU Nutech
|
||||
*
|
||||
* Logging iteractions of neutrons with components
|
||||
*
|
||||
* %D
|
||||
* Start of the trace-region in which SCATTER events should be logged.
|
||||
* Whenever a SCATTER occurs in components between this one and its counterpart Scatter_logger_stop the neutron state is logged.
|
||||
* The log is kept in memory - but only for one single
|
||||
* numerical neutron at a time, so there should be no real danger of memory overrun.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* (none)
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_logger
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
|
||||
// m_local_refl is to be used in shielding applications and accessed by shielding logger
|
||||
#ifndef MVALUE_LOCAL_IS_DEF
|
||||
#define MVALUE_LOCAL_IS_DEF 1
|
||||
double m_local_refl=-1;
|
||||
#endif
|
||||
|
||||
struct Generalized_State_t {
|
||||
double _x,_y,_z,_vx,_vy,_vz;
|
||||
double _p,_t,_sx,_sy,_sz;
|
||||
double _mvalue; // coating m-value in the place of the hit
|
||||
long long int _nid;
|
||||
int _comp, _idx;
|
||||
};
|
||||
|
||||
#define SCATTER_LOG do { \
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE-1){\
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);\
|
||||
if( (bp-1)->_p!=p || (bp-1)->_vx!=vx || (bp-1)->_vy!=vy || (bp-1)->_vz!=vz ){\
|
||||
Coords ctmp=POS_A_CURRENT_COMP;\
|
||||
Coords _r = coords_set(x,y,z);\
|
||||
Coords _v = coords_set(vx,vy,vz);\
|
||||
Coords _s = coords_set(sx,sy,sz);\
|
||||
Coords _rg,_vg,_sg;\
|
||||
Rotation _Rt;\
|
||||
rot_transpose(ROT_A_CURRENT_COMP,_Rt);\
|
||||
_rg=coords_add(rot_apply(_Rt,_r),ctmp);\
|
||||
_vg=rot_apply(_Rt,_v);\
|
||||
_sg=rot_apply(_Rt,_s);\
|
||||
coords_get(_rg,&(bp->_x),&(bp->_y),&(bp->_z));\
|
||||
coords_get(_vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));\
|
||||
coords_get(_sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));\
|
||||
bp->_t=t;\
|
||||
bp->_p=p;\
|
||||
bp->_nid=mcget_run_num();\
|
||||
bp->_comp=INDEX_CURRENT_COMP;\
|
||||
bp->_idx=bounce_store_index;\
|
||||
/* New: registering m-value at reflection. it is set to -1 by component the if we missed coating */\
|
||||
bp->_mvalue=m_local_refl;\
|
||||
/* printf("Recording scattering, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p:%g\n",\
|
||||
bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);*/\
|
||||
bounce_store_index++;\
|
||||
}\
|
||||
}else if(bounce_store_index==(BOUNCE_LOG_SIZE-1) && !bounce_store_overrun){\
|
||||
printf("Warning (%s): Scatter_log overrun at %llu - not logging any more events\n",NAME_CURRENT_COMP,mcget_run_num());\
|
||||
bounce_store_overrun=1;\
|
||||
}\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0);\
|
||||
} while(0)
|
||||
|
||||
#define ABSORB_LOG do { /*printf("DOING ABSORB_LOG\n");*/ \
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE-1){\
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);\
|
||||
/* if( (bp-1)->_p!=p || (bp-1)->_vx!=vx || (bp-1)->_vy!=vy || (bp-1)->_vz!=vz )*//*<-- Check removed, works wrong if neutron is absorbed at the entrance*/{\
|
||||
Coords ctmp=POS_A_CURRENT_COMP;\
|
||||
Coords _r = coords_set(x,y,z);\
|
||||
Coords _v = coords_set(vx,vy,vz);\
|
||||
Coords _s = coords_set(sx,sy,sz);\
|
||||
Coords _rg,_vg,_sg;\
|
||||
Rotation _Rt;\
|
||||
rot_transpose(ROT_A_CURRENT_COMP,_Rt);\
|
||||
_rg=coords_add(rot_apply(_Rt,_r),ctmp);\
|
||||
_vg=rot_apply(_Rt,_v);\
|
||||
_sg=rot_apply(_Rt,_s);\
|
||||
coords_get(_rg,&(bp->_x),&(bp->_y),&(bp->_z));\
|
||||
coords_get(_vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));\
|
||||
/*bp->_vx=0.; bp->_vy=0.; bp->_vz=0.;*/\
|
||||
/*vx=0; vy=0;vz=0;*/\
|
||||
coords_get(_sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));\
|
||||
bp->_t=t;\
|
||||
bp->_p=0.;\
|
||||
bp->_nid=mcget_run_num();\
|
||||
bp->_comp=INDEX_CURRENT_COMP;\
|
||||
bp->_idx=bounce_store_index;\
|
||||
bp->_mvalue=m_local_refl;\
|
||||
/* printf("Recording absorption event, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p:%g\n",\
|
||||
bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);*/\
|
||||
bounce_store_index++;\
|
||||
}\
|
||||
}else if(bounce_store_index==(BOUNCE_LOG_SIZE-1) && !bounce_store_overrun){\
|
||||
printf("Warning (%s): Scatter_log overrun at %llu - not logging any more events\n",NAME_CURRENT_COMP,mcget_run_num());\
|
||||
bounce_store_overrun=1;\
|
||||
}\
|
||||
absorbed_in_optics=1;\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0);\
|
||||
} while(0)
|
||||
|
||||
|
||||
#define SCATTER0\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0)
|
||||
|
||||
#define ABSORB0\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0)
|
||||
|
||||
|
||||
|
||||
const int BOUNCE_LOG_SIZE=10000000;
|
||||
struct Generalized_State_t *Bounce_store;
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int bounce_store_index;
|
||||
int bounce_store_overrun;
|
||||
int absorbed_in_optics;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
bounce_store_index=0;
|
||||
absorbed_in_optics=0;
|
||||
Bounce_store=malloc(sizeof(*Bounce_store)*BOUNCE_LOG_SIZE);
|
||||
Bounce_store[BOUNCE_LOG_SIZE-1]._p=-1;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
#undef SCATTER
|
||||
#define SCATTER SCATTER_LOG
|
||||
|
||||
#undef ABSORB
|
||||
#define ABSORB ABSORB_LOG
|
||||
|
||||
/*
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
|
||||
#define scatter_logger_stop logger
|
||||
*/
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_logger_stop
|
||||
|
||||
bounce_store_index=0;/*we are now starting logging so we should start afresh*/
|
||||
absorbed_in_optics=0;/*we are now starting logging so we should start afresh*/
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE){
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);
|
||||
Coords ctmp=POS_A_CURRENT_COMP;
|
||||
Coords r = coords_set(x,y,z);
|
||||
Coords v = coords_set(vx,vy,vz);
|
||||
Coords s = coords_set(sx,sy,sz);
|
||||
|
||||
Coords rg,vg,sg;
|
||||
|
||||
Rotation Rt;
|
||||
rot_transpose(ROT_A_CURRENT_COMP,Rt);
|
||||
rg=coords_add(rot_apply(Rt,r),ctmp);
|
||||
vg=rot_apply(Rt,v);
|
||||
sg=rot_apply(Rt,s);
|
||||
coords_get(rg,&(bp->_x),&(bp->_y),&(bp->_z));
|
||||
coords_get(vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));
|
||||
coords_get(sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));
|
||||
bp->_t=t;
|
||||
bp->_p=p;
|
||||
bp->_nid=mcget_run_num();
|
||||
bp->_comp=INDEX_CURRENT_COMP;
|
||||
bp->_idx=bounce_store_index;
|
||||
bp->_mvalue=0.0; // This is the incoming particle. No collision, setting mvalue=0.
|
||||
// printf("Starting scatter logger, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p: %g\n", bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);
|
||||
bounce_store_index++;
|
||||
}else if(bounce_store_index==BOUNCE_LOG_SIZE && !bounce_store_overrun){
|
||||
printf("Warning (%s): Scatter_log overrun - not logging any more SCATTER events\n",NAME_CURRENT_COMP);
|
||||
bounce_store_overrun=1;
|
||||
}
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,121 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_logger_stop.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen, Peter K Willendrup & Esben Klinkby
|
||||
* Date: January 2013
|
||||
* Version: $Revision: 1.12 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics / DTU Nutech
|
||||
*
|
||||
* Stop logging iteractions of neutrons with components
|
||||
*
|
||||
* %D
|
||||
* Component which marks the end of the region where SCATTER events should be logged.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* logger: The Scatter_logger.comp which began the logging region. This is necessary to allow communication between the components. ( )
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_logger_stop
|
||||
DEFINITION PARAMETERS (logger)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#define SCATTER0\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0)
|
||||
#define ABSORB0\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0)
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int bounce_store_index;
|
||||
int bounce_store_overrun;
|
||||
int logger_buffer_cleared;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
#ifndef logger
|
||||
fprintf(stderr,"Error(%s): Logger undefined - can't stop noexisting logger\n", NAME_CURRENT_COMP);
|
||||
#endif
|
||||
logger_buffer_cleared=0;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
#undef SCATTER
|
||||
#define SCATTER SCATTER0
|
||||
|
||||
#undef ABSORB
|
||||
#define ABSORB ABSORB0
|
||||
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
|
||||
#define scatter_logger_stop logger
|
||||
scatter_logger_stop:
|
||||
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE){
|
||||
struct Generalized_State_t *bp;
|
||||
bp=&(Bounce_store[bounce_store_index]);
|
||||
bp->_x=0;
|
||||
bp->_y=0;
|
||||
bp->_z=0;
|
||||
bp->_vx=0;
|
||||
bp->_vy=0;
|
||||
bp->_vz=0;
|
||||
bp->_sx=0;
|
||||
bp->_sy=0;
|
||||
bp->_sz=0;
|
||||
bp->_t=0;
|
||||
if(absorbed_in_optics==1) bp->_p=-2; else bp->_p=-1;
|
||||
bp->_nid=0;
|
||||
bp->_comp=0;
|
||||
bounce_store_index++;
|
||||
}else if(bounce_store_index==BOUNCE_LOG_SIZE && !bounce_store_overrun){
|
||||
printf("Warning (%s): Scatter_log overrun - cannot set stop bit. Aborting\n",NAME_CURRENT_COMP);
|
||||
exit(1);
|
||||
}
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
/*If we are a logger stop we should also free the memory*/
|
||||
int i;
|
||||
struct Generalized_State_t *bp;
|
||||
for (i=0;i<BOUNCE_LOG_SIZE;i++){
|
||||
bp=&(Bounce_store[i]);
|
||||
/* maybe this should be in the share block - must be able to discern between different loggers*/
|
||||
/* printf("SCATTERLOG: %d %lld %g %g %g %g %g %g %g %g %g %g %g %d\n", \*/
|
||||
/* i,bp->_nid,bp->_x,bp->_y,bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_t, \*/
|
||||
/* bp->_sx, bp->_sy, bp->_sz, bp->_p, bp->_comp);*/
|
||||
}
|
||||
if (!logger_buffer_cleared) { free(Bounce_store); logger_buffer_cleared=1;}
|
||||
|
||||
/*do nothing if the memory has already been freed*/
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -1,6 +1,4 @@
|
||||
#!/bin/bash
|
||||
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/afs/psi.ch/project/sinq/sl6-64/mcstas2.4/mcstas/2.4/libs
|
||||
|
||||
|
||||
if [ Estia_baseline.instr -nt Estia_baseline.out ] || [ ! -f Estia_baseline.out ] \
|
||||
|| [ Estia_feeder.instr -nt Estia_baseline.out ] \
|
||||
@@ -8,6 +6,6 @@ if [ Estia_baseline.instr -nt Estia_baseline.out ] || [ ! -f Estia_baseline.out
|
||||
|| [ Estia_mf.instr -nt Estia_baseline.out ] \
|
||||
|| [ Estia_selene2.instr -nt Estia_baseline.out ]; then
|
||||
rm Estia_baseline.c Estia_baseline.out
|
||||
mcstas -o Estia_baseline.c Estia_baseline.instr --trace
|
||||
mpicc -O2 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
|
||||
mcstas -o Estia_baseline.c Estia_baseline.instr
|
||||
mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI
|
||||
fi
|
||||
|
||||
@@ -1,3 +0,0 @@
|
||||
mcc05:72
|
||||
mcc06:72
|
||||
|
||||
@@ -1,4 +0,0 @@
|
||||
foreach z ( `seq -0.5 0.125 0.5` )
|
||||
echo $z
|
||||
./run_simu.sh $z
|
||||
end
|
||||
Binary file not shown.
Executable
+52
@@ -0,0 +1,52 @@
|
||||
#!/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
|
||||
@@ -0,0 +1,31 @@
|
||||
#!/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 \
|
||||
--dir=../results/shielding --ncount=1e9 \
|
||||
sample=0
|
||||
|
||||
|
||||
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