gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-26 12:08:52 -05:00
Code Issues Projects Releases Wiki Activity

Merge pull request #72 from hrosiak/u3

Browse Source 
u3
This commit is contained in:
Andrej Prochazka 2019-11-21 14:00:43 +01:00 committed by GitHub
parent fed6d9f57d a727b3f93d
commit 8e062b0bf9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 93 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
calculations.cpp
View File

@ -532,10 +532,6 @@ double precalculated_lindhard(const Projectile &p){
double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z); double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z);
z = z-1; z = z-1;
//catima::Interpolator ls_a(ls_coefficients::ls_energy_points,ls_coefficients::ls_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
//catima::Interpolator ls_ahi(ls_coefficients::ls_energy_points,ls_coefficients::ls_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
//catima::Interpolator ls_a(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::cspline);
//catima::Interpolator ls_ahi(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::cspline);
double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_a[z]); double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_a[z]);
double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_ahi[z]); double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_ahi[z]);
@ -551,13 +547,10 @@ double precalculated_lindhard_X(const Projectile &p){
int z = (int)p.Z ; int z = (int)p.Z ;
if(z>LS_MAX_Z)z=LS_MAX_Z; if(z>LS_MAX_Z)z=LS_MAX_Z;
//if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0); //if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0);
if(p.T<ls_coefficients::ls_energy_table(0)) if(p.T<ls_coefficients::ls_energy_table(0))return 1.0;
return 1.0;
double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z); double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z);
z = z-1; z = z-1;
//catima::Interpolator ls_X_a(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_X_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
//catima::Interpolator ls_X_ahi(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_X_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_a[z]); double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_a[z]);
double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_ahi[z]); double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_ahi[z]);

1
catima.cpp
View File

@ -371,7 +371,6 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
return 1.0/dedx(p(x),t,c); return 1.0/dedx(p(x),t,c);
}; };
auto fomega = [&](double x)->double{ auto fomega = [&](double x)->double{
//return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t,c),3);
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);
}; };

10
pymodule/pycatima.cpp
View File

@ -120,6 +120,13 @@ PYBIND11_MODULE(pycatima,m){
//.def_readwrite("T", &Projectile::T) //.def_readwrite("T", &Projectile::T)
//.def_readwrite("Q", &Projectile::Q); //.def_readwrite("Q", &Projectile::Q);
py::class_<Target>(m,"Target")
.def(py::init<>(),"constructor")
.def_readwrite("A",&Target::A)
.def_readwrite("Z",&Target::Z)
.def_readwrite("stn",&Target::stn);
py::class_<Material>(m,"Material") py::class_<Material>(m,"Material")
.def(py::init<>(),"constructor") .def(py::init<>(),"constructor")
.def(py::init<const Material&>(),"constructor") .def(py::init<const Material&>(),"constructor")
@ -210,7 +217,7 @@ PYBIND11_MODULE(pycatima,m){
.value("hubert", z_eff_type::hubert) .value("hubert", z_eff_type::hubert)
.value("winger", z_eff_type::winger) .value("winger", z_eff_type::winger)
.value("schiwietz", z_eff_type::schiwietz) .value("schiwietz", z_eff_type::schiwietz)
.value("global", z_eff_type::global) .value("cglobal", z_eff_type::global)
.value("atima14", z_eff_type::atima14); .value("atima14", z_eff_type::atima14);
py::enum_<corrections>(m,"corrections") py::enum_<corrections>(m,"corrections")
@ -295,6 +302,7 @@ PYBIND11_MODULE(pycatima,m){
m.def("storage_info",&storage_info); m.def("storage_info",&storage_info);
m.def("get_energy_table",&get_energy_table); m.def("get_energy_table",&get_energy_table);
m.def("energy_table",[](int i){return energy_table(i);}); m.def("energy_table",[](int i){return energy_table(i);});
m.def("z_effective",&z_effective);
m.attr("max_datapoints") = max_datapoints; m.attr("max_datapoints") = max_datapoints;
m.attr("max_storage_data") = max_storage_data; m.attr("max_storage_data") = max_storage_data;
m.attr("logEmin")=logEmin; m.attr("logEmin")=logEmin;

29
srim.cpp
View File

@ -82,7 +82,6 @@ double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
return se; return se;
} }
else{ // heavy ion else{ // heavy ion
double h1;
double a,q,b; double a,q,b;
double l1,l0,l; double l1,l0,l;
double YRmin = 0.130; // YRmin = VR / ZP**0.67 <= 0.13 OR VR <= 1.0 double YRmin = 0.130; // YRmin = VR / ZP**0.67 <= 0.13 OR VR <= 1.0
@ -93,6 +92,8 @@ double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
double zeta = 0; double zeta = 0;
double se; double se;
int i; int i;
double pZ13_inv = 1.0/std::pow(pZ,1.0/3.0);
double pZ23_inv = pZ13_inv*pZ13_inv;
i = tZ - 1; i = tZ - 1;
if(tZ>92){ if(tZ>92){
@ -103,21 +104,18 @@ double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
v = sqrt(e/25.0)/vfermi; v = sqrt(e/25.0)/vfermi;
double v2=v*v; double v2=v*v;
double vr = (v >= 1)? v*vfermi*(1.+ 1./(5.*v2)) : 3.0*vfermi/4.0*(1.0+v2*(2.0/3.0-v2/15.0)); double vr = (v >= 1)? v*vfermi*(1.+ 0.2/v2) : 3.0*vfermi/4.0*(1.0+v2*(2.0/3.0-v2/15.0));
h1= 1./std::pow(pZ,0.6667); yr = std::max(YRmin,vr * pZ23_inv);
yr = std::max(YRmin,vr*h1); yr = std::max(yr,VRmin * pZ23_inv);
yr = std::max(yr, VRmin*h1);
//-- CALCULATE ZEFF //-- CALCULATE ZEFF
a = -0.803*std::pow(yr,0.3) + 1.3167*std::pow(yr,0.6) + 0.38157*yr + 0.008983*yr*yr; a = -0.803*std::pow(yr,0.3) + 1.3167*std::pow(yr,0.6) + 0.38157*yr + 0.008983*yr*yr;
q = std::min(1.0, std::max(0.0 , (1.0 - exp(-std::min(a, 50.0))))); //-- Q = IONIZATION LEVEL OF THE ION AT RELATIVE VELOCITY YR q = std::min(1.0, std::max(0.0 , (1.0 - exp(-std::min(a, 50.0))))); //-- Q = IONIZATION LEVEL OF THE ION AT RELATIVE VELOCITY YR
//-- IONIZATION LEVEL TO EFFECTIVE CHARGE //-- IONIZATION LEVEL TO EFFECTIVE CHARGE
h1 = 1./ std::pow(pZ,0.3333); b = (std::min(0.43, std::max(0.32,0.12 + 0.025*pZ)))*pZ13_inv;
l0 = (.8 - q * std::min(1.2,0.6 +pZ/30.0))*pZ13_inv;
b = (std::min(0.43, std::max(0.32,0.12 + 0.025*pZ)))*h1;
l0 = (.8 - q * std::min(1.2,0.6 +pZ/30.0))*h1;
if(q < 0.2){ if(q < 0.2){
l1 = 0; l1 = 0;
} }
@ -133,14 +131,14 @@ double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
// calculate screening // calculate screening
i = (pZ>92)?91:pZ-1; i = (pZ>92)?91:pZ-1;
l = std::max(l1,l0*atima_lambda_screening[i]); l = std::max(l1,l0*atima_lambda_screening[i]);
h1 =4.0*l*vfermi/1.919; double h1 =4.0*l*vfermi/1.919;
zeta = q + (1./(2.*(vfermi*vfermi)))*(1. - q)* log1p(h1*h1); zeta = q + (1./(2.*(vfermi*vfermi)))*(1. - q)* log1p(h1*h1);
// ZP**3 EFFECT AS IN REF. 779? // ZP**3 EFFECT AS IN REF. 779?
a = 7.6 - std::max(0.0, log(e)); a = 7.6 - std::max(0.0, log(e));
zeta = zeta*(1. + (1./(pZ*pZ))*(0.18 + .0015*tZ)*exp(-a*a)); zeta = zeta*(1. + (1./(pZ*pZ))*(0.18 + .0015*tZ)*exp(-a*a));
h1= 1./std::pow(pZ,0.6667); double c1 = zeta *pZ;
if (yr <= ( std::max(YRmin, VRmin*h1))){ if (yr <= ( std::max(YRmin, VRmin * pZ23_inv))){
VRmin=std::max(VRmin, YRmin/h1); VRmin=std::max(VRmin, YRmin/pZ23_inv);
//--C ..CALCULATE VELOCITY STOPPING FOR YR < YRmin //--C ..CALCULATE VELOCITY STOPPING FOR YR < YRmin
double vmin =.5*(VRmin + sqrt(std::max(0.0,VRmin*VRmin - .8*vfermi*vfermi))); double vmin =.5*(VRmin + sqrt(std::max(0.0,VRmin*VRmin - .8*vfermi*vfermi)));
double eee = 25.0*vmin*vmin; double eee = 25.0*vmin*vmin;
@ -149,13 +147,12 @@ double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
if((tZ == 6) || (((tZ == 14) || (tZ == 32)) && (pZ <= 19))) eval = 0.35; if((tZ == 6) || (((tZ == 14) || (tZ == 32)) && (pZ <= 19))) eval = 0.35;
else eval = 0.5; else eval = 0.5;
h1 = zeta *pZ;
double h4 = std::pow(e / eee,eval); double h4 = std::pow(e / eee,eval);
se = (*fp_se)(tZ, eee*0.001) * h1*h1*h4; se = (*fp_se)(tZ, eee*0.001) * c1*c1*h4;
return se; return se;
} }
else { else {
return (*fp_se)(tZ,energy)*std::pow(zeta*pZ,2.0); return (*fp_se)(tZ,energy)*std::pow(c1,2.0);
} }
return 0; return 0;
} }