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

Merge pull request #18 from hrosiak/omega

Browse Source 
Omega
This commit is contained in:
Andrej Prochazka 2018-01-12 01:35:58 +01:00 committed by GitHub
parent d9fc231883 397473360f
commit bb133397c2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 27 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
calculations.cpp
View File

@ -532,8 +532,6 @@ double beta_from_T(double T){
return sqrt(1.0-1.0/(gamma*gamma)); return sqrt(1.0-1.0/(gamma*gamma));
} }
double energy_straggling_firsov(double z1,double energy, double z2, double m2){ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
double gamma = gamma_from_T(energy); double gamma = gamma_from_T(energy);
double beta2=1.0-1.0/(gamma*gamma); double beta2=1.0-1.0/(gamma*gamma);
@ -615,8 +613,9 @@ double precalculated_lindhard_X(const Projectile &p){
double T = p.T; double T = p.T;
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_points[0])T=ls_coefficients::ls_energy_points[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))T=ls_coefficients::ls_energy_table(0); if(p.T<ls_coefficients::ls_energy_table(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;
@ -634,10 +633,20 @@ double precalculated_lindhard_X(const Projectile &p){
double dedx_variance(Projectile &p, Target &t, const Config &c){ double dedx_variance(Projectile &p, Target &t, const Config &c){
double zp_eff = z_effective(p,t,c); double zp_eff = z_effective(p,t,c);
double gamma = gamma_from_T(p.T); double gamma = gamma_from_T(p.T);
double f = domega2dx_constant*pow(zp_eff,2)*t.Z*gamma*gamma; double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
double cor;
double beta2 = pow(beta_from_T(p.T),2);
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 );
//double X = bethek_lindhard_X(p); //double X = bethek_lindhard_X(p);
double X = precalculated_lindhard_X(p); double X = precalculated_lindhard_X(p);
return f*X/t.A; X *= gamma*gamma;
if(p.T<1.0)
return std::min(f*(X+cor), energy_straggling_firsov(p.Z, p.T, t.Z,t.A));
else
return f*(X+cor);
} }
double z_effective(const Projectile &p,const Target &t, const Config &c){ double z_effective(const Projectile &p,const Target &t, const Config &c){

17
catima.cpp
View File

@ -336,15 +336,16 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
double res; double res;
//calculate 1st point to have i-1 element ready for loop //calculate 1st point to have i-1 element ready for loop
res = integrator.integrate(fdedx,Ezero,energy_table(0)); //res = integrator.integrate(fdedx,Ezero,energy_table(0));
res = p.A*res; //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] = res;
res = integrator.integrate(fomega,Ezero,energy_table(0)); //res = da2dx(p,energy_table(0),t)*res;
res = p.A*res; dp.angular_variance[0] = 0.0;
dp.range_straggling[0]=res;
//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++){ for(int i=1;i<max_datapoints;i++){
res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i)); res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));

2
constants.h
View File

@ -3,7 +3,7 @@
namespace catima { namespace catima {
constexpr double Ezero = 9E-4; // lowest E to calculate, below taken as 0 constexpr double Ezero = 1E-3; // lowest E to calculate, below taken as 0
constexpr double logEmin = -3; // log of minimum energy constexpr double logEmin = -3; // log of minimum energy
constexpr double logEmax = 5.0; // log of max energy constexpr double logEmax = 5.0; // log of max energy
constexpr int max_datapoints = 500; // how many datapoints between logEmin and logEmax constexpr int max_datapoints = 500; // how many datapoints between logEmin and logEmax