Commit 0e30b1b0 authored by jannick.wolters's avatar jannick.wolters
Browse files

added ICRU77 fortran code converted to c++

parent 5db67d32
Pipeline #110468 failed with stages
in 28 minutes and 52 seconds
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
**** Material composition, oscillator data and
stopping powers for electrons and positrons.
Material # 278, WATER, LIQUID
Mass density = 1.00000000E+00 g/cm**3
Number of elements in the molecule = 2
atomic number = 1, atoms/molecule = 2.00000000E+00
atomic number = 8, atoms/molecule = 1.00000000E+00
Mean excitation energy = 7.50000000E+01 eV
Number of oscillators = 3
Fi Ui (eV) Wi (eV)
---------------------------------------------------
1 6.00000000E+00 1.36133333E+01 3.13708109E+01
2 2.00000000E+00 2.84800000E+01 6.40292871E+01
3 2.00000000E+00 5.38000000E+02 1.20046753E+03
Stopping powers. mtu=g/cm**2
Energy Scol,e- Srad,e- Scol,e+ Srad,e+
(eV) (MeV/mtu) (MeV/mtu) (MeV/mtu) (MeV/mtu)
------------------------------------------------------------
5.00E+01 3.70852E+02 1.36585E-03 4.92645E+02 1.73951E-04
6.00E+01 3.69967E+02 1.42600E-03 5.11158E+02 1.95865E-04
7.00E+01 4.13339E+02 1.47902E-03 5.49996E+02 2.16318E-04
8.00E+01 4.27409E+02 1.52679E-03 5.58552E+02 2.35605E-04
9.00E+01 4.26017E+02 1.57017E-03 5.52407E+02 2.53874E-04
1.00E+02 4.18657E+02 1.61021E-03 5.40424E+02 2.71305E-04
1.25E+02 3.92215E+02 1.69862E-03 5.02921E+02 3.11773E-04
1.50E+02 3.68877E+02 1.77471E-03 4.65656E+02 3.48695E-04
1.75E+02 3.46353E+02 1.84203E-03 4.32255E+02 3.82856E-04
2.00E+02 3.25742E+02 1.90258E-03 4.03025E+02 4.14753E-04
2.50E+02 2.90852E+02 2.00847E-03 3.55191E+02 4.73136E-04
3.00E+02 2.63017E+02 2.09988E-03 3.18078E+02 5.25980E-04
3.50E+02 2.40447E+02 2.18058E-03 2.88540E+02 5.74506E-04
4.00E+02 2.21797E+02 2.25339E-03 2.64472E+02 6.19653E-04
4.50E+02 2.06121E+02 2.31960E-03 2.44465E+02 6.61928E-04
5.00E+02 1.92746E+02 2.38084E-03 2.27551E+02 7.01905E-04
6.00E+02 1.71091E+02 2.49090E-03 2.00463E+02 7.76038E-04
7.00E+02 1.54264E+02 2.58825E-03 1.79664E+02 8.43986E-04
8.00E+02 1.40771E+02 2.67599E-03 1.63141E+02 9.07065E-04
9.00E+02 1.29681E+02 2.75615E-03 1.49665E+02 9.66183E-04
1.00E+03 1.20387E+02 2.83015E-03 1.38442E+02 1.02199E-03
1.25E+03 1.03967E+02 2.99172E-03 1.18648E+02 1.14897E-03
1.50E+03 9.25403E+01 3.12146E-03 1.05352E+02 1.25919E-03
1.75E+03 8.31983E+01 3.22408E-03 9.45899E+01 1.35470E-03
2.00E+03 7.56610E+01 3.30698E-03 8.59347E+01 1.43866E-03
2.50E+03 6.43691E+01 3.43525E-03 7.29213E+01 1.58193E-03
3.00E+03 5.64172E+01 3.52937E-03 6.35888E+01 1.70068E-03
3.50E+03 5.03563E+01 3.60082E-03 5.65473E+01 1.80148E-03
4.00E+03 4.55809E+01 3.65660E-03 5.10302E+01 1.88866E-03
4.50E+03 4.17146E+01 3.70083E-03 4.65809E+01 1.96505E-03
5.00E+03 3.85149E+01 3.73660E-03 4.29102E+01 2.03283E-03
6.00E+03 3.35140E+01 3.79027E-03 3.71941E+01 2.14847E-03
7.00E+03 2.97717E+01 3.82778E-03 3.29339E+01 2.24410E-03
8.00E+03 2.68565E+01 3.85514E-03 2.96262E+01 2.32525E-03
9.00E+03 2.45157E+01 3.87499E-03 2.69777E+01 2.39498E-03
1.00E+04 2.25910E+01 3.89031E-03 2.48054E+01 2.45624E-03
1.25E+04 1.89924E+01 3.91535E-03 2.07574E+01 2.58175E-03
1.50E+04 1.64805E+01 3.92847E-03 1.79434E+01 2.67913E-03
1.75E+04 1.46201E+01 3.93520E-03 1.58660E+01 2.75763E-03
2.00E+04 1.31828E+01 3.94039E-03 1.42652E+01 2.82421E-03
2.50E+04 1.10997E+01 3.95251E-03 1.19524E+01 2.93548E-03
3.00E+04 9.65695E+00 3.96525E-03 1.03559E+01 3.02585E-03
3.50E+04 8.59518E+00 3.97694E-03 9.18408E+00 3.10082E-03
4.00E+04 7.77929E+00 3.98823E-03 8.28550E+00 3.16496E-03
4.50E+04 7.13175E+00 3.99954E-03 7.57359E+00 3.22120E-03
5.00E+04 6.60477E+00 4.01265E-03 6.99509E+00 3.27279E-03
6.00E+04 5.79786E+00 4.04514E-03 6.11088E+00 3.36762E-03
7.00E+04 5.20822E+00 4.08256E-03 5.46605E+00 3.45378E-03
8.00E+04 4.75809E+00 4.12356E-03 4.97457E+00 3.53407E-03
9.00E+04 4.40304E+00 4.16748E-03 4.58741E+00 3.61040E-03
1.00E+05 4.11580E+00 4.21392E-03 4.27453E+00 3.68405E-03
1.25E+05 3.59175E+00 4.33880E-03 3.70450E+00 3.86111E-03
1.50E+05 3.23794E+00 4.47486E-03 3.32024E+00 4.03432E-03
1.75E+05 2.98383E+00 4.62060E-03 3.04455E+00 4.20769E-03
2.00E+05 2.79321E+00 4.77605E-03 2.83788E+00 4.38423E-03
2.50E+05 2.52822E+00 5.11547E-03 2.55070E+00 4.75293E-03
3.00E+05 2.35500E+00 5.48815E-03 2.36296E+00 5.14416E-03
3.50E+05 2.23483E+00 5.89073E-03 2.23263E+00 5.55872E-03
4.00E+05 2.14799E+00 6.31882E-03 2.13830E+00 5.99460E-03
4.50E+05 2.08317E+00 6.76912E-03 2.06777E+00 6.44990E-03
5.00E+05 2.03313E+00 7.23911E-03 2.01324E+00 6.92297E-03
6.00E+05 1.96227E+00 8.23473E-03 1.93581E+00 7.92086E-03
7.00E+05 1.91608E+00 9.29369E-03 1.88507E+00 8.97908E-03
8.00E+05 1.88491E+00 1.04109E-02 1.85058E+00 1.00938E-02
9.00E+05 1.86341E+00 1.15735E-02 1.82656E+00 1.12532E-02
1.00E+06 1.84842E+00 1.27827E-02 1.80960E+00 1.24587E-02
1.25E+06 1.82807E+00 1.59725E-02 1.78578E+00 1.56392E-02
1.50E+06 1.82099E+00 1.93754E-02 1.77646E+00 1.90338E-02
1.75E+06 1.82040E+00 2.29647E-02 1.77430E+00 2.26164E-02
2.00E+06 1.82318E+00 2.67121E-02 1.77591E+00 2.63589E-02
2.50E+06 1.83327E+00 3.45981E-02 1.78439E+00 3.42405E-02
3.00E+06 1.84547E+00 4.29315E-02 1.79550E+00 4.25753E-02
3.50E+06 1.85788E+00 5.16306E-02 1.80712E+00 5.12801E-02
4.00E+06 1.86983E+00 6.06341E-02 1.81849E+00 6.02927E-02
4.50E+06 1.88112E+00 6.98826E-02 1.82931E+00 6.95527E-02
5.00E+06 1.89170E+00 7.93339E-02 1.83952E+00 7.90171E-02
6.00E+06 1.91083E+00 9.87446E-02 1.85807E+00 9.84566E-02
7.00E+06 1.92761E+00 1.18788E-01 1.87444E+00 1.18529E-01
8.00E+06 1.94245E+00 1.39368E-01 1.88896E+00 1.39138E-01
9.00E+06 1.95569E+00 1.60384E-01 1.90196E+00 1.60181E-01
1.00E+07 1.96762E+00 1.81765E-01 1.91369E+00 1.81586E-01
1.25E+07 1.99296E+00 2.36496E-01 1.93866E+00 2.36367E-01
1.50E+07 2.01357E+00 2.92656E-01 1.95902E+00 2.92565E-01
1.75E+07 2.03082E+00 3.49939E-01 1.97609E+00 3.49874E-01
2.00E+07 2.04556E+00 4.08148E-01 1.99069E+00 4.08103E-01
2.50E+07 2.06964E+00 5.26832E-01 2.01458E+00 5.26809E-01
3.00E+07 2.08871E+00 6.47665E-01 2.03353E+00 6.47653E-01
3.50E+07 2.10439E+00 7.70024E-01 2.04911E+00 7.70018E-01
4.00E+07 2.11764E+00 8.93826E-01 2.06229E+00 8.93823E-01
4.50E+07 2.12909E+00 1.01898E+00 2.07368E+00 1.01898E+00
5.00E+07 2.13915E+00 1.14502E+00 2.08370E+00 1.14502E+00
6.00E+07 2.15621E+00 1.39879E+00 2.10069E+00 1.39879E+00
7.00E+07 2.17034E+00 1.65536E+00 2.11478E+00 1.65536E+00
8.00E+07 2.18241E+00 1.91442E+00 2.12681E+00 1.91442E+00
9.00E+07 2.19294E+00 2.17525E+00 2.13731E+00 2.17525E+00
1.00E+08 2.20229E+00 2.43750E+00 2.14663E+00 2.43750E+00
1.25E+08 2.22190E+00 3.09741E+00 2.16620E+00 3.09741E+00
1.50E+08 2.23778E+00 3.76212E+00 2.18206E+00 3.76212E+00
1.75E+08 2.25114E+00 4.43048E+00 2.19540E+00 4.43048E+00
2.00E+08 2.26267E+00 5.10164E+00 2.20691E+00 5.10164E+00
2.50E+08 2.28188E+00 6.44895E+00 2.22610E+00 6.44895E+00
3.00E+08 2.29753E+00 7.80056E+00 2.24174E+00 7.80056E+00
3.50E+08 2.31074E+00 9.15441E+00 2.25494E+00 9.15441E+00
4.00E+08 2.32217E+00 1.05108E+01 2.26636E+00 1.05108E+01
4.50E+08 2.33224E+00 1.18697E+01 2.27643E+00 1.18697E+01
5.00E+08 2.34125E+00 1.32305E+01 2.28544E+00 1.32305E+01
6.00E+08 2.35683E+00 1.59560E+01 2.30101E+00 1.59560E+01
7.00E+08 2.37000E+00 1.86873E+01 2.31417E+00 1.86873E+01
8.00E+08 2.38140E+00 2.14227E+01 2.32556E+00 2.14227E+01
9.00E+08 2.39145E+00 2.41602E+01 2.33561E+00 2.41602E+01
1.00E+09 2.40044E+00 2.68990E+01 2.34460E+00 2.68990E+01
#ifndef ICRU_H
#define ICRU_H
#include <cmath>
#include <fstream>
#include <sstream>
#include <vector>
#include "common/typedef.h"
#include "interpolation.h"
#include "toolboxes/errormessages.h"
enum ParticleType { ELECTRON, POSITRON };
class ICRU
{
private:
const unsigned materialID = 278;
const ParticleType particle = ParticleType::ELECTRON;
const double PI = 3.1415926535897932e0;
const double R4PI = 1.0e0 / ( 4.0e0 * PI );
// const double FOURPI = 4.0e0 * PI;
// const double AVOG = 6.02214199e23;
const double EMC2 = 510.998902e3;
const double A0B = 5.291772083e-9;
const double HREV = 27.2113834e0;
const double TEMC2 = EMC2 + EMC2;
const double CONS = 2.0e0 * PI * ( HREV * A0B ) * ( HREV * A0B ) / EMC2;
const std::vector<double> EGRD{
1.00e0, 1.25e0, 1.50e0, 1.75e0, 2.00e0, 2.50e0, 3.00e0, 3.50e0, 4.00e0, 4.50e0, 5.00e0, 6.00e0, 7.00e0, 8.00e0, 9.00e0, 1.00e1 };
const std::vector<double> ELAW{
1.007900e0, 4.002600e0, 6.941000e0, 9.012200e0, 1.081100e1, 1.201070e1, 1.400670e1, 1.599940e1, 1.899840e1, 2.017970e1, 2.298980e1,
2.430500e1, 2.698150e1, 2.808550e1, 3.097380e1, 3.206600e1, 3.545270e1, 3.994800e1, 3.909830e1, 4.007800e1, 4.495590e1, 4.786700e1,
5.094150e1, 5.199610e1, 5.493800e1, 5.584500e1, 5.893320e1, 5.869340e1, 6.354600e1, 6.539000e1, 6.972300e1, 7.261000e1, 7.492160e1,
7.896000e1, 7.990400e1, 8.380000e1, 8.546780e1, 8.762000e1, 8.890590e1, 9.122400e1, 9.290640e1, 9.594000e1, 9.890630e1, 1.010700e2,
1.029055e2, 1.064200e2, 1.078682e2, 1.124110e2, 1.148180e2, 1.187100e2, 1.217600e2, 1.276000e2, 1.269045e2, 1.312900e2, 1.329054e2,
1.373270e2, 1.389055e2, 1.401160e2, 1.409076e2, 1.442400e2, 1.449127e2, 1.503600e2, 1.519640e2, 1.572500e2, 1.589253e2, 1.625000e2,
1.649303e2, 1.672600e2, 1.689342e2, 1.730400e2, 1.749670e2, 1.784900e2, 1.809479e2, 1.838400e2, 1.862070e2, 1.902300e2, 1.922170e2,
1.950780e2, 1.969666e2, 2.005900e2, 2.043833e2, 2.072000e2, 2.089804e2, 2.089824e2, 2.099871e2, 2.220176e2, 2.230197e2, 2.260254e2,
2.270277e2, 2.320381e2, 2.310359e2, 2.380289e2, 2.370482e2, 2.440642e2, 2.430614e2, 2.470000e2, 2.470000e2, 2.510000e2, 2.520000e2,
2.570000e2, 2.580000e2, 2.590000e2, 2.620000e2 };
unsigned _NE;
unsigned _NA;
unsigned _NELEM;
unsigned _NOSC;
double _CSI;
double _TCS1I, _TCS2I;
double _RMU1, _RMU2;
double _B1, _B2, _B3, _B4;
double _C0, _C1, _C2, _C3, _C4;
double _BETA2;
double _R1;
Vector _E, _QMU;
std::vector<double> _ET, _ETL;
std::vector<double> _XMU, _XMUL;
std::vector<double> _ECS, _PCS;
std::vector<double> _ETCS1, _ETCS2;
std::vector<double> _PTCS1, _PTCS2;
std::vector<double> _DCSI, _DCSIL;
std::vector<double> _STF;
std::vector<double> _IZ;
std::vector<double> _F, _WRI;
Matrix _EDCS, _PDCS;
void ELINIT( std::vector<double> IZ, std::vector<double> STF );
void DCSEL0( double E );
void SPLINE( const std::vector<double>& X,
const std::vector<double>& Y,
std::vector<double>& A,
std::vector<double>& B,
std::vector<double>& C,
std::vector<double>& D,
double S1,
double SN,
unsigned N );
void FINDI( std::vector<double> X, double XC, unsigned N, unsigned& I );
double DCSIN3( double E, double RMU );
double DCSEL( double RMU );
double DODCSD( double RMU );
double DODCSC( double RMU, double EL );
void DOICSS( double WR, double E, double RMU, double& DCSD, double& DCSC );
void angdcs( unsigned ncor, double e, std::vector<double>& dxse, std::vector<double>& dxsi, std::vector<double>& dxs );
ICRU() = delete;
public:
ICRU( const Vector& mu, const Vector& energy ) : _NE( 96 ), _NA( 606 ), _E( 1e6 * energy ), _QMU( mu ) {
std::vector<double> TH( _NA ), THR( _NA );
_XMU.resize( _NA );
for( unsigned i = 0; i < _NA; ++i ) {
if( i == 0 ) {
TH[i] = 0.0e0;
THR[i] = TH[i] * M_PI / 180.0e0;
_XMU[i] = ( 1.0e0 - std::cos( THR[i] ) ) / 2.0e0;
}
else if( i == 1 ) {
TH[i] = 1.0e-4;
THR[i] = TH[i] * M_PI / 180.0e0;
_XMU[i] = ( 1.0e0 - std::cos( THR[i] ) ) / 2.0e0;
}
else {
if( TH[i - 1] < 0.9999e-3 )
TH[i] = TH[i - 1] + 2.5e-5;
else if( TH[i - 1] < 0.9999e-2 )
TH[i] = TH[i - 1] + 2.5e-4;
else if( TH[i - 1] < 0.9999e-1 )
TH[i] = TH[i - 1] + 2.5e-3;
else if( TH[i - 1] < 0.9999e+0 )
TH[i] = TH[i - 1] + 2.5e-2;
else if( TH[i - 1] < 0.9999e+1 )
TH[i] = TH[i - 1] + 1.0e-1;
else if( TH[i - 1] < 2.4999e+1 )
TH[i] = TH[i - 1] + 2.5e-1;
else
TH[i] = TH[i - 1] + 5.0e-1;
THR[i] = TH[i] * M_PI / 180.0e0;
_XMU[i] = std::max( ( 1.0e0 - std::cos( THR[i] ) ) / 2.0e0, 1.0e-35 );
}
}
_ET.resize( _NE );
_ETL.resize( _NE );
_XMUL.resize( _NA );
_ECS.resize( _NE );
_ETCS1.resize( _NE );
_ETCS2.resize( _NE );
_PCS.resize( _NE );
_PTCS1.resize( _NE );
_PTCS2.resize( _NE );
_DCSI.resize( _NA );
_DCSIL.resize( _NA );
_EDCS.resize( _NE, _NA );
_PDCS.resize( _NE, _NA );
}
~ICRU(){};
void GetAngularScatteringXS( Matrix& elastic, Matrix& inelastic, Matrix& total );
};
#endif // ICRU_H
#include "icru.h"
void ICRU::angdcs( unsigned ncor, double e, std::vector<double>& dxse, std::vector<double>& dxsi, std::vector<double>& dxs ) {
unsigned NM = materialID;
if( NM < 1 || NM > 280 ) {
ErrorMessages::Error( "The material number " + std::to_string( NM ) + " is not defined.", CURRENT_FUNCTION );
}
std::string PATHMT = "../data/material/";
std::stringstream ss;
ss << std::setw( 3 ) << std::setfill( '0' ) << materialID;
std::string FILE1 = PATHMT + "mater" + ss.str() + ".stp";
std::ifstream f( FILE1 );
std::string line;
for( unsigned i = 0; i < 5; ++i ) std::getline( f, line );
// double RHO = std::stod( line.substr( line.find_first_of( "=" ) + 1, 16 ) );
std::getline( f, line );
_NELEM = std::stoul( line.substr( line.find_first_of( "=" ) + 1, 3 ) );
double ATW = 0.0, ZTOT = 0.0;
_IZ.resize( _NELEM );
_STF.resize( _NELEM );
for( unsigned IEL = 0; IEL < _NELEM; ++IEL ) {
std::getline( f, line );
_IZ[IEL] = std::stoul( line.substr( line.find_first_of( "=" ) + 1, 3 ) );
_STF[IEL] = std::stod( line.substr( line.find_last_of( "=" ) + 1, 3 ) );
ZTOT += _STF[IEL] * _IZ[IEL];
ATW += ELAW[_IZ[IEL]] * _STF[IEL];
}
std::getline( f, line );
std::getline( f, line );
double EXPOT = std::stod( line.substr( line.find_first_of( "=" ) + 1, 16 ) );
std::getline( f, line );
_NOSC = std::stoul( line.substr( line.find_first_of( "=" ) + 1, line.back() ) );
std::getline( f, line );
std::getline( f, line );
std::getline( f, line );
double EXPOTR = 0.0;
_F.resize( _NOSC );
_WRI.resize( _NOSC );
double buffer;
for( unsigned I = 0; I < _NOSC; ++I ) {
f >> buffer >> _F[I] >> buffer >> _WRI[I];
EXPOTR += _F[I] * std::log( _WRI[I] );
}
f.close();
EXPOTR = std::exp( EXPOTR / ZTOT );
if( std::fabs( EXPOT - EXPOTR ) > 1.0e-5 * EXPOT ) {
ErrorMessages::Error( "Inconsistent oscillator data", CURRENT_FUNCTION );
}
unsigned NDAT = _NA;
ELINIT( _IZ, _STF );
DCSEL0( e );
dxse.resize( NDAT );
dxsi.resize( NDAT );
dxs.resize( NDAT );
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxse[IA] = DCSEL( _XMU[IA] );
}
if( ncor == 1u || ncor == 2u ) { // 1: Inelastic DCS = elastic DCS/Z | 2: Morse's incoherent scattering function.
ErrorMessages::Error( "Unsupported scattering model", CURRENT_FUNCTION );
}
else if( ncor == 3u ) { // Sternheimer -Liljequist GOS model.
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxsi[IA] = DCSIN3( e, _XMU[IA] );
}
}
else { // No inelastic scattering correction
for( unsigned IA = 0; IA < NDAT; ++IA ) {
dxsi[IA] = 0.0e0;
}
}
for( unsigned IA = 0; IA < _NA; ++IA ) {
dxs[IA] = dxse[IA] + dxsi[IA];
}
}
void ICRU::ELINIT( std::vector<double> IZ, std::vector<double> STF ) {
unsigned IE = 0;
double FGRID = 10.0e0;
unsigned IGRID = 10;
while( IE < _NE ) {
double EV = EGRD[IGRID] * FGRID;
if( IGRID == EGRD.size() - 1 ) {
IGRID = 0;
FGRID *= 10.0e0;
}
_ET[IE] = EV;
_ETL[IE] = std::log( _ET[IE] );
IE++;
IGRID++;
}
for( unsigned i = 0; i < _NA; ++i ) {
_XMUL[i] = std::log( std::max( _XMU[i], 1.0e-35 ) );
}
for( unsigned IE = 0; IE < _NE; ++IE ) {
_ECS[IE] = 0.0e0;
_ETCS1[IE] = 0.0e0;
_ETCS2[IE] = 0.0e0;
_PCS[IE] = 0.0e0;
_PTCS1[IE] = 0.0e0;
_PTCS2[IE] = 0.0e0;
for( unsigned IA = 0; IA < _NA; ++IA ) {
_EDCS( IE, IA ) = 0.0e0;
_PDCS( IE, IA ) = 0.0e0;
}
}
for( unsigned IEL = 0; IEL < _NELEM; ++IEL ) {
int IZZ = IZ[IEL];
double STFF = STF[IEL];
unsigned NS = IZ[IEL];
if( NS > 999 ) NS = 999;
// Electrons.
std::string PATHEL = "../data/elsep-database/";
std::stringstream ss;
ss << std::setw( 3 ) << std::setfill( '0' ) << NS;
std::string FILE1 = PATHEL + "eeldx" + ss.str() + ".tab";
_DCSI.resize( _NA );
std::ifstream file( FILE1 );
for( unsigned IE = 0; IE < _NE; ++IE ) {
int IELEC, IZR;
double ENR, CSE, TCS1E, TCS2E;
file >> IELEC >> IZR >> ENR >> CSE >> TCS1E >> TCS2E;
if( particle != ParticleType::ELECTRON || IZR != IZZ || std::fabs( ENR - _ET[IE] ) > 1.0e-3 ) {
ErrorMessages::Error( "Corrupt data file", CURRENT_FUNCTION );
}
for( unsigned IA = 0; IA < _NA; ++IA ) file >> _DCSI[IA];
_ECS[IE] = _ECS[IE] + STFF * CSE;
_ETCS1[IE] = _ETCS1[IE] + STFF * TCS1E;
_ETCS2[IE] = _ETCS2[IE] + STFF * TCS2E;
for( unsigned IA = 0; IA < _NA; ++IA ) {
_EDCS( IE, IA ) = _EDCS( IE, IA ) + STFF * _DCSI[IA];
}
}
file.close();
FILE1 = PATHEL + "peldx" + ss.str() + ".tab";
file.open( FILE1 );
for( unsigned IE = 0; IE < _NE; ++IE ) {
int IELEC, IZR;
double ENR, CSP, TCS1P, TCS2P;
file >> IELEC >> IZR >> ENR >> CSP >> TCS1P >> TCS2P;
if( IELEC != +1 || IZR != IZZ || std::fabs( ENR - _ET[IE] ) > 1.0e-3 ) {
ErrorMessages::Error( "Corrupt data file", CURRENT_FUNCTION );
}
for( unsigned IA = 0; IA < _NA; ++IA ) file >> _DCSI[IA];
_PCS[IE] = _PCS[IE] + STFF * CSP;
_PTCS1[IE] = _PTCS1[IE] + STFF * TCS1P;
_PTCS2[IE] = _PTCS2[IE] + STFF * TCS2P;
for( unsigned IA = 0; IA < _NA; ++IA ) {
_PDCS( IE, IA ) += STFF * _DCSI[IA];
}
}
file.close();
}
}
void ICRU::DCSEL0( double E ) {
if( E < 49.9999e0 || E > 1.0001e8 ) {
ErrorMessages::Error( "Outside of supported energy range", CURRENT_FUNCTION );
}
double EL = std::log( E );
unsigned JE;
FINDI( _ETL, EL, _NE, JE );
std::vector<double> Y( _NE );
for( unsigned IA = 0; IA < _NA; ++IA ) {
for( unsigned IE = 0; IE < _NE; ++IE ) {
if( particle == ParticleType::ELECTRON )
Y[IE] = std::log( _EDCS( IE, IA ) );
else if( particle == ParticleType::POSITRON )
Y[IE] = std::log( _PDCS( IE, IA ) );
else {
ErrorMessages::Error( "Unsupported particle type", CURRENT_FUNCTION );
}
}
std::vector<double> A, B, C, D;
SPLINE( _ETL, Y, A, B, C, D, 0.0e0, 0.0e0, _NE );
_DCSIL[IA] = A[JE] + EL * ( B[JE] + EL * ( C[JE] + EL * D[JE] ) );
_DCSI[IA] = std::exp( _DCSIL[IA] );
}
std::vector<double> AA, BB, CC, DD;
SPLINE( _XMUL, _DCSIL, AA, BB, CC, DD, 0.0e0, 0.0e0, _NA );
for( unsigned IE = 0; IE < _NE; ++IE ) {
if( particle == ParticleType::ELECTRON )
Y[IE] = std::log( _ECS[IE] );
else if( particle == ParticleType::POSITRON )
Y[IE] = std::log( _PCS[IE] );
else {
ErrorMessages::Error( "Unsupported particle type", CURRENT_FUNCTION );
}
}
std::vector<double> A, B, C, D;
SPLINE( _ETL, Y, A, B, C, D, 0.0e0, 0.0e0, _NE );
_CSI = std::exp( A[JE] + EL * ( B[JE] + EL * ( C[JE] + EL * D[JE] ) ) );
for( unsigned IE = 0; IE < _NE; ++IE ) {
if( particle == ParticleType::ELECTRON )
Y[IE] = std::log( _ETCS1[IE] );
else if( particle == ParticleType::POSITRON )
Y[IE] = std::log( _PTCS1[IE] );
else {
std::cerr << "Unsupported particle type" << std::endl;
exit( EXIT_FAILURE );
}
}
SPLINE( _ETL, Y, A, B, C, D, 0.0e0, 0.0e0, _NE );
_TCS1I = std::exp( A[JE] + EL * ( B[JE] + EL * ( C[JE] + EL * D[JE] ) ) );
for( unsigned IE = 0; IE < _NE; ++IE ) {
if( particle == ParticleType::ELECTRON )
Y[IE] = std::log( _ETCS2[IE] );
else if( particle == ParticleType::POSITRON )
Y[IE] = std::log( _PTCS2[IE] );
else {
ErrorMessages::Error( "Unsupported particle type", CURRENT_FUNCTION );
}
SPLINE( _ETL, Y, A, B, C, D, 0.0e0, 0.0e0, _NE );
_TCS2I = std::exp( A[JE] + EL * ( B[JE] + EL * ( C[JE] + EL * D[JE] ) ) );
}
}
double ICRU::DCSEL( double RMU ) {
double RMUL = std::log( std::max( RMU, 1.0e-35 ) );
unsigned I;
FINDI( _XMUL, RMUL, _NA, I );
return std::exp( _DCSIL[I] + ( _DCSIL[I + 1] - _DCSIL[I] ) * ( ( RMUL - _XMUL[I] ) / ( _XMUL[I + 1] - _XMUL[I] ) ) );
}
void ICRU::SPLINE( const std::vector<double>& X,
const std::vector<double>& Y,
std::vector<double>& A,
std::vector<double>& B,
std::vector<double>& C,
std::vector<double>& D,
double S1,
double SN,
unsigned N ) {
if( N < 4 ) {
ErrorMessages::Error( "Too few points", CURRENT_FUNCTION );
}
unsigned N1 = N - 2;
unsigned N2 = N - 3;
A.resize( N );
B.resize( N );
C.resize( N );
D.resize( N );
for( unsigned I = 0; I < N1; ++I ) {
if( X[I + 1] - X[I] < 1.0e-13 ) {
std::cerr << "x values not increasing" << std::endl;
exit( EXIT_FAILURE );
}
A[I] = X[I + 1] - X[I];
D[I] = ( Y[I + 1] - Y[I] ) / A[I];
}
for( unsigned I = 0; I < N2; ++I ) {
B[I] = 2.0e0 * ( A[I] + A[I + 1] );
unsigned K = N1 - I;
D[K] = 6.0e0 * ( D[K] - D[K - 1] );
}
D[1] -= A[0] * S1;
D[N1] = D[N1] - A[N1] * SN;
for( unsigned I = 1; I < N2 + 1; ++I ) {
double R = A[I] / B[I - 1];
B[I] -= R * A[I];
D[I + 1] -= R * D[I];
}
D[N1] = D[N1] / B[N2];
for( unsigned I = 1; I < N2; ++I ) {
unsigned K = N1 - I;
D[K] = ( D[K] - A[K] * D[K + 1] ) / B[K - 1];
}
D[N - 1] = SN;
double SI1 = S1;
for( unsigned I = 0; I < N1; ++I ) {
double SI = SI1;
SI1 = D[I + 1];
double H = A[I];
double HI = 1.0e0 / H;
A[I] = ( HI / 6.0e0 ) * ( SI * X[I + 1] * X[I + 1] * X[I + 1] - SI1 * X[I] * X[I] * X[I] ) + HI * ( Y[I] * X[I + 1] - Y[I + 1] * X[I] ) +
( H / 6.0e0 ) * ( SI1 * X[I] - SI * X[I + 1] );
B[I] = ( HI / 2.0e0 ) * ( SI1 * X[I] * X[I] - SI * X[I + 1] * X[I + 1] ) + HI * ( Y[I + 1] - Y[I] ) + ( H / 6.0e0 ) * ( SI - SI1 );