gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-22 18:28:51 -05:00
Code Issues Projects Releases Wiki Activity

Config refactoring

Browse Source
This commit is contained in:
hrocho 2019-05-10 21:02:51 +02:00
parent 75b6d15fec
commit 2068c037b2
7 changed files with 64 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
calculations.cpp
View File

@ -88,7 +88,7 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
double f2 = log(2.0*electron_mass*1000000*beta2/Ipot);
double eta = beta*gamma;
if(!(c.dedx&corrections::no_shell_correction) && eta>=0.13){ //shell corrections
if(!(c.corrections&corrections::no_shell_correction) && eta>=0.13){ //shell corrections
double cor = (+0.422377*pow(eta,-2)
+0.0304043*pow(eta,-4)
-0.00038106*pow(eta,-6))*1e-6*pow(Ipot,2)
@ -100,21 +100,21 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
f2+=2*log(gamma) -beta2;
double barkas=1.0;
if(!(c.dedx&corrections::no_barkas)){
if(!(c.corrections&corrections::no_barkas)){
barkas = bethek_barkas(zp_eff,eta,t.Z);
}
double delta = bethek_density_effect(beta, t.Z);
double LS = 0.0;
if(!(c.dedx&corrections::no_lindhard)){
if(!(c.corrections&corrections::no_lindhard)){
//double LS = bethek_lindhard(p);
LS = precalculated_lindhard(p);
}
double result = (f2)*barkas + LS - delta/2.;
result *=f1;
if( (p.T>50000.0) && !(c.dedx&corrections::no_highenergy)){
if( (p.T>50000.0) && !(c.corrections&corrections::no_highenergy)){
result += pair_production(p,t);
result += bremsstrahlung(p,t);
}
@ -444,23 +444,16 @@ double bremsstrahlung(const Projectile &p, const Target &t){
return 16.0*C*gamma*p.Z*p.Z*p.Z*p.Z*t.Z*t.Z*Lbs/(t.A*p.A*3.0*4.0*M_PI);
};
double sezi_p_se(double energy,const Target &t){
return 100*p_se(t.Z, energy)*Avogadro/t.A;
}
double sezi_dedx_e(const Projectile &p, const Target &t){
return 100*srim_dedx_e(p.Z,t.Z,p.T)*Avogadro/t.A;
}
double sezi_dedx_e(const Projectile &p, const Material &mat){
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){
double w;
double sum=0.0;
bool use95 = config_lowenergy(c) == low_energy_types::srim_95;
for(int i=0;i<mat.ncomponents();i++){
auto t = mat.get_element(i);
w = mat.weight_fraction(i);
sum += w*sezi_dedx_e(p,t);
sum += w*srim_dedx_e(p.Z,t.Z,p.T, use95)/t.A;
}
return sum;
return 100*sum*Avogadro; // returning MeV/g/cm2
}
@ -583,7 +576,7 @@ double dedx_variance(Projectile &p, Target &t, const Config &c){
double zp_eff = z_effective(p,t,c);
double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
if(c.dedx_straggling == omega::atima){
if(config_omega(c) == omega_types::atima){
cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364)));
cor = std::max(cor, 0.0 );

14
calculations.h
View File

@ -87,21 +87,11 @@ namespace catima{
* this function is not used and is not tested
*/
double energy_straggling_firsov(double z1,double energy, double z2, double m2);
/**
* electronic energy loss for low energy, should be like SRIM
*/
double sezi_dedx_e(const Projectile &p, const Target &t);
/**
* electronic energy loss for low energy, should be like SRIM
*/
double sezi_dedx_e(const Projectile &p, const Material &mat);
/**
* electronic energy loss of protons for low energy, should be like SRIM
*/
double sezi_p_se(double energy,const Target &t);
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
double angular_scattering_variance(Projectile &p, Target &t);

42
catima.pyx
View File

@ -305,47 +305,53 @@ cdef class Config:
#self.cbase = catimac.Config()
self.cbase.z_effective = z_eff_type.pierce_blann
self.cbase.skip = 0
self.cbase.dedx = 0
self.cbase.dedx_straggling = omega_type.atima
self.cbase.calculation = 1
self.cbase.corrections = 0
def z_effective(self, val=None):
if(val is None):
return self.cbase.z_effective
else:
self.cbase.z_effective = val
def skip_calculation(self, val=None):
if(val is None):
return self.cbase.skip
else:
self.cbase.skip = val
def dedx(self, val=None):
def corrections(self, val=None):
if(val is None):
return self.cbase.dedx
return self.cbase.corrections
else:
self.cbase.dedx = val
def dedx_straggling(self, val=None):
self.cbase.corrections = val
def calculation(self, val=None):
if(val is None):
return self.cbase.dedx_straggling
return self.cbase.calculation
else:
self.cbase.dedx_straggling = val
self.cbase.calculation = val
def set(self,other):
if("z_effective" in other):
self.cbase.z_effective = other["z_effective"]
if("dedx" in other):
self.cbase.dedx = other["dedx"]
if("dedx_straggling" in other):
self.cbase.dedx_straggling = other["dedx_straggling"]
if("calculation" in other):
self.cbase.calculation = other["calculation"]
if("corrections" in other):
self.cbase.corrections = other["corrections"]
def get(self):
res = {}
res["z_effective"] = self.cbase.z_effective
res["dedx"] = self.cbase.dedx
res["dedx_straggling"] = self.cbase.dedx_straggling
res["corrections"] = self.cbase.corrections
res["calculation"] = self.cbase.calculation
res["skip"] = self.cbase.skip
return res
def print_info(self):
print("z_effective = %s"%z_eff_type(self.cbase.z_effective))
print("dedx_straggling = %s"%omega_type(self.cbase.dedx_straggling))
print("calculation = %s"%omega_type(self.cbase.dedx_straggling))
default_config = Config()
@ -464,9 +470,6 @@ def w_magnification(Projectile projectile, Material material, energy = None, Con
energy = projectile.T()
return catimac.w_magnification(projectile.cbase,energy, material.cbase, config.cbase)
def sezi_dedx_e(Projectile projectile, Target t):
return catimac.sezi_dedx_e(projectile.cbase, t.cbase)
def bethek_dedx_e(Projectile projectile, Target t, Config c = default_config, Ipot=0.0):
return catimac.bethek_dedx_e(projectile.cbase, t.cbase,c.cbase,Ipot)
@ -506,9 +509,6 @@ def gamma_from_T(double T):
def beta_from_T(double T):
return catimac.beta_from_T(T);
def sezi_p_se(double T, Target t):
return catimac.sezi_p_se(T, t.cbase)
def get_data(Projectile projectile, Material material, Config config = default_config):
data = catimac.get_data(projectile.cbase, material.cbase, config.cbase)
return [data.range,data.range_straggling,data.angular_variance]

6
catimac.pxd
View File

@ -68,8 +68,8 @@ cdef extern from "catima/config.h" namespace "catima":
cdef struct Config:
char z_effective;
char skip;
char dedx;
char dedx_straggling
char corrections;
char calculation;
cdef extern from "catima/catima.h" namespace "catima":
cdef double dedx(Projectile &p, double T, const Material &t,const Config &c)
@ -96,7 +96,6 @@ cdef extern from "catima/calculations.h" namespace "catima":
cdef double bethek_lindhard(const Projectile &p);
cdef double bethek_lindhard_X(const Projectile &p);
cdef double bethek_dedx_e(Projectile &p,const Target &t, const Config &c, double I);
cdef double sezi_dedx_e(const Projectile &p, const Target &t);
cdef double z_effective(const Projectile &p, const Target &t, const Config &c);
cdef double z_eff_Pierce_Blann(double z, double beta);
cdef double z_eff_Anthony_Landford(double pz, double beta, double tz);
@ -107,7 +106,6 @@ cdef extern from "catima/calculations.h" namespace "catima":
cdef double z_eff_Schiwietz(double pz, double beta, double tz);
cdef double gamma_from_T(double T);
cdef double beta_from_T(double T);
cdef double sezi_p_se(double, const Target&);
cdef extern from "catima/constants.h" namespace "catima":
int max_datapoints "catima::max_datapoints"

23
srim.cpp
View File

@ -8,7 +8,7 @@ namespace catima{
* @param energy - energy per nuclein in MeV/u
*/
double p_se95(int Z, double energy){
double p_se(int Z, double energy){
double sp = -1;
double e = 1000*energy; //e in keV/u
int i = Z - 1;
@ -16,8 +16,8 @@ double p_se95(int Z, double energy){
i = 91;
}
if(e<=25)e=25;
double sl = (pse[i][0]*std::pow(e,pse[i][1])) + (pse[i][2]*std::pow(e,pse[i][3]));
double sh = pse[i][4]/std::pow(e,pse[i][5]) * std::log( (pse[i][6]/e) + (pse[i][7]*e));
double sl = (pse_95[i][0]*std::pow(e,pse_95[i][1])) + (pse_95[i][2]*std::pow(e,pse_95[i][3]));
double sh = pse_95[i][4]/std::pow(e,pse_95[i][5]) * std::log( (pse_95[i][6]/e) + (pse_95[i][7]*e));
sp = sl*sh/(sl+sh);
e=1000*energy;
if(e<=25){
@ -31,7 +31,7 @@ double p_se95(int Z, double energy){
* @param Z - proton number of material
* @param energy - energy per nuclein in MeV/u
*/
double p_se(int Z, double energy){
double p_se85(int Z, double energy){
double sp = -1;
double e = 1000*energy; //e in keV/u
int i = Z - 1;
@ -56,12 +56,13 @@ double p_se(int Z, double energy){
* @param pZ - material Z
* @param energy - projectile energy in MeV/u unit
*/
double srim_dedx_e(int pZ, int tZ, double energy){
double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
double e=energy*1000; // e in keV/u
double se = 0;
double (*const fp_se)(int, double) = (use_new)?p_se:p_se85;
// double (*const fp_se)(int, double) = p_se85;
if(pZ==1){
return p_se(tZ, energy);
return (*fp_se)(tZ, energy);
}
else if(pZ == 2){
double a=0;
@ -76,7 +77,7 @@ double srim_dedx_e(int pZ, int tZ, double energy){
a = (1.0 + (0.007 + 0.00005*tZ)*exp(- b*b ));
heh *= a*a;
//zeta = sqrt(heh);
se = p_se(tZ, energy)*heh*4.0; //scale proton stopping
se = (*fp_se)(tZ, energy)*heh*4.0; //scale proton stopping
if(e==1)se*= sqrt(e); //vel proportional
return se;
}
@ -150,17 +151,17 @@ double srim_dedx_e(int pZ, int tZ, double energy){
h1 = zeta *pZ;
h4 = std::pow(e / eee,eval);
se = p_se(tZ, eee*0.001) * h1*h1*h4;
se = (*fp_se)(tZ, eee*0.001) * h1*h1*h4;
return se;
}
else {
return p_se(tZ,energy)*std::pow(zeta*pZ,2.0);
return (*fp_se)(tZ,energy)*std::pow(zeta*pZ,2.0);
}
return 0;
}
};
const double pse[92][8] = {
const double pse_95[92][8] = {
//H
{0.0128116,0.00533047,0.651042,0.531902,1959.01,1.1887,598.263,0.00954514},
//He

5
srim.h
View File

@ -6,8 +6,9 @@ namespace catima{
* @param pZ - projectile Z
* @param pZ - material Z
* @param energy - projectile energy in MeV/u unit
* @param use_v95 - use srim proton coefficient from version 95, otherwise version 85 will be used
*/
double srim_dedx_e(int pZ, int tZ, double energy);
double srim_dedx_e(int pZ, int tZ, double energy, bool use_v95=1);
/**
* return SRIM proton stopping power
@ -17,7 +18,7 @@ double srim_dedx_e(int pZ, int tZ, double energy);
double p_se(int Z, double energy);
double p_se95(int Z, double energy);
extern const double pse[92][8];
extern const double pse_95[92][8];
extern const double atima_lambda_screening[92];
extern const double atima_vfermi[92];
extern const double proton_stopping_coef[92][8];

39
tests/test_calculations.cpp
View File

@ -34,36 +34,27 @@ const lest::test specification[] =
CASE("proton stopping power from srim"){
catima::Projectile p{1,1,1,1};
catima::Target he{4.002600,2};
catima::Target carbon{12.0107,6};
EXPECT( catima::sezi_p_se(1,he) == approx(283,1));
p.T = 1;
EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
EXPECT(catima::sezi_p_se(10,he)==approx(45.6,1));
p.T = 10;
EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
EXPECT(catima::sezi_p_se(30,he) == approx(18.38,1));
p.T = 30;
EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
EXPECT( catima::sezi_p_se(1,carbon) == approx(229.5,1));
p.T = 1;
EXPECT( catima::sezi_p_se(p.T,carbon) == approx(catima::sezi_dedx_e(p,carbon)).R(1e-6));
auto he = catima::get_material(2);
auto carbon = catima::get_material(6);
p.T = 1;
EXPECT( catima::sezi_dedx_e(p,he) == approx(283,1));
p.T = 10;
EXPECT( catima::sezi_dedx_e(p,he) == approx(45.6,1));
EXPECT( catima::sezi_p_se(10,carbon) == approx(40.8,1));
EXPECT(catima::sezi_p_se(30,carbon) == approx(16.8,1));
p.T = 30;
EXPECT( catima::sezi_dedx_e(p,he) == approx(18.38,1));
p.T = 1;
EXPECT( catima::sezi_dedx_e(p,carbon) == approx(229.5,1));
p.T = 10;
EXPECT( catima::sezi_dedx_e(p,carbon) == approx(40.8,1));
p.T = 30;
EXPECT( catima::sezi_p_se(p.T,carbon) == approx(catima::sezi_dedx_e(p,carbon)).R(1e-6));
EXPECT( catima::sezi_dedx_e(p,carbon) == approx(16.8,1));
},
CASE("dedx, low energy, from sezi"){
catima::Projectile p{4,2,2,1};
catima::Target carbon{12.0107,6};
auto carbon = catima::get_material(6);
// He projectile case
p.T = 1;