mirror of
https://github.com/gwm17/catima.git
synced 2024-11-26 12:08:52 -05:00
Andrej Prochazka
GitHub
commit
f21646b132
|
|
@ -72,7 +72,7 @@ if(GLOBAL)
|
|||
file(GLOB GLOBAL_SOURCES global/*.c)
|
||||
LIST (APPEND SOURCES ${GLOBAL_SOURCES})
|
||||
endif(GLOBAL)
|
||||
file(GLOB HEADERS *.h)
|
||||
file(GLOB HEADERS *.h libs/*.h)
|
||||
|
||||
add_library(catima ${SOURCES})
|
||||
set_target_properties(catima PROPERTIES
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ double dedx_n(const Projectile &p, const Target &t){
|
|||
return sn;
|
||||
}
|
||||
|
||||
double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
|
||||
double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c){
|
||||
double w;
|
||||
double sum=0.0;
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
|
|
@ -72,7 +72,7 @@ double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
|
|||
return sum;
|
||||
}
|
||||
|
||||
double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
|
||||
double bethek_dedx_e(const Projectile &p, const Target &t, const Config &c, double I){
|
||||
assert(t.Z>0 && p.Z>0);
|
||||
assert(t.A>0 && p.A>0);
|
||||
assert(p.T>0.0);
|
||||
|
|
@ -447,10 +447,11 @@ double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){
|
|||
double w;
|
||||
double sum=0.0;
|
||||
bool use95 = c.low_energy == low_energy_types::srim_95;
|
||||
double T = p.T;
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto t = mat.get_element(i);
|
||||
w = mat.weight_fraction(i);
|
||||
sum += w*srim_dedx_e(p.Z,t.Z,p.T, use95)/t.A;
|
||||
sum += w*srim_dedx_e(p.Z,t.Z,T, use95)/t.A;
|
||||
}
|
||||
return 100*sum*Avogadro; // returning MeV/g/cm2
|
||||
}
|
||||
|
|
@ -476,7 +477,7 @@ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
|
|||
return factor*beta2/fine_structure/fine_structure;
|
||||
}
|
||||
|
||||
double angular_scattering_variance(Projectile &p, Target &t){
|
||||
double angular_scattering_variance(const Projectile &p, const Target &t){
|
||||
if(p.T<=0)return 0.0;
|
||||
double e=p.T;
|
||||
double _p = p_from_T(e,p.A);
|
||||
|
|
@ -485,7 +486,7 @@ double angular_scattering_variance(Projectile &p, Target &t){
|
|||
return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
|
||||
}
|
||||
|
||||
/// radioation lengths are taken frm Particle Data Group 2014
|
||||
/// radiation lengths are taken from Particle Data Group 2014
|
||||
double radiation_length(int z, double m){
|
||||
double lr = 0;
|
||||
if(z==1){return 63.04;}
|
||||
|
|
@ -560,7 +561,7 @@ double precalculated_lindhard_X(const Projectile &p){
|
|||
return v1+(dif*da/ls_coefficients::a_rel_increase);
|
||||
}
|
||||
|
||||
double dedx_variance(Projectile &p, Target &t, const Config &c){
|
||||
double dedx_variance(const Projectile &p, const Target &t, const Config &c){
|
||||
double gamma = gamma_from_T(p.T);
|
||||
double cor=0;
|
||||
double beta = beta_from_T(p.T);
|
||||
|
|
|
|||
|
|
@ -31,7 +31,7 @@ namespace catima{
|
|||
/**
|
||||
* returns energy loss straggling
|
||||
*/
|
||||
double dedx_variance(Projectile &p, Target &t, const Config &c=default_config);
|
||||
double dedx_variance(const Projectile &p, const Target &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns reduced energy loss unit for projectile-target combination
|
||||
|
|
@ -43,10 +43,17 @@ namespace catima{
|
|||
* @brief bethek_dedx_e - electronics stopping power
|
||||
* @return stopping power
|
||||
*/
|
||||
double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
|
||||
double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config);
|
||||
double bethek_dedx_e(const Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
|
||||
double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates barkas effect
|
||||
*/
|
||||
double bethek_barkas(double zp_eff,double eta, double zt);
|
||||
|
||||
/**
|
||||
* calculates density effect
|
||||
*/
|
||||
double bethek_density_effect(double beta, int zt);
|
||||
|
||||
/**
|
||||
|
|
@ -94,7 +101,7 @@ namespace catima{
|
|||
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
|
||||
|
||||
|
||||
double angular_scattering_variance(Projectile &p, Target &t);
|
||||
double angular_scattering_variance(const Projectile &p, const Target &t);
|
||||
|
||||
/**
|
||||
* returns radiation length of the (M,Z) material
|
||||
|
|
|
|||
162
catima.cpp
162
catima.cpp
|
|
@ -26,17 +26,17 @@ bool operator==(const Config &a, const Config&b){
|
|||
}
|
||||
|
||||
|
||||
double dedx(Projectile &p, const Material &mat, const Config &c){
|
||||
double dedx(const Projectile &p, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
if(p.T<=0)return 0.0;
|
||||
sum += dedx_n(p,mat);
|
||||
double se=0;
|
||||
if(p.T<=10){
|
||||
se = sezi_dedx_e(p,mat);
|
||||
se = sezi_dedx_e(p,mat,c );
|
||||
}
|
||||
else if(p.T>10 && p.T<30){
|
||||
double factor = 0.05 * ( p.T - 10.0 );
|
||||
se = (1-factor)*sezi_dedx_e(p,mat) + factor*bethek_dedx_e(p,mat,c);
|
||||
se = (1-factor)*sezi_dedx_e(p,mat,c) + factor*bethek_dedx_e(p,mat,c);
|
||||
}
|
||||
else {
|
||||
se = bethek_dedx_e(p,mat,c);
|
||||
|
|
@ -46,46 +46,44 @@ double dedx(Projectile &p, const Material &mat, const Config &c){
|
|||
return sum;
|
||||
}
|
||||
|
||||
double domega2dx(Projectile &p, double T, const Material &mat, const Config &c){
|
||||
double domega2dx(const Projectile &p, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto t= mat.get_element(i);
|
||||
double w = mat.weight_fraction(i);
|
||||
p.T = T;
|
||||
sum += w*dedx_variance(p,t,c);
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
|
||||
double da2dx(const Projectile &p, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto t = mat.get_element(i);
|
||||
double w = mat.weight_fraction(i);
|
||||
p.T = T;
|
||||
sum += w*angular_scattering_variance(p,t);
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
|
||||
double range(Projectile &p, const Material &t, const Config &c){
|
||||
double range(const Projectile &p, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
return range_spline(p.T);
|
||||
}
|
||||
|
||||
double dedx_from_range(Projectile &p, const Material &t, const Config &c){
|
||||
double dedx_from_range(const Projectile &p, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
return p.A/range_spline.derivative(p.T);
|
||||
}
|
||||
|
||||
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
|
||||
std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
|
|
@ -102,42 +100,42 @@ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T,
|
|||
return dedx;
|
||||
}
|
||||
|
||||
double range_straggling(Projectile &p, double T, const Material &t, const Config &c){
|
||||
double range_straggling(const Projectile &p, double T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
return sqrt(range_straggling_spline(T));
|
||||
}
|
||||
|
||||
double range_variance(Projectile &p, double T, const Material &t, const Config &c){
|
||||
double range_variance(const Projectile &p, double T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
return range_straggling_spline(T);
|
||||
}
|
||||
|
||||
double domega2de(Projectile &p, double T, const Material &t, const Config &c){
|
||||
double domega2de(const Projectile &p, double T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
return range_straggling_spline.derivative(T);
|
||||
}
|
||||
|
||||
double da2de(Projectile &p, double T, const Material &t, const Config &c){
|
||||
double da2de(const Projectile &p, double T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
return angular_variance_spline.derivative(T);
|
||||
}
|
||||
|
||||
double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
|
||||
}
|
||||
|
||||
double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){
|
||||
double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
|
|
@ -173,14 +171,14 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){
|
|||
return -1;
|
||||
}
|
||||
|
||||
double energy_out(Projectile &p, const Material &t, const Config &c){
|
||||
double energy_out(const Projectile &p, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
return energy_out(p.T,t.thickness(),range_spline);
|
||||
}
|
||||
|
||||
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
|
||||
std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
|
|
@ -199,7 +197,24 @@ std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, cons
|
|||
return eout;
|
||||
}
|
||||
|
||||
Result calculate(Projectile &p, const Material &t, const Config &c){
|
||||
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
|
||||
double res;
|
||||
std::vector<double> values;
|
||||
values.reserve(max_datapoints);
|
||||
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
|
||||
res = integrator.integrate(function,Ezero,energy_table(0));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
values.push_back(res);
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
}
|
||||
return values;
|
||||
}
|
||||
|
||||
Result calculate(Projectile p, const Material &t, const Config &c){
|
||||
Result res;
|
||||
double T = p.T;
|
||||
if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
|
||||
|
|
@ -215,7 +230,7 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
|
|||
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
if(res.Eout<Ezero){
|
||||
res.dEdxo = 0.0;
|
||||
res.sigma_a = 0.0;
|
||||
|
|
@ -224,24 +239,18 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
|
|||
}
|
||||
else{
|
||||
res.dEdxo = p.A/range_spline.derivative(res.Eout);
|
||||
|
||||
#ifdef THIN_TARGET_APPROXIMATION
|
||||
if(thin_target_limit*res.Ein<res.Eout){
|
||||
double edif = (res.Ein-res.Eout);
|
||||
double s1 = range_straggling_spline.derivative(T);
|
||||
double s2 = range_straggling_spline.derivative(res.Eout);
|
||||
res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A;
|
||||
|
||||
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
s1 = angular_variance_spline.derivative(T);
|
||||
s2 = angular_variance_spline.derivative(res.Eout);
|
||||
res.sigma_a = sqrt(0.5*(s1+s2)*edif);
|
||||
}
|
||||
else{
|
||||
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
|
||||
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
|
||||
}
|
||||
|
||||
|
|
@ -260,26 +269,50 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
|
|||
}
|
||||
res.sigma_r = sqrt(range_straggling_spline(T));
|
||||
res.Eloss = (res.Ein - res.Eout)*p.A;
|
||||
|
||||
double rrange = std::min(res.range/t.density(), t.thickness_cm());
|
||||
auto fx2p = [&](double x)->double{
|
||||
double range = range_spline(T);
|
||||
double e =energy_out(T,x*t.density(),range_spline);
|
||||
return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
|
||||
};
|
||||
|
||||
//res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3, 1e-3,1);
|
||||
res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange);
|
||||
res.sigma_x = sqrt(res.sigma_x);
|
||||
|
||||
auto fx1p = [&](double x)->double{
|
||||
double e =energy_out(T,x*t.density(),range_spline);
|
||||
return (rrange-x)*da2dx(p(e), t, c)*t.density();
|
||||
};
|
||||
|
||||
//res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,1);
|
||||
res.cov = integrator.integrate(fx1p,0, rrange);
|
||||
#ifdef REACTIONS
|
||||
res.sp = nonreaction_rate(p,t,c);
|
||||
#endif
|
||||
return res;
|
||||
}
|
||||
|
||||
MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
|
||||
MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c){
|
||||
MultiResult res;
|
||||
double e = p.T;
|
||||
res.total_result.Ein = e;
|
||||
res.results.reserve(layers.num());
|
||||
|
||||
for(auto&m:layers.get_materials()){
|
||||
Result r = calculate(p,m,e,c);
|
||||
e = r.Eout;
|
||||
res.total_result.sigma_a += r.sigma_a*r.sigma_a;
|
||||
res.total_result.Eloss += r.Eloss;
|
||||
res.total_result.sigma_E += r.sigma_E*r.sigma_E;
|
||||
res.total_result.tof += r.tof;
|
||||
res.total_result.Eout = r.Eout;
|
||||
double a2 = res.total_result.sigma_a;
|
||||
res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov)
|
||||
+ (a2*m.thickness_cm()*m.thickness_cm())
|
||||
+ r.sigma_x*r.sigma_x;
|
||||
//res.total_result.sigma_x += (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x;
|
||||
res.total_result.cov += a2*m.thickness_cm() + r.cov;
|
||||
res.total_result.sigma_a += r.sigma_a*r.sigma_a;
|
||||
#ifdef REACTIONS
|
||||
res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1;
|
||||
#endif
|
||||
|
|
@ -288,10 +321,13 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
|
|||
if(e>Ezero){
|
||||
res.total_result.sigma_a = sqrt(res.total_result.sigma_a);
|
||||
res.total_result.sigma_E = sqrt(res.total_result.sigma_E);
|
||||
res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
|
||||
|
||||
}
|
||||
else{
|
||||
res.total_result.sigma_a = 0.0;
|
||||
res.total_result.sigma_E = 0.0;
|
||||
res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
|
@ -302,24 +338,6 @@ Result calculate(double pa, int pz, double T, double ta, double tz, double thick
|
|||
return calculate(p(T),m);
|
||||
}
|
||||
|
||||
|
||||
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
|
||||
double res;
|
||||
std::vector<double> values;
|
||||
values.reserve(max_datapoints);
|
||||
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
|
||||
res = integrator.integrate(function,Ezero,energy_table(0));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
values.push_back(res);
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
}
|
||||
return values;
|
||||
}
|
||||
/*
|
||||
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
||||
DataPoint dp(p,t,c);
|
||||
dp.range.resize(max_datapoints);
|
||||
|
|
@ -329,49 +347,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
return 1.0/dedx(p(x),t,c);
|
||||
};
|
||||
auto fomega = [&](double x)->double{
|
||||
//return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t),3);
|
||||
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
|
||||
};
|
||||
|
||||
double res;
|
||||
//calculate 1st point to have i-1 element ready for loop
|
||||
//res = integrator.integrate(fdedx,Ezero,energy_table(0));
|
||||
//res = p.A*res;
|
||||
//dp.range[0] = res;
|
||||
dp.range[0] = 0.0;
|
||||
|
||||
res = da2dx(p,energy_table(0),t)*res;
|
||||
dp.angular_variance[0] = 0.0;
|
||||
|
||||
//res = integrator.integrate(fomega,Ezero,energy_table(0));
|
||||
//res = p.A*res;
|
||||
dp.range_straggling[0]=0.0;
|
||||
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
|
||||
dp.range[i] = res + dp.range[i-1];
|
||||
res = da2dx(p,energy_table(i),t)*res;
|
||||
dp.angular_variance[i] = res + dp.angular_variance[i-1];
|
||||
|
||||
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
|
||||
res = p.A*res;
|
||||
dp.range_straggling[i] = res + dp.range_straggling[i-1];
|
||||
}
|
||||
return dp;
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
||||
DataPoint dp(p,t,c);
|
||||
dp.range.resize(max_datapoints);
|
||||
dp.range_straggling.resize(max_datapoints);
|
||||
dp.angular_variance.resize(max_datapoints);
|
||||
auto fdedx = [&](double x)->double{
|
||||
return 1.0/dedx(p(x),t,c);
|
||||
};
|
||||
auto fomega = [&](double x)->double{
|
||||
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
|
||||
return domega2dx(p(x),t,c)/catima::power(dedx(p(x),t,c),3);
|
||||
};
|
||||
|
||||
//double res=0.0;
|
||||
|
|
@ -390,7 +366,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
for(int i=1;i<max_datapoints;i++){
|
||||
double res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
|
||||
dp.range[i] = res + dp.range[i-1];
|
||||
res = da2dx(p,energy_table(i),t)*res;
|
||||
res = da2dx(p(energy_table(i)),t)*res;
|
||||
dp.angular_variance[i] = res + dp.angular_variance[i-1];
|
||||
|
||||
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
|
||||
|
|
@ -402,15 +378,13 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
|
||||
double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){
|
||||
double res;
|
||||
//double beta_in = beta_from_T(p.T);
|
||||
//double beta_out = beta_from_T(Eout);
|
||||
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
|
||||
res = integrator.integrate(function,Eout,p.T);
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
return res;
|
||||
}
|
||||
|
||||
std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c){
|
||||
std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c){
|
||||
std::pair<double, double> res{1.0,1.0};
|
||||
if(t.density()<= 0.0 || t.thickness()<=0){
|
||||
return res;
|
||||
|
|
|
|||
43
catima.h
43
catima.h
|
|
@ -37,7 +37,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return dEdx
|
||||
*/
|
||||
double dedx(Projectile &p, const Material &mat, const Config &c=default_config);
|
||||
double dedx(const Projectile &p, const Material &mat, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculate energy loss straggling variance for projectile-Material combination
|
||||
|
|
@ -45,12 +45,12 @@ namespace catima{
|
|||
* @param t - Material
|
||||
* @return dOmega^2/dx
|
||||
*/
|
||||
double domega2dx(Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
double domega2dx(const Projectile &p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates variance of angular scattering of Projectile p on Material m
|
||||
*/
|
||||
double da2dx(Projectile &p, double T, const Material &m, const Config &c=default_config);
|
||||
double da2dx(const Projectile &p, const Material &m, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the range of the Projectile in Material calculated from range spline
|
||||
|
|
@ -58,7 +58,7 @@ namespace catima{
|
|||
* @param t - Material
|
||||
* @return range
|
||||
*/
|
||||
double range(Projectile &p, const Material &t, const Config &c=default_config);
|
||||
double range(const Projectile &p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the dEdx calculated from range spline as derivative
|
||||
|
|
@ -66,7 +66,7 @@ namespace catima{
|
|||
* @param t - Material
|
||||
* @return range
|
||||
*/
|
||||
double dedx_from_range(Projectile &p, const Material &t, const Config &c=default_config);
|
||||
double dedx_from_range(const Projectile &p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the dEdx calculated from range spline as derivative
|
||||
|
|
@ -75,7 +75,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return range
|
||||
*/
|
||||
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
|
||||
std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the range straggling of the Projectile in Material from spline
|
||||
|
|
@ -84,7 +84,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return range straggling
|
||||
*/
|
||||
double range_straggling(Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
double range_straggling(const Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the range variance of the Projectile in Material from spline
|
||||
|
|
@ -93,7 +93,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return range straggling
|
||||
*/
|
||||
double range_variance(Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
double range_variance(const Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the range variance per dE, calculated as derivative of range variance spline
|
||||
|
|
@ -102,7 +102,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return range variance / dE
|
||||
*/
|
||||
double domega2de(Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
double domega2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the angular variance per dE, calculated as derivative of angular variance spline
|
||||
|
|
@ -111,7 +111,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return angular variance / dE
|
||||
*/
|
||||
double da2de(Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
double da2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates angular scattering in the material from difference of incoming a nd outgoing energies
|
||||
|
|
@ -121,7 +121,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return angular straggling
|
||||
*/
|
||||
double angular_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates Energy straggling in the material from difference of incoming a nd outgoing energies
|
||||
|
|
@ -131,7 +131,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return angular straggling
|
||||
*/
|
||||
double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
|
||||
double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates outcoming energy from range spline
|
||||
|
|
@ -149,7 +149,7 @@ namespace catima{
|
|||
* @param T - incoming energy
|
||||
* @return outcoming energy after the material in Mev/u
|
||||
*/
|
||||
double energy_out(Projectile &p, const Material &t, const Config &c=default_config);
|
||||
double energy_out(const Projectile &p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates outcoming energy
|
||||
|
|
@ -158,7 +158,7 @@ namespace catima{
|
|||
* @param T - incoming energy vector
|
||||
* @return outcoming energy after the material in Mev/u
|
||||
*/
|
||||
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
|
||||
std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates all observables for projectile passing material
|
||||
|
|
@ -166,8 +166,8 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return structure of Result
|
||||
*/
|
||||
Result calculate(Projectile &p, const Material &t, const Config &c=default_config);
|
||||
inline Result calculate(Projectile &p, const Material &t, double T, const Config &c=default_config){
|
||||
Result calculate(Projectile p, const Material &t, const Config &c=default_config);
|
||||
inline Result calculate(Projectile p, const Material &t, double T, const Config &c=default_config){
|
||||
p.T = T;
|
||||
return calculate(p, t, c);
|
||||
}
|
||||
|
|
@ -186,15 +186,14 @@ namespace catima{
|
|||
* @return results stored in MultiResult structure
|
||||
*
|
||||
*/
|
||||
MultiResult calculate(Projectile &p, const Layers &layers, const Config &c=default_config);
|
||||
inline MultiResult calculate(Projectile &p, double T, const Layers &layers, const Config &c=default_config){
|
||||
p.T = T;
|
||||
return calculate(p, layers, c);
|
||||
MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
|
||||
inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
|
||||
return calculate(p(T), layers, c);
|
||||
}
|
||||
|
||||
|
||||
/// this calculate tof spline, at the moment it is not used
|
||||
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c=default_config);
|
||||
std::vector<double> calculate_tof(const Projectile p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates TOF of the Projectile in Material
|
||||
|
|
@ -206,7 +205,7 @@ namespace catima{
|
|||
/**
|
||||
* returns energy magnification after passing material t
|
||||
*/
|
||||
std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c=default_config);
|
||||
std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c=default_config);
|
||||
|
||||
class DataPoint;
|
||||
/**
|
||||
|
|
|
|||
|
|
@ -7,6 +7,8 @@
|
|||
|
||||
namespace catima{
|
||||
integrator_type integrator;
|
||||
integrator_adaptive_type integrator_adaptive;
|
||||
|
||||
#ifdef GSL_INTEGRATION
|
||||
double funcwrapper3(double x, void *_c){
|
||||
std::function<double(double)> *f = (std::function<double(double)> *)_c;
|
||||
|
|
|
|||
132
integrator.h
132
integrator.h
|
|
@ -23,11 +23,16 @@
|
|||
//#ifdef USE_THREADS
|
||||
//#include <mutex>
|
||||
//#endif
|
||||
|
||||
#include "catima/glq_integrator.h"
|
||||
#include "catima/gkq_integrator.h"
|
||||
#ifdef GSL_INTEGRATION
|
||||
#include "gsl/gsl_integration.h"
|
||||
#endif
|
||||
|
||||
|
||||
using integrators::GaussLegendreIntegration;
|
||||
using integrators::GaussKronrodIntegration;
|
||||
|
||||
namespace catima{
|
||||
|
||||
#ifdef GSL_INTEGRATION
|
||||
|
|
@ -52,136 +57,15 @@ class IntegratorGSL{
|
|||
};
|
||||
#endif
|
||||
|
||||
template<int order>
|
||||
struct GL_data{
|
||||
};
|
||||
|
||||
template<int order>
|
||||
class GaussLegendreIntegration{
|
||||
public:
|
||||
template<typename F>
|
||||
double integrate(F f, double a, double b) const;
|
||||
template<typename F>
|
||||
double operator()(F f, double a, double b) const {return integrate(f, a, b);}
|
||||
double w(int i) const {return GL_data<order>::w()[i];}
|
||||
double x(int i) const {return GL_data<order>::x()[i];}
|
||||
int n() const {return order;}
|
||||
std::array<double,order> get_points(double a = -1.0, double b = 1.0)const;
|
||||
};
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussLegendreIntegration<order>::integrate(F f, double a, double b) const{
|
||||
double res=0.0;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
if(order%2){res+= w(0) * (f(p*x(0) + q));} // in case odd-order
|
||||
for(int i=order%2;i<order/2 + order%2;i++){
|
||||
res += w(i) * (f(p*x(i) + q) + f(-p*x(i) + q));
|
||||
}
|
||||
return p*res;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
std::array<double,order> GaussLegendreIntegration<order>::get_points(double a, double b)const{
|
||||
std::array<double,order> points;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
int num = (order/2);
|
||||
for(int i=0;i< num;i++){
|
||||
points[num-i-1] = -p*x(i) + q;
|
||||
points[num+i] = p*x(i) + q;
|
||||
}
|
||||
return points;
|
||||
}
|
||||
|
||||
// order = 8
|
||||
template<>
|
||||
struct GL_data<8>{
|
||||
static const std::array<double,4>& x(){
|
||||
static const std::array<double,4> _x =
|
||||
{0.1834346424956498049394761,0.5255324099163289858177390,0.7966664774136267395915539,0.9602898564975362316835609};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,4>& w(){
|
||||
static const std::array<double,4> _w =
|
||||
{0.3626837833783619829651504,0.3137066458778872873379622,0.2223810344533744705443560,0.1012285362903762591525314};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 10
|
||||
template<>
|
||||
struct GL_data<10>{
|
||||
static const std::array<double,5>& x(){
|
||||
static const std::array<double,5> _x =
|
||||
{0.1488743389816312108848260,0.4333953941292471907992659,0.6794095682990244062343274,0.8650633666889845107320967,0.9739065285171717200779640};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,5>& w(){
|
||||
static const std::array<double,5> _w =
|
||||
{0.2955242247147528701738930,0.2692667193099963550912269,0.2190863625159820439955349,0.1494513491505805931457763,0.0666713443086881375935688};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 12
|
||||
template<>
|
||||
struct GL_data<12>{
|
||||
static const std::array<double,6>& x(){
|
||||
static const std::array<double,6> _x =
|
||||
{0.1252334085114689154724414,0.3678314989981801937526915,0.5873179542866174472967024,0.7699026741943046870368938,0.9041172563704748566784659,0.9815606342467192506905491};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,6>& w(){
|
||||
static const std::array<double,6> _w =
|
||||
{0.2491470458134027850005624,0.2334925365383548087608499,0.2031674267230659217490645,0.1600783285433462263346525,0.1069393259953184309602547,0.0471753363865118271946160};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 14
|
||||
template<>
|
||||
struct GL_data<14>{
|
||||
static const std::array<double,7>& x(){
|
||||
static const std::array<double,7> _x =
|
||||
{0.1080549487073436620662447,0.3191123689278897604356718,0.5152486363581540919652907,0.6872929048116854701480198,0.8272013150697649931897947,0.9284348836635735173363911,0.9862838086968123388415973};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,7>& w(){
|
||||
static const std::array<double,7> _w =
|
||||
{0.2152638534631577901958764,0.2051984637212956039659241,0.1855383974779378137417166,0.1572031671581935345696019,0.1215185706879031846894148,0.0801580871597602098056333,0.0351194603317518630318329};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 16
|
||||
template<>
|
||||
struct GL_data<16>{
|
||||
static const std::array<double,8>& x(){
|
||||
static const std::array<double,8> _x =
|
||||
{0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,8>& w(){
|
||||
static const std::array<double,8> _w =
|
||||
{0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
#ifdef GSL_INTEGRATION
|
||||
using integrator_type = IntegratorGSL;
|
||||
#else
|
||||
using integrator_type = GaussLegendreIntegration<8>;
|
||||
#endif
|
||||
using integrator_adaptive_type = GaussKronrodIntegration<15>;
|
||||
|
||||
extern integrator_type integrator;
|
||||
extern integrator_adaptive_type integrator_adaptive;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
|
|||
440
libs/gkq_integrator.h
Normal file
440
libs/gkq_integrator.h
Normal file
|
|
@ -0,0 +1,440 @@
|
|||
/*
|
||||
* Copyright(C) 2017, Andrej Prochazka
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU Affero General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU Affero General Public License for more details.
|
||||
|
||||
* You should have received a copy of the GNU Affero General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef GKQ_INTEGRATOR_H
|
||||
#define GKQ_INTEGRATOR_H
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace integrators{
|
||||
|
||||
template<int order>
|
||||
struct GK_data{
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief adaptive integrator
|
||||
* following orders are supported: 7,15,21,31,61
|
||||
* @tparam order Kronrod order
|
||||
*/
|
||||
template<int order>
|
||||
class GaussKronrodIntegration{
|
||||
public:
|
||||
template<typename F>
|
||||
static double integrate(F& f, double a, double b, double eps = 1e-3, double reps=1e-6, int N = 1);
|
||||
|
||||
template<typename F>
|
||||
static double integrate_intervals(F& f, const std::vector<std::pair<double,double>>& intervals, double eps = 1e-3, double reps=1e-6);
|
||||
|
||||
template<typename F>
|
||||
static std::pair<double,double> integrate_nonadaptive(F& f, double a, double b);
|
||||
|
||||
template<typename F>
|
||||
static double integrate_adaptive(F& f, double a, double b, double eps = 1e-3, double reps=1e-6, int level=49);
|
||||
|
||||
static double w(int i) {return GK_data<order>::w()[i];}
|
||||
static double wg(int i) {return GK_data<order>::wg()[i];}
|
||||
static double x(int i) {return GK_data<order>::x()[i];}
|
||||
static int n() {return order;}
|
||||
static std::array<double,order> get_points(double a = -1.0, double b = 1.0);
|
||||
|
||||
private:
|
||||
static const size_t ngauss = (order-1)/2;
|
||||
static const size_t xsize = order/2 + 1;
|
||||
static const bool oddgauss = ((order-1)/2)&1;
|
||||
};
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
std::pair<double,double> GaussKronrodIntegration<order>::integrate_nonadaptive(F& f, double a, double b){
|
||||
double res;
|
||||
double gres = 0.0;
|
||||
double val;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
unsigned int gi = 1; //by default assume odd-order gauss
|
||||
unsigned int ki = 2;
|
||||
|
||||
// 1st kronrod point
|
||||
val = f(p*x(0) + q);
|
||||
res = w(0) * val;
|
||||
if(oddgauss){ // 1st gaus point if odd order gauss
|
||||
gres = wg(0) * val;
|
||||
gi = 2;
|
||||
ki = 1;
|
||||
}
|
||||
|
||||
for(unsigned int i=gi;i<xsize;i+=2){
|
||||
val = f(p*x(i) + q);
|
||||
res += w(i) * val;
|
||||
gres += wg(i/2)*val;
|
||||
|
||||
val = f(-p*x(i) + q);
|
||||
res += w(i) * val;
|
||||
gres += wg(i/2)*val;
|
||||
}
|
||||
|
||||
for(unsigned int i=ki;i<xsize;i+=2){
|
||||
res += w(i) * (f(p*x(i) + q) + f(-p*x(i) + q));
|
||||
}
|
||||
|
||||
double err = std::max( std::abs(gres-res), std::numeric_limits<double>::epsilon());
|
||||
|
||||
return std::make_pair(p*res, p*err);
|
||||
}
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussKronrodIntegration<order>::integrate_adaptive(F& f, double a, double b, double eps, double rprec, int level){
|
||||
double result = 0.0;
|
||||
double err = 0.0;
|
||||
const double numlimit = 10*std::numeric_limits<double>::epsilon();
|
||||
|
||||
auto r = integrate_nonadaptive(f, a, b);
|
||||
result = r.first;
|
||||
err = r.second;
|
||||
//printf("level = %d, I=%lf, e=%lf, %lf - %lf\n",level,result,err,a,b);
|
||||
if( (std::abs(result) < numlimit) || ((b-a)<numlimit)){
|
||||
return result;
|
||||
}
|
||||
|
||||
double aeps = std::max(rprec*std::abs(result), eps);
|
||||
if((aeps<numlimit) || (std::abs(result)<aeps)){
|
||||
return result;
|
||||
}
|
||||
|
||||
if( level && (err > aeps)){
|
||||
double mid = 0.5*(a+b);
|
||||
result = integrate_adaptive(f, a, mid, 0.707*aeps, 0, level-1);
|
||||
result += integrate_adaptive(f, mid, b, 0.707*aeps, 0, level-1);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussKronrodIntegration<order>::integrate(F& f, double a, double b, double eps, double reps, int N){
|
||||
double step = (b-a)/N;
|
||||
double result = 0;
|
||||
for(int i=0; i<N;i++){
|
||||
double m =a+(i*step);
|
||||
result+=integrate_adaptive(f,m,m+step, eps/N, reps);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussKronrodIntegration<order>::integrate_intervals(F& f, const std::vector<std::pair<double,double>>& intervals, double eps, double reps){
|
||||
double result = 0;
|
||||
for(auto& i:intervals){
|
||||
result+=integrate_adaptive(f,i.first,i.second, eps/intervals.size(), reps);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
std::array<double,order> GaussKronrodIntegration<order>::get_points(double a, double b){
|
||||
std::array<double,order> points;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
int num = (order/2);
|
||||
for(int i=0;i< num;i++){
|
||||
points[num-i-1] = -p*x(i) + q;
|
||||
points[num+i] = p*x(i) + q;
|
||||
}
|
||||
return points;
|
||||
}
|
||||
|
||||
|
||||
/// weights and abscissas
|
||||
// order = 7
|
||||
template<>
|
||||
struct GK_data<7>{
|
||||
static std::array<double,4> const & x(){
|
||||
static const std::array<double,4> _x =
|
||||
{0.0, 0.4342437493468025580021, 0.774596669241483377036,0.9604912687080202834235};
|
||||
return _x;
|
||||
}
|
||||
static std::array<double,4> const & w(){
|
||||
static const std::array<double,4> _w =
|
||||
{0.4509165386584741423451, 0.4013974147759622229051, 0.2684880898683334407286, 0.1046562260264672651938};
|
||||
return _w;
|
||||
}
|
||||
static std::array<double,4> const & wg(){
|
||||
static const std::array<double,4> _wg =
|
||||
{0.8888888888888888888889, 0.555555555555555555556};
|
||||
return _wg;
|
||||
}
|
||||
};
|
||||
|
||||
//order = 15
|
||||
template<>
|
||||
struct GK_data<15>{
|
||||
static std::array<double,8> const & x(){
|
||||
static const std::array<double,8> _x =
|
||||
{0.0,
|
||||
0.207784955007898467600689403773244913,
|
||||
0.405845151377397166906606412076961463,
|
||||
0.586087235467691130294144838258729598,
|
||||
0.741531185599394439863864773280788407,
|
||||
0.864864423359769072789712788640926201,
|
||||
0.949107912342758524526189684047851262,
|
||||
0.991455371120812639206854697526328517};
|
||||
return _x;
|
||||
}
|
||||
static std::array<double,8> const & w(){
|
||||
static const std::array<double,8> _w =
|
||||
{
|
||||
0.209482141084727828012999174891714264,
|
||||
0.204432940075298892414161999234649085,
|
||||
0.190350578064785409913256402421013683,
|
||||
0.169004726639267902826583426598550284,
|
||||
0.140653259715525918745189590510237920,
|
||||
0.104790010322250183839876322541518017,
|
||||
0.0630920926299785532907006631892042867,
|
||||
0.0229353220105292249637320080589695920};
|
||||
return _w;
|
||||
}
|
||||
static std::array<double,4> const & wg(){
|
||||
static const std::array<double,4> _wg =
|
||||
{ 0.417959183673469387755102040816326531,
|
||||
0.381830050505118944950369775488975134,
|
||||
0.279705391489276667901467771423779582,
|
||||
0.129484966168869693270611432679082018};
|
||||
return _wg;
|
||||
}
|
||||
};
|
||||
|
||||
//order = 21
|
||||
template<>
|
||||
struct GK_data<21>{
|
||||
static std::array<double,11> const & x(){
|
||||
static const std::array<double,11> _x = {
|
||||
0.00000000000000000e+00,
|
||||
1.48874338981631211e-01,
|
||||
2.94392862701460198e-01,
|
||||
4.33395394129247191e-01,
|
||||
5.62757134668604683e-01,
|
||||
6.79409568299024406e-01,
|
||||
7.80817726586416897e-01,
|
||||
8.65063366688984511e-01,
|
||||
9.30157491355708226e-01,
|
||||
9.73906528517171720e-01,
|
||||
9.95657163025808081e-01};
|
||||
return _x;
|
||||
}
|
||||
static std::array<double,11> const & w(){
|
||||
static const std::array<double,11> _w =
|
||||
{
|
||||
1.49445554002916906e-01,
|
||||
1.47739104901338491e-01,
|
||||
1.42775938577060081e-01,
|
||||
1.34709217311473326e-01,
|
||||
1.23491976262065851e-01,
|
||||
1.09387158802297642e-01,
|
||||
9.31254545836976055e-02,
|
||||
7.50396748109199528e-02,
|
||||
5.47558965743519960e-02,
|
||||
3.25581623079647275e-02,
|
||||
1.16946388673718743e-02,
|
||||
};
|
||||
return _w;
|
||||
}
|
||||
static std::array<double,5> const & wg(){
|
||||
static const std::array<double,5> _wg =
|
||||
{
|
||||
2.95524224714752870e-01,
|
||||
2.69266719309996355e-01,
|
||||
2.19086362515982044e-01,
|
||||
1.49451349150580593e-01,
|
||||
6.66713443086881376e-02,};
|
||||
return _wg;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
//order = 31
|
||||
template<>
|
||||
struct GK_data<31>{
|
||||
static std::array<double,16> const & x(){
|
||||
static const std::array<double,16> _x = {
|
||||
0.00000000000000000e+00,
|
||||
1.01142066918717499e-01,
|
||||
2.01194093997434522e-01,
|
||||
2.99180007153168812e-01,
|
||||
3.94151347077563370e-01,
|
||||
4.85081863640239681e-01,
|
||||
5.70972172608538848e-01,
|
||||
6.50996741297416971e-01,
|
||||
7.24417731360170047e-01,
|
||||
7.90418501442465933e-01,
|
||||
8.48206583410427216e-01,
|
||||
8.97264532344081901e-01,
|
||||
9.37273392400705904e-01,
|
||||
9.67739075679139134e-01,
|
||||
9.87992518020485428e-01,
|
||||
9.98002298693397060e-01
|
||||
};
|
||||
return _x;
|
||||
}
|
||||
static std::array<double,16> const & w(){
|
||||
static const std::array<double,16> _w =
|
||||
{
|
||||
1.01330007014791549e-01,
|
||||
1.00769845523875595e-01,
|
||||
9.91735987217919593e-02,
|
||||
9.66427269836236785e-02,
|
||||
9.31265981708253212e-02,
|
||||
8.85644430562117706e-02,
|
||||
8.30805028231330210e-02,
|
||||
7.68496807577203789e-02,
|
||||
6.98541213187282587e-02,
|
||||
6.20095678006706403e-02,
|
||||
5.34815246909280873e-02,
|
||||
4.45897513247648766e-02,
|
||||
3.53463607913758462e-02,
|
||||
2.54608473267153202e-02,
|
||||
1.50079473293161225e-02,
|
||||
5.37747987292334899e-03
|
||||
};
|
||||
return _w;
|
||||
}
|
||||
static std::array<double,8> const & wg(){
|
||||
static const std::array<double,8> _wg =
|
||||
{
|
||||
2.02578241925561273e-01,
|
||||
1.98431485327111576e-01,
|
||||
1.86161000015562211e-01,
|
||||
1.66269205816993934e-01,
|
||||
1.39570677926154314e-01,
|
||||
1.07159220467171935e-01,
|
||||
7.03660474881081247e-02,
|
||||
3.07532419961172684e-02
|
||||
};
|
||||
return _wg;
|
||||
}
|
||||
};
|
||||
|
||||
//order = 61
|
||||
template<>
|
||||
struct GK_data<61>{
|
||||
static std::array<double,31> const & x(){
|
||||
static const std::array<double,31> _x = {
|
||||
0.00000000000000000e+00,
|
||||
5.14718425553176958e-02,
|
||||
1.02806937966737030e-01,
|
||||
1.53869913608583547e-01,
|
||||
2.04525116682309891e-01,
|
||||
2.54636926167889846e-01,
|
||||
3.04073202273625077e-01,
|
||||
3.52704725530878113e-01,
|
||||
4.00401254830394393e-01,
|
||||
4.47033769538089177e-01,
|
||||
4.92480467861778575e-01,
|
||||
5.36624148142019899e-01,
|
||||
5.79345235826361692e-01,
|
||||
6.20526182989242861e-01,
|
||||
6.60061064126626961e-01,
|
||||
6.97850494793315797e-01,
|
||||
7.33790062453226805e-01,
|
||||
7.67777432104826195e-01,
|
||||
7.99727835821839083e-01,
|
||||
8.29565762382768397e-01,
|
||||
8.57205233546061099e-01,
|
||||
8.82560535792052682e-01,
|
||||
9.05573307699907799e-01,
|
||||
9.26200047429274326e-01,
|
||||
9.44374444748559979e-01,
|
||||
9.60021864968307512e-01,
|
||||
9.73116322501126268e-01,
|
||||
9.83668123279747210e-01,
|
||||
9.91630996870404595e-01,
|
||||
9.96893484074649540e-01,
|
||||
9.99484410050490638e-01
|
||||
};
|
||||
return _x;
|
||||
}
|
||||
static std::array<double,31> const & w(){
|
||||
static const std::array<double,31> _w =
|
||||
{
|
||||
5.14947294294515676e-02,
|
||||
5.14261285374590259e-02,
|
||||
5.12215478492587722e-02,
|
||||
5.08817958987496065e-02,
|
||||
5.04059214027823468e-02,
|
||||
4.97956834270742064e-02,
|
||||
4.90554345550297789e-02,
|
||||
4.81858617570871291e-02,
|
||||
4.71855465692991539e-02,
|
||||
4.60592382710069881e-02,
|
||||
4.48148001331626632e-02,
|
||||
4.34525397013560693e-02,
|
||||
4.19698102151642461e-02,
|
||||
4.03745389515359591e-02,
|
||||
3.86789456247275930e-02,
|
||||
3.68823646518212292e-02,
|
||||
3.49793380280600241e-02,
|
||||
3.29814470574837260e-02,
|
||||
3.09072575623877625e-02,
|
||||
2.87540487650412928e-02,
|
||||
2.65099548823331016e-02,
|
||||
2.41911620780806014e-02,
|
||||
2.18280358216091923e-02,
|
||||
1.94141411939423812e-02,
|
||||
1.69208891890532726e-02,
|
||||
1.43697295070458048e-02,
|
||||
1.18230152534963417e-02,
|
||||
9.27327965951776343e-03,
|
||||
6.63070391593129217e-03,
|
||||
3.89046112709988405e-03,
|
||||
1.38901369867700762e-03
|
||||
};
|
||||
return _w;
|
||||
}
|
||||
static std::array<double,15> const & wg(){
|
||||
static const std::array<double,15> _wg =
|
||||
{
|
||||
1.02852652893558840e-01,
|
||||
1.01762389748405505e-01,
|
||||
9.95934205867952671e-02,
|
||||
9.63687371746442596e-02,
|
||||
9.21225222377861287e-02,
|
||||
8.68997872010829798e-02,
|
||||
8.07558952294202154e-02,
|
||||
7.37559747377052063e-02,
|
||||
6.59742298821804951e-02,
|
||||
5.74931562176190665e-02,
|
||||
4.84026728305940529e-02,
|
||||
3.87991925696270496e-02,
|
||||
2.87847078833233693e-02,
|
||||
1.84664683110909591e-02,
|
||||
7.96819249616660562e-03
|
||||
};
|
||||
return _wg;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
} // end of namespace
|
||||
#endif
|
||||
356
libs/glq_integrator.h
Normal file
356
libs/glq_integrator.h
Normal file
|
|
@ -0,0 +1,356 @@
|
|||
/*
|
||||
* Copyright(C) 2017, Andrej Prochazka
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU Affero General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU Affero General Public License for more details.
|
||||
|
||||
* You should have received a copy of the GNU Affero General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef GLQ_INTEGRATOR_H
|
||||
#define GLQ_INTEGRATOR_H
|
||||
#include <array>
|
||||
|
||||
namespace integrators{
|
||||
|
||||
template<int order>
|
||||
struct GL_data{
|
||||
};
|
||||
|
||||
template<int order>
|
||||
class GaussLegendreIntegration{
|
||||
public:
|
||||
template<typename F>
|
||||
double integrate(F& f, double a, double b) const;
|
||||
template<typename F>
|
||||
double operator()(F& f, double a, double b) const {return integrate(f, a, b);}
|
||||
double w(int i) const {return GL_data<order>::w()[i];}
|
||||
double x(int i) const {return GL_data<order>::x()[i];}
|
||||
int n() const {return order;}
|
||||
std::array<double,order> get_points(double a = -1.0, double b = 1.0)const;
|
||||
};
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussLegendreIntegration<order>::integrate(F& f, double a, double b) const{
|
||||
double res=0.0;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
if(order%2){res+= w(0) * (f(p*x(0) + q));} // in case odd-order
|
||||
for(int i=order%2;i<order/2 + order%2;i++){
|
||||
res += w(i) * (f(p*x(i) + q) + f(-p*x(i) + q));
|
||||
}
|
||||
return p*res;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
std::array<double,order> GaussLegendreIntegration<order>::get_points(double a, double b)const{
|
||||
std::array<double,order> points;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
int num = (order/2);
|
||||
for(int i=0;i< num;i++){
|
||||
points[num-i-1] = -p*x(i) + q;
|
||||
points[num+i] = p*x(i) + q;
|
||||
}
|
||||
return points;
|
||||
}
|
||||
|
||||
|
||||
template<int order>
|
||||
class GaussLegendreIntegration2D{
|
||||
public:
|
||||
template<typename F>
|
||||
double integrate(F& f, double a, double b, double c, double d) const;
|
||||
template<typename F>
|
||||
double operator()(F& f, double a, double b, double c, double d) const {return integrate(f, a,b,c,d);}
|
||||
private:
|
||||
GaussLegendreIntegration<order> integrator;
|
||||
};
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussLegendreIntegration2D<order>::integrate(F& f, double a, double b, double c, double d) const{
|
||||
double res=0.0;
|
||||
double p = (b-a)/2.0;
|
||||
double q = (b+a)/2.0;
|
||||
double r = (d-c)/2.0;
|
||||
double s = (d+c)/2.0;
|
||||
|
||||
double xx, yy, sum;
|
||||
for(int i=0;i<order/2;i++){
|
||||
xx = (p*integrator.x(i));
|
||||
for(int j=0;j<order/2;j++){
|
||||
yy = (r*integrator.x(j));
|
||||
sum = (f(xx + q, yy + s) + f(xx + q, -yy+s));
|
||||
sum += (f(-xx + q, yy + s) + f(-xx + q, -yy+s));
|
||||
res += integrator.w(i) * integrator.w(j) * sum;
|
||||
}
|
||||
}
|
||||
return p*r*res;
|
||||
}
|
||||
|
||||
template<int order, int order2 = order>
|
||||
class GaussLegendreIntegration2DA{
|
||||
public:
|
||||
template<typename F>
|
||||
double integrate(F& f, double a, double b, double c, double d) const;
|
||||
template<typename F>
|
||||
double operator()(F& f, double a, double b, double c, double d) const {return integrate(f, a,b,c,d);}
|
||||
private:
|
||||
GaussLegendreIntegration<order> integrator;
|
||||
GaussLegendreIntegration<order2> integrator2;
|
||||
};
|
||||
|
||||
template<int order, int order2>
|
||||
template<typename F>
|
||||
double GaussLegendreIntegration2DA<order, order2>::integrate(F& f, double a, double b, double c, double d) const{
|
||||
double res=0.0;
|
||||
double p = (b-a)/2.0;
|
||||
double q = (b+a)/2.0;
|
||||
double r = (d-c)/2.0;
|
||||
double s = (d+c)/2.0;
|
||||
|
||||
double xx, yy, sum;
|
||||
for(int i=0;i<order/2;i++){
|
||||
xx = (p*integrator.x(i));
|
||||
for(int j=0;j<order2/2;j++){
|
||||
yy = (r*integrator2.x(j));
|
||||
sum = (f(xx + q, yy + s) + f(xx + q, -yy+s));
|
||||
sum += (f(-xx + q, yy + s) + f(-xx + q, -yy+s));
|
||||
res += integrator.w(i) * integrator2.w(j) * sum;
|
||||
}
|
||||
}
|
||||
return p*r*res;
|
||||
}
|
||||
|
||||
/// weights
|
||||
//order = 3
|
||||
|
||||
template<>
|
||||
struct GL_data<3>{
|
||||
static const std::array<double,2>& x(){
|
||||
static const std::array<double,2> _x = {0.0,0.7745966692414833770359};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,2>& w(){
|
||||
static const std::array<double,2> _w = {0.8888888888888888888889,0.555555555555555555556};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 4
|
||||
template<>
|
||||
struct GL_data<4>{
|
||||
static const std::array<double,2>& x(){
|
||||
static const std::array<double,2> _x =
|
||||
{0.3399810435848562648026658,0.8611363115940525752239465};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,2>& w(){
|
||||
static const std::array<double,2> _w =
|
||||
{0.6521451548625461426269361,0.3478548451374538573730639};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 5
|
||||
template<>
|
||||
struct GL_data<5>{
|
||||
static const std::array<double,3>& x(){
|
||||
static const std::array<double,3> _x =
|
||||
{0.0, 0.5384693101056830910363, 0.9061798459386639927976};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,3>& w(){
|
||||
static const std::array<double,3> _w =
|
||||
{0.568888888888888888888, 0.4786286704993664680413, 0.2369268850561890875143};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 6
|
||||
template<>
|
||||
struct GL_data<6>{
|
||||
static const std::array<double,3>& x(){
|
||||
static const std::array<double,3> _x =
|
||||
{0.2386191860831969086305017,0.6612093864662645136613996,0.9324695142031520278123016};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,3>& w(){
|
||||
static const std::array<double,3> _w =
|
||||
{0.4679139345726910473898703,0.3607615730481386075698335,0.1713244923791703450402961};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 7
|
||||
template<>
|
||||
struct GL_data<7>{
|
||||
static const std::array<double,4>& x(){
|
||||
static const std::array<double,4> _x =
|
||||
{0.0, 0.4058451513773971669066, 0.7415311855993944398639,0.9491079123427585245262};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,4>& w(){
|
||||
static const std::array<double,4> _w =
|
||||
{0.417959183673469387755, 0.38183005050511894495, 0.279705391489276667901, 0.129484966168869693271};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 8
|
||||
template<>
|
||||
struct GL_data<8>{
|
||||
static const std::array<double,4>& x(){
|
||||
static const std::array<double,4> _x =
|
||||
{0.1834346424956498049394761,0.5255324099163289858177390,0.7966664774136267395915539,0.9602898564975362316835609};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,4>& w(){
|
||||
static const std::array<double,4> _w =
|
||||
{0.3626837833783619829651504,0.3137066458778872873379622,0.2223810344533744705443560,0.1012285362903762591525314};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 10
|
||||
template<>
|
||||
struct GL_data<10>{
|
||||
static const std::array<double,5>& x(){
|
||||
static const std::array<double,5> _x =
|
||||
{0.1488743389816312108848260,0.4333953941292471907992659,0.6794095682990244062343274,0.8650633666889845107320967,0.9739065285171717200779640};
|
||||
return _x;
|
||||
}
|
||||
static const std::array<double,5>& w(){
|
||||
static const std::array<double,5> _w =
|
||||
{0.2955242247147528701738930,0.2692667193099963550912269,0.2190863625159820439955349,0.1494513491505805931457763,0.0666713443086881375935688};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 12
|
||||
template<>
|
||||
struct GL_data<12>{
|
||||
static const std::array<double,6>& x(){
|
||||
static const std::array<double,6> _x =
|
||||
{0.1252334085114689154724414,0.3678314989981801937526915,0.5873179542866174472967024,0.7699026741943046870368938,0.9041172563704748566784659,0.9815606342467192506905491};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,6>& w(){
|
||||
static const std::array<double,6> _w =
|
||||
{0.2491470458134027850005624,0.2334925365383548087608499,0.2031674267230659217490645,0.1600783285433462263346525,0.1069393259953184309602547,0.0471753363865118271946160};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 14
|
||||
template<>
|
||||
struct GL_data<14>{
|
||||
static const std::array<double,7>& x(){
|
||||
static const std::array<double,7> _x =
|
||||
{0.1080549487073436620662447,0.3191123689278897604356718,0.5152486363581540919652907,0.6872929048116854701480198,0.8272013150697649931897947,0.9284348836635735173363911,0.9862838086968123388415973};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,7>& w(){
|
||||
static const std::array<double,7> _w =
|
||||
{0.2152638534631577901958764,0.2051984637212956039659241,0.1855383974779378137417166,0.1572031671581935345696019,0.1215185706879031846894148,0.0801580871597602098056333,0.0351194603317518630318329};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 16
|
||||
template<>
|
||||
struct GL_data<16>{
|
||||
static const std::array<double,8>& x(){
|
||||
static const std::array<double,8> _x =
|
||||
{0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,8>& w(){
|
||||
static const std::array<double,8> _w =
|
||||
{0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 20
|
||||
template<>
|
||||
struct GL_data<20>{
|
||||
static const std::array<double,10>& x(){
|
||||
static const std::array<double,10> _x =
|
||||
{0.0765265211334973337546404,0.2277858511416450780804962,0.3737060887154195606725482,0.5108670019508270980043641,0.6360536807265150254528367,0.7463319064601507926143051,0.8391169718222188233945291,0.9122344282513259058677524,0.9639719272779137912676661,0.9931285991850949247861224};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,10>& w(){
|
||||
static const std::array<double,10> _w =
|
||||
{0.1527533871307258506980843,0.1491729864726037467878287,0.1420961093183820513292983,0.1316886384491766268984945,0.1181945319615184173123774,0.1019301198172404350367501,0.0832767415767047487247581,0.0626720483341090635695065,0.0406014298003869413310400,0.0176140071391521183118620};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
// order = 32
|
||||
template<>
|
||||
struct GL_data<32>{
|
||||
static const std::array<double,16>& x(){
|
||||
static const std::array<double,16> _x =
|
||||
{0.0483076656877383162348126,0.1444719615827964934851864,0.2392873622521370745446032,0.3318686022821276497799168,0.4213512761306353453641194,0.5068999089322293900237475,0.5877157572407623290407455,0.6630442669302152009751152,0.7321821187402896803874267,0.7944837959679424069630973,0.8493676137325699701336930,0.8963211557660521239653072,0.9349060759377396891709191,0.9647622555875064307738119,0.9856115115452683354001750,0.9972638618494815635449811};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,16>& w(){
|
||||
static const std::array<double,16> _w =
|
||||
{0.0965400885147278005667648,0.0956387200792748594190820,0.0938443990808045656391802,0.0911738786957638847128686,0.0876520930044038111427715,0.0833119242269467552221991,0.0781938957870703064717409,0.0723457941088485062253994,0.0658222227763618468376501,0.0586840934785355471452836,0.0509980592623761761961632,0.0428358980222266806568786,0.0342738629130214331026877,0.0253920653092620594557526,0.0162743947309056706051706,0.0070186100094700966004071};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
//order = 64
|
||||
template<>
|
||||
struct GL_data<64>{
|
||||
static const std::array<double,32>& x(){
|
||||
static const std::array<double,32> _x =
|
||||
{0.0243502926634244325089558,0.0729931217877990394495429,0.1214628192961205544703765,0.1696444204239928180373136,0.2174236437400070841496487,0.2646871622087674163739642,0.3113228719902109561575127,0.3572201583376681159504426,0.4022701579639916036957668,0.4463660172534640879849477,0.4894031457070529574785263,0.5312794640198945456580139,0.5718956462026340342838781,0.6111553551723932502488530,0.6489654712546573398577612,0.6852363130542332425635584,0.7198818501716108268489402,0.7528199072605318966118638,0.7839723589433414076102205,0.8132653151227975597419233,0.8406292962525803627516915,0.8659993981540928197607834,0.8893154459951141058534040,0.9105221370785028057563807,0.9295691721319395758214902,0.9464113748584028160624815,0.9610087996520537189186141,0.9733268277899109637418535,0.9833362538846259569312993,0.9910133714767443207393824,0.9963401167719552793469245,0.9993050417357721394569056};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,32>& w(){
|
||||
static const std::array<double,32> _w =
|
||||
{0.0486909570091397203833654,0.0485754674415034269347991,0.0483447622348029571697695,0.0479993885964583077281262,0.0475401657148303086622822,0.0469681828162100173253263,0.0462847965813144172959532,0.0454916279274181444797710,0.0445905581637565630601347,0.0435837245293234533768279,0.0424735151236535890073398,0.0412625632426235286101563,0.0399537411327203413866569,0.0385501531786156291289625,0.0370551285402400460404151,0.0354722132568823838106931,0.0338051618371416093915655,0.0320579283548515535854675,0.0302346570724024788679741,0.0283396726142594832275113,0.0263774697150546586716918,0.0243527025687108733381776,0.0222701738083832541592983,0.0201348231535302093723403,0.0179517157756973430850453,0.0157260304760247193219660,0.0134630478967186425980608,0.0111681394601311288185905,0.0088467598263639477230309,0.0065044579689783628561174,0.0041470332605624676352875,0.0017832807216964329472961};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
//order = 128
|
||||
template<>
|
||||
struct GL_data<128>{
|
||||
static const std::array<double,64>& x(){
|
||||
static const std::array<double,64> _x =
|
||||
{0.0122236989606157641980521,0.0366637909687334933302153,0.0610819696041395681037870,0.0854636405045154986364980,0.1097942311276437466729747,0.1340591994611877851175753,0.1582440427142249339974755,0.1823343059853371824103826,0.2063155909020792171540580,0.2301735642266599864109866,0.2538939664226943208556180,0.2774626201779044028062316,0.3008654388776772026671541,0.3240884350244133751832523,0.3471177285976355084261628,0.3699395553498590266165917,0.3925402750332674427356482,0.4149063795522750154922739,0.4370245010371041629370429,0.4588814198335521954490891,0.4804640724041720258582757,0.5017595591361444642896063,0.5227551520511754784539479,0.5434383024128103634441936,0.5637966482266180839144308,0.5838180216287630895500389,0.6034904561585486242035732,0.6228021939105849107615396,0.6417416925623075571535249,0.6602976322726460521059468,0.6784589224477192593677557,0.6962147083695143323850866,0.7135543776835874133438599,0.7304675667419088064717369,0.7469441667970619811698824,0.7629743300440947227797691,0.7785484755064119668504941,0.7936572947621932902433329,0.8082917575079136601196422,0.8224431169556438424645942,0.8361029150609068471168753,0.8492629875779689691636001,0.8619154689395484605906323,0.8740527969580317986954180,0.8856677173453972174082924,0.8967532880491581843864474,0.9073028834017568139214859,0.9173101980809605370364836,0.9267692508789478433346245,0.9356743882779163757831268,0.9440202878302201821211114,0.9518019613412643862177963,0.9590147578536999280989185,0.9656543664319652686458290,0.9717168187471365809043384,0.9771984914639073871653744,0.9820961084357185360247656,0.9864067427245862088712355,0.9901278184917343833379303,0.9932571129002129353034372,0.9957927585349811868641612,0.9977332486255140198821574,0.9990774599773758950119878,0.9998248879471319144736081};
|
||||
return _x;
|
||||
}
|
||||
|
||||
static const std::array<double,64>& w(){
|
||||
static const std::array<double,64> _w =
|
||||
{0.0244461801962625182113259,0.0244315690978500450548486,0.0244023556338495820932980,0.0243585572646906258532685,0.0243002001679718653234426,0.0242273192228152481200933,0.0241399579890192849977167,0.0240381686810240526375873,0.0239220121367034556724504,0.0237915577810034006387807,0.0236468835844476151436514,0.0234880760165359131530253,0.0233152299940627601224157,0.0231284488243870278792979,0.0229278441436868469204110,0.0227135358502364613097126,0.0224856520327449668718246,0.0222443288937997651046291,0.0219897106684604914341221,0.0217219495380520753752610,0.0214412055392084601371119,0.0211476464682213485370195,0.0208414477807511491135839,0.0205227924869600694322850,0.0201918710421300411806732,0.0198488812328308622199444,0.0194940280587066028230219,0.0191275236099509454865185,0.0187495869405447086509195,0.0183604439373313432212893,0.0179603271850086859401969,0.0175494758271177046487069,0.0171281354231113768306810,0.0166965578015892045890915,0.0162550009097851870516575,0.0158037286593993468589656,0.0153430107688651440859909,0.0148731226021473142523855,0.0143943450041668461768239,0.0139069641329519852442880,0.0134112712886163323144890,0.0129075627392673472204428,0.0123961395439509229688217,0.0118773073727402795758911,0.0113513763240804166932817,0.0108186607395030762476596,0.0102794790158321571332153,0.0097341534150068058635483,0.0091830098716608743344787,0.0086263777986167497049788,0.0080645898904860579729286,0.0074979819256347286876720,0.0069268925668988135634267,0.0063516631617071887872143,0.0057726375428656985893346,0.0051901618326763302050708,0.0046045842567029551182905,0.0040162549837386423131943,0.0034255260409102157743378,0.0028327514714579910952857,0.0022382884309626187436221,0.0016425030186690295387909,0.0010458126793403487793129,0.0004493809602920903763943};
|
||||
return _w;
|
||||
}
|
||||
};
|
||||
|
||||
} // end of namespace
|
||||
#endif
|
||||
|
|
@ -139,7 +139,7 @@ PYBIND11_MODULE(pycatima,m){
|
|||
.def("molar_mass",py::overload_cast<>(&Material::M, py::const_), "get mass")
|
||||
.def("thickness",py::overload_cast<>(&Material::thickness, py::const_), "get thickness")
|
||||
.def("thickness",py::overload_cast<double>(&Material::thickness), "set thickness")
|
||||
.def("thickness_cm",&Material::thickness_cm,"set thickness in cm unit")
|
||||
.def("thickness_cm",py::overload_cast<double>(&Material::thickness_cm),"set thickness in cm unit")
|
||||
.def("I",py::overload_cast<>(&Material::I, py::const_), "get I")
|
||||
.def("I",py::overload_cast<double>(&Material::I), "set I")
|
||||
.def("__str__",&material_to_string);
|
||||
|
|
@ -281,14 +281,14 @@ PYBIND11_MODULE(pycatima,m){
|
|||
|
||||
m.def("srim_dedx_e",&srim_dedx_e);
|
||||
m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e", py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("calculate",py::overload_cast<Projectile&, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("calculate_layers",py::overload_cast<Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx_from_range",py::overload_cast<Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx_from_range",py::overload_cast<Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx",py::overload_cast<Projectile&, const Material&, const Config&>(&dedx), "dedx",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("range",py::overload_cast<Projectile&, const Material&, const Config&>(&range), "range",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("energy_out",py::overload_cast<Projectile&, const std::vector<double>&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("energy_out",py::overload_cast<Projectile&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("calculate",py::overload_cast<Projectile, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("calculate_layers",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx_from_range",py::overload_cast<const Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx_from_range",py::overload_cast<const Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("dedx",py::overload_cast<const Projectile&, const Material&, const Config&>(&dedx), "dedx",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("range",py::overload_cast<const Projectile&, const Material&, const Config&>(&range), "range",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("energy_out",py::overload_cast<const Projectile&, const std::vector<double>&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("energy_out",py::overload_cast<const Projectile&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("lindhard",&bethek_lindhard);
|
||||
m.def("lindhard_X",&bethek_lindhard_X);
|
||||
m.def("get_material",py::overload_cast<int>(&get_material));
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
#include "structures.h"
|
||||
#include "catima/nucdata.h"
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
|
||||
|
||||
namespace catima{
|
||||
|
|
@ -14,7 +15,7 @@ bool operator==(const Projectile &a, const Projectile&b){
|
|||
}
|
||||
|
||||
bool operator==(const Material &a, const Material&b){
|
||||
if(a.density() != b.density())return false;
|
||||
if(std::fabs(a.density() - b.density())> 1e-6)return false;
|
||||
if(a.ncomponents() != b.ncomponents())return false;
|
||||
if(a.I() != b.I())return false;
|
||||
for(int i=0;i<a.ncomponents();i++){
|
||||
|
|
@ -66,6 +67,22 @@ void Layers::add(Material m){
|
|||
materials.push_back(m);
|
||||
}
|
||||
|
||||
double Layers::thickness() const {
|
||||
double sum = 0;
|
||||
for(auto &m : materials){
|
||||
sum += m.thickness();
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
double Layers::thickness_cm() const {
|
||||
double sum = 0;
|
||||
for(auto &m : materials){
|
||||
sum += m.thickness_cm();
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
Layers operator+(const Layers &a, const Layers&b){
|
||||
Layers res;
|
||||
for(auto &e:a.materials){
|
||||
|
|
|
|||
26
structures.h
26
structures.h
|
|
@ -69,7 +69,7 @@ namespace catima{
|
|||
*/
|
||||
class Material{
|
||||
private:
|
||||
double rho=0;
|
||||
double rho=1e-5;
|
||||
double th=0;
|
||||
double molar_mass=0;
|
||||
double i_potential=0;
|
||||
|
|
@ -84,7 +84,7 @@ namespace catima{
|
|||
* @param _rho - density of the material in g/cm3, default 0.0
|
||||
* @param _th - thickness of the material in g/cm2, default 0.0
|
||||
*/
|
||||
Material(double _a, int _z, double _rho=0.0, double _th=0.0, double _ipot = 0.0);
|
||||
Material(double _a, int _z, double _rho=1e-5, double _th=0.0, double _ipot = 0.0);
|
||||
|
||||
|
||||
/**
|
||||
|
|
@ -97,7 +97,7 @@ namespace catima{
|
|||
* });
|
||||
* \endcode
|
||||
*/
|
||||
Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0, double mass=0.0);
|
||||
Material(std::initializer_list<std::array<double,3>>list,double _density=1e-5, double ipot = 0.0, double mass=0.0);
|
||||
|
||||
/**
|
||||
* calculates internal variables if needed
|
||||
|
|
@ -161,6 +161,11 @@ namespace catima{
|
|||
*/
|
||||
double thickness() const {return th;};
|
||||
|
||||
/**
|
||||
* @return returns thickness in cm
|
||||
*/
|
||||
double thickness_cm() const {return th/rho;};
|
||||
|
||||
/**
|
||||
* sets thickness in g/cm^2
|
||||
*/
|
||||
|
|
@ -212,9 +217,6 @@ namespace catima{
|
|||
return number_density_cm2()*molar_fraction(i);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
friend bool operator==(const Material &a, const Material&b);
|
||||
};
|
||||
|
||||
|
|
@ -233,6 +235,8 @@ namespace catima{
|
|||
double sigma_E=0.0;
|
||||
double sigma_a=0.0;
|
||||
double sigma_r=0.0;
|
||||
double sigma_x=0.0;
|
||||
double cov = 0.0;
|
||||
double tof=0.0;
|
||||
#ifdef REACTIONS
|
||||
double sp = 1.0;
|
||||
|
|
@ -275,6 +279,16 @@ namespace catima{
|
|||
*/
|
||||
int num()const{return materials.size();};
|
||||
|
||||
/**
|
||||
* @ return total thickness
|
||||
*/
|
||||
double thickness() const;
|
||||
|
||||
/**
|
||||
* @ return total thickness in cm
|
||||
*/
|
||||
double thickness_cm() const;
|
||||
|
||||
Material& operator[](int i){return materials[i];}
|
||||
|
||||
friend Layers operator+(const Layers &a, const Layers&b);
|
||||
|
|
|
|||
|
|
@ -138,6 +138,7 @@ using namespace std;
|
|||
catima::Material graphite({
|
||||
{12.011,6,1},
|
||||
});
|
||||
graphite.density(2);
|
||||
|
||||
auto res = catima::calculate(p(1000),graphite);
|
||||
CHECK(catima::dedx(p(1000), graphite) == approx(res.dEdxi).R(0.001) );
|
||||
|
|
@ -183,6 +184,93 @@ using namespace std;
|
|||
dif = res.tof - 0.038;
|
||||
CHECK( fabs(dif) < 0.01);
|
||||
}
|
||||
TEST_CASE("angular scattering"){
|
||||
catima::Projectile p{1,1,1,158.6};
|
||||
catima::Material cu = catima::get_material(29);
|
||||
|
||||
|
||||
cu.thickness_cm(0.02963);
|
||||
auto res = catima::calculate(p,cu);
|
||||
CHECK( 1000*res.sigma_a == approx(7.2,0.5));
|
||||
|
||||
cu.thickness_cm(0.2963);
|
||||
res = catima::calculate(p,cu);
|
||||
CHECK( 1000*res.sigma_a == approx(23.5,3));
|
||||
|
||||
catima::Layers ll;
|
||||
cu.thickness_cm(0.2963/2.);
|
||||
ll.add(cu);
|
||||
ll.add(cu);
|
||||
auto res2 = catima::calculate(p,ll);
|
||||
CHECK( 1000*res2.total_result.sigma_a == approx(23.5,3));
|
||||
CHECK( 1000*res2.total_result.sigma_a == approx(1000*res.sigma_a,0.05));
|
||||
}
|
||||
|
||||
TEST_CASE("displacement test"){
|
||||
catima::Projectile p{1,1,1,215};
|
||||
catima::Material water({
|
||||
{1.00794,1,2},
|
||||
{15.9994,8,1}
|
||||
});
|
||||
water.density(1.0);
|
||||
catima::Material water2({
|
||||
{1.00794,1,2},
|
||||
{15.9994,8,1}
|
||||
});
|
||||
water2.thickness_cm(9.6);
|
||||
water2.density(2.0);
|
||||
|
||||
water.thickness_cm(9.6);
|
||||
auto res = catima::calculate(p,water);
|
||||
CHECK( res.sigma_x == approx(0.1,0.03));
|
||||
auto resb = catima::calculate(p,water2);
|
||||
CHECK( res.sigma_x > resb.sigma_x);
|
||||
|
||||
water.thickness_cm(17.0);
|
||||
auto res17 = catima::calculate(p,water);
|
||||
CHECK( res17.sigma_x == approx(0.25,0.05));
|
||||
|
||||
water.thickness_cm(29.4);
|
||||
auto res29 = catima::calculate(p,water);
|
||||
CHECK( res29.sigma_x == approx(0.66,0.08));
|
||||
|
||||
catima::Layers ll;
|
||||
water.thickness_cm(9.6);
|
||||
ll.add(water);
|
||||
auto res2 = catima::calculate(p,ll);
|
||||
CHECK(res2.total_result.sigma_x == approx(0.1,0.03));
|
||||
water.thickness_cm(7.4);
|
||||
ll.add(water);
|
||||
res2 = catima::calculate(p,ll);
|
||||
CHECK(ll.thickness_cm() == approx(17,0.00001));
|
||||
CHECK(res2.total_result.sigma_x == approx(0.25,0.05));
|
||||
CHECK(res2.total_result.sigma_x == approx(res17.sigma_x).R(0.03));
|
||||
|
||||
catima::Layers l;
|
||||
water.thickness_cm(9.6/3);
|
||||
l.add(water);
|
||||
l.add(water);
|
||||
l.add(water);
|
||||
res2 = catima::calculate(p(215),l);
|
||||
CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.03));
|
||||
|
||||
for(int n=2;n<10;n++){
|
||||
catima::Layers lll;
|
||||
water.thickness_cm(29.4/n);
|
||||
for(int i=0;i<n;i++)lll.add(water);
|
||||
res2 = catima::calculate(p(215),lll);
|
||||
CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
|
||||
}
|
||||
|
||||
|
||||
/// position sigma is taken from range thickness if particle stoppped inside
|
||||
water.thickness_cm(30);
|
||||
auto r1 = catima::calculate(p,water);
|
||||
water.thickness_cm(40);
|
||||
auto r2 = catima::calculate(p,water);
|
||||
CHECK(r1.sigma_x == approx(r2.sigma_x).R(0.01));
|
||||
|
||||
}
|
||||
TEST_CASE("result from stopped material"){
|
||||
catima::Projectile p{12,6,6,1000};
|
||||
catima::Material water({
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ using namespace std;
|
|||
graphite.add_element(18,40,1);
|
||||
CHECK(graphite.ncomponents()==2);
|
||||
CHECK_FALSE(graphite.M()==approx(12,0.1));
|
||||
|
||||
CHECK( water.density() == approx(water2.density(),1e-6) );
|
||||
CHECK(water==water2);
|
||||
CHECK(!(water==graphite));
|
||||
water.density(1.0);
|
||||
|
|
@ -42,8 +42,15 @@ using namespace std;
|
|||
CHECK(water.thickness()==approx(1.0,0.0001));
|
||||
water.thickness_cm(1.0);
|
||||
CHECK(water.thickness()==approx(1.0,0.0001));
|
||||
CHECK(water.thickness_cm()==approx(1.0,0.0001));
|
||||
water.thickness_cm(2.0);
|
||||
CHECK(water.thickness()==approx(2.0,0.0001));
|
||||
CHECK(water.thickness_cm()==approx(2.0,0.0001));
|
||||
water.density(2.0);
|
||||
CHECK(water.thickness()==approx(2.0,0.0001));
|
||||
CHECK(water.thickness_cm()==approx(1.0,0.0001));
|
||||
water.thickness(1.0);
|
||||
CHECK(water.thickness_cm()==approx(0.5,0.0001));
|
||||
}
|
||||
TEST_CASE("Material automatic atomic weight"){
|
||||
catima::Material water({{0,1,2},{0,8,1}});
|
||||
|
|
@ -85,6 +92,7 @@ using namespace std;
|
|||
detector1.add(water2);
|
||||
detector1.add(graphite);
|
||||
CHECK(detector1.num() == 3);
|
||||
CHECK(detector1.thickness() == approx(4.0,1e-6));
|
||||
// check correct density and thickness
|
||||
CHECK(detector1[0].density()==1.8);
|
||||
CHECK(detector1[0].thickness()==1.0);
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user