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

windows build fix

Browse Source
This commit is contained in:
U-MEDAUSTRON\APR 2020-12-02 00:10:40 +01:00
parent be0e763db4
commit 8e8921c9ba
6 changed files with 29 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
calculations.cpp
View File

@ -197,7 +197,7 @@ double bethek_lindhard(const Projectile &p){
std::complex<double> cexir_den (sk,-eta); std::complex<double> cexir_den (sk,-eta);
std::complex<double> cexir = std::sqrt(cexir_n/cexir_den); std::complex<double> cexir = std::sqrt(cexir_n/cexir_den);
std::complex<double> csketa (sk + 1.0, eta); std::complex<double> csketa (sk + 1.0, eta);
std::complex<double> cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag()); std::complex<double> cpiske(0.0,(PI*(l-sk)/2.0) - lngamma(csketa).imag());
std::complex<double> cedr = cexir*std::exp(cpiske); std::complex<double> cedr = cexir*std::exp(cpiske);
double H=0; double H=0;
@ -206,7 +206,7 @@ double bethek_lindhard(const Projectile &p){
std::complex<double> cmsketa (-sk + 1.0, eta); std::complex<double> cmsketa (-sk + 1.0, eta);
std::complex<double> cexis_den (-sk,-eta); std::complex<double> cexis_den (-sk,-eta);
std::complex<double> cexis = std::sqrt(cexir_n/cexis_den); std::complex<double> cexis = std::sqrt(cexir_n/cexis_den);
std::complex<double> cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); std::complex<double> cpimske(0.0,(PI*(l+sk)/2.0) - lngamma(cmsketa).imag());
std::complex<double> ceds = cexis*std::exp(cpimske); std::complex<double> ceds = cexis*std::exp(cpimske);
std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R); std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R);
std::complex<double> c2beta_gamma_R(0,2.0*beta_gamma_R); std::complex<double> c2beta_gamma_R(0,2.0*beta_gamma_R);
@ -315,7 +315,7 @@ double bethek_lindhard_X(const Projectile &p){
std::complex<double> cexir_den (sk,-eta); std::complex<double> cexir_den (sk,-eta);
std::complex<double> cexir = std::sqrt(cexir_n/cexir_den); std::complex<double> cexir = std::sqrt(cexir_n/cexir_den);
std::complex<double> csketa (sk + 1.0, eta); std::complex<double> csketa (sk + 1.0, eta);
std::complex<double> cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag()); std::complex<double> cpiske(0.0,(PI*(l-sk)/2.0) - lngamma(csketa).imag());
std::complex<double> cedr = cexir*std::exp(cpiske); std::complex<double> cedr = cexir*std::exp(cpiske);
double H=0; double H=0;
@ -324,7 +324,7 @@ double bethek_lindhard_X(const Projectile &p){
std::complex<double> cmsketa (-sk + 1.0, eta); std::complex<double> cmsketa (-sk + 1.0, eta);
std::complex<double> cexis_den (-sk,-eta); std::complex<double> cexis_den (-sk,-eta);
std::complex<double> cexis = std::sqrt(cexir_n/cexis_den); std::complex<double> cexis = std::sqrt(cexir_n/cexis_den);
std::complex<double> cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag()); std::complex<double> cpimske(0.0,(PI*(l+sk)/2.0) - lngamma(cmsketa).imag());
std::complex<double> ceds = cexis*std::exp(cpimske); std::complex<double> ceds = cexis*std::exp(cpimske);
std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R); std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R);
std::complex<double> c2beta_gamma_R(0,2.0*beta_gamma_R); std::complex<double> 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 l0 = log(2.0*gamma);
double Zt13 = 183.0*power(t.Z,-1.0/3.0); 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 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) - (21.0/27.0)
- (248.0*l0/(27.0*M_PI*M_PI)) - (248.0*l0/(27.0*PI*PI))
+ ( ((28.0*l0/(9.0)) - 446.0/27.0)*logd2/(M_PI*M_PI)) + ( ((28.0*l0/(9.0)) - 446.0/27.0)*logd2/(PI*PI))
+ (14.0*logd2*logd2/(9.0*M_PI*M_PI)); + (14.0*logd2*logd2/(9.0*PI*PI));
L1 *= d; 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; 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; 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 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 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); 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){ double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){
@ -817,7 +817,7 @@ std::complex<double> lngamma( const std::complex<double> &z )
double aterm1=y*log(r); double aterm1=y*log(r);
double aterm2=(x+0.5)*atan2(y,(x+5.5))-y; double aterm2=(x+0.5)*atan2(y,(x+5.5))-y;
double lterm1=(x+0.5)*log(r); 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 num=0.0;
double denom=1.000000000190015; double denom=1.000000000190015;
for(int j=1;j<7;j++){ for(int j=1;j<7;j++){
@ -831,8 +831,8 @@ std::complex<double> lngamma( const std::complex<double> &z )
double lterm3 = 0.5*log(num*num + denom*denom); double lterm3 = 0.5*log(num*num + denom*denom);
std::complex<double> result(lterm1+lterm2+lterm3,aterm1+aterm2+aterm3); std::complex<double> result(lterm1+lterm2+lterm3,aterm1+aterm2+aterm3);
if(z.real() < 0){ if(z.real() < 0){
std::complex<double> lpi(log(M_PI), 0.0); std::complex<double> lpi(log(PI), 0.0);
result = lpi - (result + std::log(std::sin(M_PI*z))); result = lpi - (result + std::log(std::sin(PI*z)));
} }
return(result); return(result);
} }

1
constants.h
View File

@ -24,6 +24,7 @@ constexpr bool reactions = false;
// constants // constants
constexpr double PI = 3.1415926535897932384626433832795; constexpr double PI = 3.1415926535897932384626433832795;
constexpr double LN10 = 2.30258509299404568401799145468;
constexpr double Avogadro = 6.022140857; // 10^23 constexpr double Avogadro = 6.022140857; // 10^23
constexpr double electron_mass = 0.510998928; // MeV/c^2 constexpr double electron_mass = 0.510998928; // MeV/c^2
constexpr double atomic_mass_unit = 931.4940954; // MeV/c^2 constexpr double atomic_mass_unit = 931.4940954; // MeV/c^2

6
storage.h
View File

@ -40,7 +40,7 @@ namespace catima{
EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){ EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){
step = (logmax-logmin)/(N - 1.0); step = (logmax-logmin)/(N - 1.0);
for(auto i=0;i<N;i++){ for(auto i=0;i<N;i++){
values[i]=exp(M_LN10*(logmin + ((double)i)*step)); values[i]=exp(LN10*(logmin + ((double)i)*step));
} }
} }
double operator()(int i)const{return values[i];} double operator()(int i)const{return values[i];}
@ -51,7 +51,7 @@ namespace catima{
double* begin(){return values;} double* begin(){return values;}
double* end(){return &values[num];} double* end(){return &values[num];}
int index(double v)const noexcept{ int index(double v)const noexcept{
double lxval = (log(v/values[0])/M_LN10); double lxval = (log(v/values[0])/LN10);
if(v<values[0] || step==0.0)return -1; if(v<values[0] || step==0.0)return -1;
if(v>=values[N-1]-numeric_epsilon)return N-1; if(v>=values[N-1]-numeric_epsilon)return N-1;
int i = static_cast<int> (std::floor(lxval/step)); int i = static_cast<int> (std::floor(lxval/step));
@ -66,7 +66,7 @@ namespace catima{
template<int N> template<int N>
int EnergyTable_index(const EnergyTable<N> &table, double val){ int EnergyTable_index(const EnergyTable<N> &table, double val){
if(val<table.values[0] || val>table.values[table.num-1])return -1; if(val<table.values[0] || val>table.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<int>( std::floor(lxval/table.step)); int i = static_cast<int>( std::floor(lxval/table.step));
if(val >= table.values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision if(val >= table.values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision
return i; return i;

2
tests/test_calculations.cpp
View File

@ -478,6 +478,6 @@ using namespace std;
TEST_CASE("constants"){ TEST_CASE("constants"){
using namespace catima; using namespace catima;
CHECK(0.1*hbar*c_light/atomic_mass_unit == approx(0.21183,0.0001)); 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));
} }

3
tests/test_dedx_range.cpp
View File

@ -8,6 +8,7 @@
using std::cout; using std::cout;
using std::endl; using std::endl;
using catima::LN10;
double logEmax = catima::logEmax-0.01; double logEmax = catima::logEmax-0.01;
double logEmin = catima::logEmin+0.001; double logEmin = catima::logEmin+0.001;
@ -16,7 +17,7 @@ template<int N>
struct EnergyTable{ struct EnergyTable{
constexpr EnergyTable():values(),num(N){ constexpr EnergyTable():values(),num(N){
for(auto i=0;i<N;i++){ for(auto i=0;i<N;i++){
values[i]=exp(M_LN10*(logEmin + ((double)i)*(logEmax-logEmin)/(N - 1.0))); values[i]=exp(LN10*(logEmin + ((double)i)*(logEmax-logEmin)/(N - 1.0)));
} }
} }
double operator()(int i)const{return values[i];} double operator()(int i)const{return values[i];}

18
tests/test_storage.cpp
View File

@ -3,10 +3,10 @@
#include "doctest.h" #include "doctest.h"
#include <math.h> #include <math.h>
#include "testutils.h" #include "testutils.h"
using namespace std;
#include "catima/catima.h" #include "catima/catima.h"
#include "catima/storage.h" #include "catima/storage.h"
using namespace std;
using catima::LN10;
TEST_CASE("datapoints equal operator"){ TEST_CASE("datapoints equal operator"){
catima::Projectile p{12,6,6,1000}; catima::Projectile p{12,6,6,1000};
@ -116,13 +116,13 @@ using namespace std;
TEST_CASE("energy table"){ TEST_CASE("energy table"){
double step = (catima::logEmax - catima::logEmin)/(catima::max_datapoints-1); double step = (catima::logEmax - catima::logEmin)/(catima::max_datapoints-1);
CHECK(catima::energy_table.step==step); CHECK(catima::energy_table.step==step);
CHECK(catima::energy_table.values[0]==exp(M_LN10*(catima::logEmin))); CHECK(catima::energy_table.values[0]==approx(exp(LN10*(catima::logEmin))).R(1e-9));
CHECK(catima::energy_table.values[1]==exp(M_LN10*(catima::logEmin+step))); CHECK(catima::energy_table.values[1]==approx(exp(LN10*(catima::logEmin+step))).R(1e-9));
CHECK(catima::energy_table.values[2]==exp(M_LN10*(catima::logEmin+2.0*step))); CHECK(catima::energy_table.values[2]==approx(exp(LN10*(catima::logEmin+2.0*step))).R(1e-9));
CHECK(catima::energy_table.values[3]==exp(M_LN10*(catima::logEmin+3.0*step))); CHECK(catima::energy_table.values[3]==approx(exp(LN10*(catima::logEmin+3.0*step))).R(1e-9));
CHECK(catima::energy_table.values[4]==exp(M_LN10*(catima::logEmin+4.0*step))); CHECK(catima::energy_table.values[4]==approx(exp(LN10*(catima::logEmin+4.0*step))).R(1e-9));
CHECK(catima::energy_table.values[5]==exp(M_LN10*(catima::logEmin+5.0*step))); 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(M_LN10*(catima::logEmax))).epsilon(1e-6)); CHECK(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(LN10*(catima::logEmax))).epsilon(1e-6));
} }
TEST_CASE("indexing"){ TEST_CASE("indexing"){
double val, dif; double val, dif;