950 lines
31 KiB
Plaintext
950 lines
31 KiB
Plaintext
/*******************************************************************************
|
|
*
|
|
* 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
|