diff --git a/calculations.cpp b/calculations.cpp index e91653b..955ac4d 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -197,7 +197,7 @@ double bethek_lindhard(const Projectile &p){ std::complex cexir_den (sk,-eta); std::complex cexir = std::sqrt(cexir_n/cexir_den); std::complex csketa (sk + 1.0, eta); - std::complex cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag()); + std::complex cpiske(0.0,(PI*(l-sk)/2.0) - lngamma(csketa).imag()); std::complex cedr = cexir*std::exp(cpiske); double H=0; @@ -206,7 +206,7 @@ double bethek_lindhard(const Projectile &p){ std::complex cmsketa (-sk + 1.0, eta); std::complex cexis_den (-sk,-eta); std::complex cexis = std::sqrt(cexir_n/cexis_den); - std::complex cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); + std::complex cpimske(0.0,(PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); std::complex ceds = cexis*std::exp(cpimske); std::complex cmbeta_gamma_R(0,-beta_gamma_R); std::complex c2beta_gamma_R(0,2.0*beta_gamma_R); @@ -315,7 +315,7 @@ double bethek_lindhard_X(const Projectile &p){ std::complex cexir_den (sk,-eta); std::complex cexir = std::sqrt(cexir_n/cexir_den); std::complex csketa (sk + 1.0, eta); - std::complex cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag()); + std::complex cpiske(0.0,(PI*(l-sk)/2.0) - lngamma(csketa).imag()); std::complex cedr = cexir*std::exp(cpiske); double H=0; @@ -324,7 +324,7 @@ double bethek_lindhard_X(const Projectile &p){ std::complex cmsketa (-sk + 1.0, eta); std::complex cexis_den (-sk,-eta); std::complex cexis = std::sqrt(cexir_n/cexis_den); - std::complex cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); + std::complex cpimske(0.0,(PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); std::complex ceds = cexis*std::exp(cpimske); std::complex cmbeta_gamma_R(0,-beta_gamma_R); std::complex c2beta_gamma_R(0,2.0*beta_gamma_R); @@ -425,11 +425,11 @@ double pair_production(const Projectile &p, const Target &t){ double l0 = log(2.0*gamma); double Zt13 = 183.0*power(t.Z,-1.0/3.0); double L0screen = (19.0/9.0)*log(Zt13/(1.0 + 25.018803808*Zt13/gamma)); - double L1 = (4178.0/(81.0*M_PI*M_PI)) + double L1 = (4178.0/(81.0*PI*PI)) - (21.0/27.0) - - (248.0*l0/(27.0*M_PI*M_PI)) - + ( ((28.0*l0/(9.0)) - 446.0/27.0)*logd2/(M_PI*M_PI)) - + (14.0*logd2*logd2/(9.0*M_PI*M_PI)); + - (248.0*l0/(27.0*PI*PI)) + + ( ((28.0*l0/(9.0)) - 446.0/27.0)*logd2/(PI*PI)) + + (14.0*logd2*logd2/(9.0*PI*PI)); L1 *= d; double s = 0.25*dedx_constant*fine_structure*fine_structure*p.Z*p.Z*t.Z*t.Z*gamma*(1.0+(1.0/t.Z))*(L0screen + L1)/t.A; return (s<0.0)?0.0:s; @@ -440,7 +440,7 @@ double bremsstrahlung(const Projectile &p, const Target &t){ double R = 1.18*(catima::power(p.A, 1.0/3.0) + catima::power(t.A, 1.0/3.0)); double Lbs = log1p(2.0*gamma*0.1*hbar*c_light/atomic_mass_unit/R/p.A); double C = dedx_constant*fine_structure*(electron_mass/atomic_mass_unit); - 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); + 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*PI); }; double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){ @@ -817,7 +817,7 @@ std::complex lngamma( const std::complex &z ) double aterm1=y*log(r); double aterm2=(x+0.5)*atan2(y,(x+5.5))-y; double lterm1=(x+0.5)*log(r); - double lterm2=-y*atan2(y,(x+5.5)) - (x+5.5) + 0.5*log(2.0*M_PI); + double lterm2=-y*atan2(y,(x+5.5)) - (x+5.5) + 0.5*log(2.0*PI); double num=0.0; double denom=1.000000000190015; for(int j=1;j<7;j++){ @@ -831,8 +831,8 @@ std::complex lngamma( const std::complex &z ) double lterm3 = 0.5*log(num*num + denom*denom); std::complex result(lterm1+lterm2+lterm3,aterm1+aterm2+aterm3); if(z.real() < 0){ - std::complex lpi(log(M_PI), 0.0); - result = lpi - (result + std::log(std::sin(M_PI*z))); + std::complex lpi(log(PI), 0.0); + result = lpi - (result + std::log(std::sin(PI*z))); } return(result); } diff --git a/constants.h b/constants.h index 5d18b97..391d27a 100644 --- a/constants.h +++ b/constants.h @@ -24,6 +24,7 @@ constexpr bool reactions = false; // constants constexpr double PI = 3.1415926535897932384626433832795; +constexpr double LN10 = 2.30258509299404568401799145468; constexpr double Avogadro = 6.022140857; // 10^23 constexpr double electron_mass = 0.510998928; // MeV/c^2 constexpr double atomic_mass_unit = 931.4940954; // MeV/c^2 diff --git a/storage.h b/storage.h index 60fd651..3005536 100644 --- a/storage.h +++ b/storage.h @@ -40,7 +40,7 @@ namespace catima{ EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){ step = (logmax-logmin)/(N - 1.0); for(auto i=0;i=values[N-1]-numeric_epsilon)return N-1; int i = static_cast (std::floor(lxval/step)); @@ -66,7 +66,7 @@ namespace catima{ template int EnergyTable_index(const EnergyTable &table, double val){ if(valtable.values[table.num-1])return -1; - double lxval = (log(val/table.values[0])/M_LN10); + double lxval = (log(val/table.values[0])/LN10); int i = static_cast( std::floor(lxval/table.step)); if(val >= table.values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision return i; diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 15881df..6d461ba 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -478,6 +478,6 @@ using namespace std; TEST_CASE("constants"){ using namespace catima; CHECK(0.1*hbar*c_light/atomic_mass_unit == approx(0.21183,0.0001)); - CHECK(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*M_PI) == approx(5.21721169334564e-7).R(1e-3)); + CHECK(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*PI) == approx(5.21721169334564e-7).R(1e-3)); } diff --git a/tests/test_dedx_range.cpp b/tests/test_dedx_range.cpp index cc13cf0..30974e7 100644 --- a/tests/test_dedx_range.cpp +++ b/tests/test_dedx_range.cpp @@ -8,6 +8,7 @@ using std::cout; using std::endl; +using catima::LN10; double logEmax = catima::logEmax-0.01; double logEmin = catima::logEmin+0.001; @@ -16,7 +17,7 @@ template struct EnergyTable{ constexpr EnergyTable():values(),num(N){ for(auto i=0;i #include "testutils.h" -using namespace std; - #include "catima/catima.h" #include "catima/storage.h" +using namespace std; +using catima::LN10; TEST_CASE("datapoints equal operator"){ catima::Projectile p{12,6,6,1000}; @@ -116,13 +116,13 @@ using namespace std; TEST_CASE("energy table"){ double step = (catima::logEmax - catima::logEmin)/(catima::max_datapoints-1); CHECK(catima::energy_table.step==step); - CHECK(catima::energy_table.values[0]==exp(M_LN10*(catima::logEmin))); - CHECK(catima::energy_table.values[1]==exp(M_LN10*(catima::logEmin+step))); - CHECK(catima::energy_table.values[2]==exp(M_LN10*(catima::logEmin+2.0*step))); - CHECK(catima::energy_table.values[3]==exp(M_LN10*(catima::logEmin+3.0*step))); - CHECK(catima::energy_table.values[4]==exp(M_LN10*(catima::logEmin+4.0*step))); - CHECK(catima::energy_table.values[5]==exp(M_LN10*(catima::logEmin+5.0*step))); - CHECK(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(M_LN10*(catima::logEmax))).epsilon(1e-6)); + CHECK(catima::energy_table.values[0]==approx(exp(LN10*(catima::logEmin))).R(1e-9)); + CHECK(catima::energy_table.values[1]==approx(exp(LN10*(catima::logEmin+step))).R(1e-9)); + CHECK(catima::energy_table.values[2]==approx(exp(LN10*(catima::logEmin+2.0*step))).R(1e-9)); + CHECK(catima::energy_table.values[3]==approx(exp(LN10*(catima::logEmin+3.0*step))).R(1e-9)); + CHECK(catima::energy_table.values[4]==approx(exp(LN10*(catima::logEmin+4.0*step))).R(1e-9)); + CHECK(catima::energy_table.values[5]==approx(exp(LN10*(catima::logEmin+5.0*step))).R(1e-9)); + CHECK(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(LN10*(catima::logEmax))).epsilon(1e-6)); } TEST_CASE("indexing"){ double val, dif; @@ -137,4 +137,4 @@ using namespace std; CHECK(catima::energy_table.index(val)==i); CHECK(catima::energy_table.index(val+0.5*dif)==i); } - } \ No newline at end of file + }