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

energy table indexing

Browse Source
This commit is contained in:
hrocho 2022-04-22 22:52:10 +02:00
parent 288349405e
commit 71d97e81c4
15 changed files with 224 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
CMakeLists.txt
View File

@ -13,15 +13,16 @@ option(STORE_SPLINES "store splines, if disables splines are always recreated" O
option(GSL_INTEGRATION "use GSL integration" OFF) option(GSL_INTEGRATION "use GSL integration" OFF)
option(GSL_INTERPOLATION "use GSL inteRPOLATION" OFF) option(GSL_INTERPOLATION "use GSL inteRPOLATION" OFF)
option(THIN_TARGET_APPROXIMATION "thin target approximation" ON) option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
option(ET_CALCULATED_INDEX "calculate energy table index, otherwise search" ON)
option(GENERATE_DATA "make data tables generator" OFF) option(GENERATE_DATA "make data tables generator" OFF)
option(PYTHON_WHEEL "make python wheel" OFF) option(PYTHON_WHEEL "make python wheel" OFF)
######## build type ############ ######## build type ############
if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug") if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug")
set(CMAKE_BUILD_TYPE "Release") set(CMAKE_BUILD_TYPE "Release")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") #set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -ffast-math")
MESSAGE(STATUS "Build type Release") MESSAGE(STATUS "Build type Release")
else() else()
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g") #set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")
if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang" OR ${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang" OR ${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare")
endif() endif()
@ -56,18 +57,21 @@ endif()
#endif(nurex_FOUND) #endif(nurex_FOUND)
find_package(fmt QUIET) find_package(fmt QUIET)
if(NOT fmt_FOUND) function(check_fmt)
message("fmt library not found, trying to dowload") if(NOT fmt_FOUND)
include(FetchContent) message("fmt library not found, trying to dowload")
FetchContent_Declare( include(FetchContent)
fmt FetchContent_Declare(
GIT_REPOSITORY https://github.com/fmtlib/fmt.git fmt
GIT_TAG 8.1.1 GIT_REPOSITORY https://github.com/fmtlib/fmt.git
) GIT_TAG 8.1.1
set(CMAKE_POSITION_INDEPENDENT_CODE ON) )
set(FMT_INSTALL ON) set(CMAKE_POSITION_INDEPENDENT_CODE ON)
FetchContent_MakeAvailable(fmt) set(FMT_INSTALL ON)
endif(NOT fmt_FOUND) FetchContent_MakeAvailable(fmt)
endif(NOT fmt_FOUND)
endfunction(check_fmt)
configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in" configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
"${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h") "${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h")
@ -85,8 +89,9 @@ file(GLOB HEADERS *.h libs/*.h)
add_library(catima ${SOURCES}) add_library(catima ${SOURCES})
set_target_properties(catima PROPERTIES set_target_properties(catima PROPERTIES
POSITION_INDEPENDENT_CODE ON POSITION_INDEPENDENT_CODE ON
LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
) )
target_link_libraries(catima PUBLIC ${EXTRA_LIBS}) target_link_libraries(catima PUBLIC ${EXTRA_LIBS})
@ -124,6 +129,7 @@ if(PYTHON_MODULE)
FetchContent_MakeAvailable(pybind11) FetchContent_MakeAvailable(pybind11)
endif(NOT pybind11_FOUND) endif(NOT pybind11_FOUND)
check_fmt()
#set(PYBIND11_CPP_STANDARD -std=c++14) #set(PYBIND11_CPP_STANDARD -std=c++14)
pybind11_add_module(pycatima pymodule/pycatima.cpp) pybind11_add_module(pycatima pymodule/pycatima.cpp)
target_include_directories(pycatima PUBLIC target_include_directories(pycatima PUBLIC
@ -135,6 +141,7 @@ endif(PYTHON_MODULE )
configure_file("${PROJECT_SOURCE_DIR}/pymodule/setup.py.in" "${PROJECT_BINARY_DIR}/setup.py") configure_file("${PROJECT_SOURCE_DIR}/pymodule/setup.py.in" "${PROJECT_BINARY_DIR}/setup.py")
if(PYTHON_WHEEL) if(PYTHON_WHEEL)
check_fmt()
execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py bdist_wheel) execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py bdist_wheel)
endif(PYTHON_WHEEL) endif(PYTHON_WHEEL)

1
build_config.in
View File

@ -8,5 +8,6 @@
#cmakedefine GLOBAL #cmakedefine GLOBAL
#cmakedefine REACTIONS #cmakedefine REACTIONS
#cmakedefine NUREX #cmakedefine NUREX
#cmakedefine ET_CALCULATED_INDEX
#endif #endif

43
calculations.cpp
View File

@ -1,4 +1,4 @@
#include <math.h> #include <cmath>
#include <algorithm> #include <algorithm>
#include <cassert> #include <cassert>
#include "catima/calculations.h" #include "catima/calculations.h"
@ -21,7 +21,7 @@ extern "C"
namespace catima{ namespace catima{
double reduced_energy_loss_unit(const Projectile &p, const Target &t){ double reduced_energy_loss_unit(const Projectile &p, const Target &t){
double zpowers = pow(p.Z,0.23)+pow(t.Z,0.23); double zpowers = std::pow(p.Z,0.23)+std::pow(t.Z,0.23);
double asum = p.A + t.A; double asum = p.A + t.A;
return 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV return 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
} }
@ -473,7 +473,7 @@ double p_from_T(double T, double M){
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);
double factor=4.8184E-3*pow(z1+z2,8.0/3.0)/m2; double factor=4.8184E-3*std::pow(z1+z2,8.0/3.0)/m2;
return factor*beta2/fine_structure/fine_structure; return factor*beta2/fine_structure/fine_structure;
} }
@ -485,7 +485,7 @@ double angular_scattering_power(const Projectile &p, const Target &t, double Es2
double lr = radiation_length(t.Z,t.A); double lr = radiation_length(t.Z,t.A);
//constexpr double Es2 = 198.81; //constexpr double Es2 = 198.81;
//constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass; //constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass;
return Es2 * pow(p.Z,2)/(lr*pow(_p*beta,2)); return Es2 * ipow(p.Z,2)/(lr*ipow(_p*beta,2));
} }
double angular_scattering_power(const Projectile &p, const Material &mat, double Es2){ double angular_scattering_power(const Projectile &p, const Material &mat, double Es2){
@ -496,7 +496,7 @@ double angular_scattering_power(const Projectile &p, const Material &mat, double
double X0 = radiation_length(mat); double X0 = radiation_length(mat);
//constexpr double Es2 = 198.81; //constexpr double Es2 = 198.81;
//constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass; //constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass;
return Es2 * pow(p.Z,2)/(X0*pow(_p*beta,2)); return Es2 * ipow(p.Z,2)/(X0*ipow(_p*beta,2));
} }
double angular_scattering_power_xs(const Projectile &p, const Material &mat, double p1, double beta1, double Es2){ double angular_scattering_power_xs(const Projectile &p, const Material &mat, double p1, double beta1, double Es2){
@ -507,10 +507,10 @@ double angular_scattering_power_xs(const Projectile &p, const Material &mat, dou
double pv = _p*beta; double pv = _p*beta;
double p1v1 = p1*beta1; double p1v1 = p1*beta1;
double Xs = scattering_length(mat); double Xs = scattering_length(mat);
double cl1 = log10(1-std::pow(pv/p1v1,2)); double cl1 = log10(1-ipow(pv/p1v1,2));
double cl2 = log10(pv); double cl2 = log10(pv);
double f = 0.5244 + 0.1975*cl1 + 0.2320*cl2 - (0.0098*cl2*cl1); double f = 0.5244 + 0.1975*cl1 + 0.2320*cl2 - (0.0098*cl2*cl1);
return f*Es2 * pow(p.Z,2)/(Xs*pow(_p*beta,2)); return f*Es2 * ipow(p.Z,2)/(Xs*ipow(_p*beta,2));
} }
/// radiation lengths are taken from Particle Data Group 2014 /// radiation lengths are taken from Particle Data Group 2014
@ -542,7 +542,7 @@ double radiation_length(int z, double m){
if(z==92){return 6.00;} if(z==92){return 6.00;}
double z2 = z*z; double z2 = z*z;
double z_13 = 1.0/pow(z,1./3.); double z_13 = 1.0/std::pow(z,1./3.);
double z_23 = z_13*z_13; double z_23 = z_13*z_13;
double a2 = fine_structure*fine_structure*z2; double a2 = fine_structure*fine_structure*z2;
double a4 = a2*a2; double a4 = a2*a2;
@ -619,11 +619,11 @@ double dedx_variance(const Projectile &p, const Target &t, const Config &c){
double beta = beta_from_T(p.T); double beta = beta_from_T(p.T);
double beta2 = beta*beta; double beta2 = beta*beta;
double zp_eff = z_effective(p,t,c); double zp_eff = z_effective(p,t,c);
double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A; double f = domega2dx_constant*ipow(zp_eff,2)*t.Z/t.A;
if( (c.calculation == omega_types::atima) ){ if( (c.calculation == omega_types::atima) ){
cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)* cor = 24.89 * std::pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364))); log( 2.0*electron_mass*1e6*beta2/(33.05*std::pow(t.Z,1.6364)));
cor = std::max(cor, 0.0 ); cor = std::max(cor, 0.0 );
} }
//double X = bethek_lindhard_X(p); //double X = bethek_lindhard_X(p);
@ -671,13 +671,13 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){
} }
double z_eff_Pierce_Blann(double z, double beta){ double z_eff_Pierce_Blann(double z, double beta){
return z*(1.0-exp(-0.95*fine_structure_inverted*beta/pow(z,2.0/3.0))); return z*(1.0-exp(-0.95*fine_structure_inverted*beta/std::pow(z,2.0/3.0)));
} }
double z_eff_Anthony_Landford(double pz, double beta, double tz){ double z_eff_Anthony_Landford(double pz, double beta, double tz){
double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz); double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz);
double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz); double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz);
return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0)))); return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/std::pow(pz,2.0/3.0))));
} }
double z_eff_Hubert(double pz, double E, double tz){ double z_eff_Hubert(double pz, double E, double tz){
@ -705,7 +705,7 @@ double z_eff_Hubert(double pz, double E, double tz){
x4 = 0.5218 + 0.02521*lntz; x4 = 0.5218 + 0.02521*lntz;
} }
return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4))); return pz*(1-x1*exp(-x2*std::pow(E,x3)*std::pow(pz,-x4)));
} }
double z_eff_Winger(double pz, double beta, double tz){ double z_eff_Winger(double pz, double beta, double tz){
@ -733,7 +733,7 @@ double z_eff_Winger(double pz, double beta, double tz){
if(tz==2){ if(tz==2){
tz = 2.8; tz = 2.8;
} }
x = beta /0.012 /pow(pz,0.45); x = beta /0.012 /std::pow(pz,0.45);
lnz =log(pz); lnz =log(pz);
lnzt=log(tz); lnzt=log(tz);
@ -752,13 +752,14 @@ double z_eff_global(double pz, double E, double tz){
else if(E<30.0 || pz<29){ else if(E<30.0 || pz<29){
return z_eff_Pierce_Blann(pz, E); return z_eff_Pierce_Blann(pz, E);
} }
else else{
#ifdef GLOBAL #ifdef GLOBAL
return global_qmean(pz, tz, E); return global_qmean(pz, tz, E);
#else #else
assert(false); assert(false);
return -1; return -1;
#endif #endif
}
} }
double z_eff_Schiwietz(double pz, double beta, double tz){ double z_eff_Schiwietz(double pz, double beta, double tz){
@ -766,10 +767,10 @@ double z_eff_Schiwietz(double pz, double beta, double tz){
double c1, c2; double c1, c2;
double x; double x;
scaled_velocity = catima::power(pz,-0.543)*beta/bohr_velocity; scaled_velocity = std::pow(pz,-0.543)*beta/bohr_velocity;
c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9); c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9);
c2 = 1+0.030*scaled_velocity*log(tz); c2 = 1+0.030*scaled_velocity*log(tz);
x = c1*catima::power(scaled_velocity/c2/1.54,1+(1.83/pz)); x = c1*std::pow(scaled_velocity/c2/1.54,1+(1.83/pz));
return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) ); return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) );
} }
@ -789,17 +790,17 @@ double z_eff_atima14(double pz, double T, double tz){
double c2 = 0.28; double c2 = 0.28;
double c3 = 0.04; double c3 = 0.04;
qglobal = z_eff_global(pz,T,tz); qglobal = z_eff_global(pz,T,tz);
qglobal = (qglobal - qpb)*c1/catima::power(tz+1,c2)*(1-exp(-c3*T)) + qpb; qglobal = (qglobal - qpb)*c1/std::pow(tz+1,c2)*(1-exp(-c3*T)) + qpb;
} }
emax = 1500.; emax = 1500.;
emin = 1000.; emin = 1000.;
if(T>emax){ if(T>emax){
qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1))); qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
qmean = qhigh; qmean = qhigh;
} }
else if(T>=emin && T<=emax){ else if(T>=emin && T<=emax){
qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1))); qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
if(pz<=28){ if(pz<=28){
qwinger = z_eff_Winger(pz,beta,tz); qwinger = z_eff_Winger(pz,beta,tz);
qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin); qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin);

8
calculations.h
View File

@ -230,5 +230,13 @@ namespace catima{
inline double power(double x, double y){ inline double power(double x, double y){
return exp(log(x)*y); return exp(log(x)*y);
} }
inline double ipow(double x, int y){
if(y==0)return 1.0;
else if(y==1) return x;
else if(y==2) return x*x;
else return std::pow(x,y);
}
} }
#endif #endif

10
catima.cpp
View File

@ -1,5 +1,5 @@
#include <iostream> #include <iostream>
#include <math.h> #include <cmath>
#include <algorithm> #include <algorithm>
#include "catima/catima.h" #include "catima/catima.h"
#include "catima/constants.h" #include "catima/constants.h"
@ -133,7 +133,7 @@ double Tfr(const Projectile &p, double X, double Es2){
if(p.T<=0)return 0.0; if(p.T<=0)return 0.0;
double _p = p_from_T(p.T,p.A); double _p = p_from_T(p.T,p.A);
double beta = _p/((p.T+atomic_mass_unit)*p.A); double beta = _p/((p.T+atomic_mass_unit)*p.A);
return Es2 /(X*pow(_p*beta,2)); return Es2 /(X*ipow(_p*beta,2));
} }
double angular_variance(Projectile p, const Material &t, const Config &c, int order){ double angular_variance(Projectile p, const Material &t, const Config &c, int order){
@ -154,7 +154,7 @@ double angular_variance(Projectile p, const Material &t, const Config &c, int or
auto fx0p = [&](double x)->double{ auto fx0p = [&](double x)->double{
double e =energy_out(T,x*t.density(),range_spline); double e =energy_out(T,x*t.density(),range_spline);
double d = std::pow((rrange-x),order); double d = ipow((rrange-x),order);
double ff = 1; double ff = 1;
if(c.scattering == scattering_types::dhighland){ if(c.scattering == scattering_types::dhighland){
double l = x*t.density()/X0; double l = x*t.density()/X0;
@ -167,7 +167,7 @@ double angular_variance(Projectile p, const Material &t, const Config &c, int or
auto fx0p_2 = [&](double x)->double{ auto fx0p_2 = [&](double x)->double{
double e =energy_out(T,x*t.density(),range_spline); double e =energy_out(T,x*t.density(),range_spline);
double d = std::pow((rrange-x),order); double d = ipow((rrange-x),order);
return d*angular_scattering_power_xs(p(e),t,p1,beta1); return d*angular_scattering_power_xs(p(e),t,p1,beta1);
}; };
@ -175,7 +175,7 @@ double angular_variance(Projectile p, const Material &t, const Config &c, int or
if(c.scattering == scattering_types::gottschalk){ if(c.scattering == scattering_types::gottschalk){
return integrator.integrate(fx0p_2,0, rrange)*t.density(); return integrator.integrate(fx0p_2,0, rrange)*t.density();
} }
return integrator.integrate(fx0p,0, rrange)*t.density()*pow(p.Z,2)*Es2/X0; return integrator.integrate(fx0p,0, rrange)*t.density()*ipow(p.Z,2)*Es2/X0;
} }
double angular_straggling(Projectile p, const Material &t, const Config &c){ double angular_straggling(Projectile p, const Material &t, const Config &c){

6
cwrapper.cpp
View File

@ -26,6 +26,7 @@ extern "C" {
CatimaResult catima_calculate(double pa, int pz, double T, double ta, double tz, double thickness, double density){ CatimaResult catima_calculate(double pa, int pz, double T, double ta, double tz, double thickness, double density){
catima::default_config.z_effective = catima_defaults.z_effective; catima::default_config.z_effective = catima_defaults.z_effective;
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, thickness, density); catima::Material mat = make_material(ta,tz, thickness, density);
catima::Result r = catima::calculate(p(T),mat); catima::Result r = catima::calculate(p(T),mat);
@ -46,6 +47,7 @@ extern "C" {
double catima_Eout(double pa, int pz, double T, double ta, double tz, double thickness, double density){ double catima_Eout(double pa, int pz, double T, double ta, double tz, double thickness, double density){
catima::default_config.z_effective = catima_defaults.z_effective; catima::default_config.z_effective = catima_defaults.z_effective;
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, thickness, density); catima::Material mat = make_material(ta,tz, thickness, density);
return energy_out(p(T), mat); return energy_out(p(T), mat);
@ -53,6 +55,7 @@ extern "C" {
double catima_range(double pa, int pz, double T, double ta, double tz){ double catima_range(double pa, int pz, double T, double ta, double tz){
catima::default_config.z_effective = catima_defaults.z_effective; catima::default_config.z_effective = catima_defaults.z_effective;
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, 0, -1); catima::Material mat = make_material(ta,tz, 0, -1);
return range(p, mat); return range(p, mat);
@ -60,12 +63,14 @@ extern "C" {
double catima_range_straggling(double pa, int pz, double T, double ta, double tz){ double catima_range_straggling(double pa, int pz, double T, double ta, double tz){
catima::default_config.z_effective = catima_defaults.z_effective; catima::default_config.z_effective = catima_defaults.z_effective;
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, 0, -1); catima::Material mat = make_material(ta,tz, 0, -1);
return range_straggling(p, T, mat); return range_straggling(p, T, mat);
} }
double catima_angular_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){ double catima_angular_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, 0, -1); catima::Material mat = make_material(ta,tz, 0, -1);
@ -73,6 +78,7 @@ extern "C" {
} }
double catima_energy_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){ double catima_energy_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){
catima::default_config.scattering = 255;
catima::Projectile p(pa,pz); catima::Projectile p(pa,pz);
catima::Material mat = make_material(ta,tz, 0, -1); catima::Material mat = make_material(ta,tz, 0, -1);
return catima::energy_straggling_from_E(p,Tin,Tout,mat); return catima::energy_straggling_from_E(p,Tin,Tout,mat);

9
pymodule/pycatima.cpp
View File

@ -49,9 +49,14 @@ py::list storage_info(){
py::list get_energy_table(){ py::list get_energy_table(){
py::list r; py::list r;
for(auto e : energy_table){ for (size_t i = 0; i < energy_table.size(); i++)
r.append(e); {
r.append(energy_table[i]);
} }
//for(auto e : energy_table){
//r.append(e);
//}
return r; return r;
} }

5
pymodule/setup.py
View File

@ -5,10 +5,7 @@ from setuptools import setup
DIR = Path(__file__).parents[0] DIR = Path(__file__).parents[0]
SRC = [str((DIR/'pycatima.cpp').resolve())]#+[str(fname.resolve()) for fname in DIR.glob('../*.cpp')] SRC = [str((DIR/'pycatima.cpp').resolve())]
#SRCG = [str(fname.resolve()) for fname in DIR.glob('../global/*.c')]
#SRC += SRCG
print (SRC)
example_module = Pybind11Extension( example_module = Pybind11Extension(
'pycatima', 'pycatima',
SRC, SRC,

22
pymodule/setup.py.in Normal file
View File

@ -0,0 +1,22 @@
from pybind11.setup_helpers import Pybind11Extension, build_ext
from setuptools import setup
SRC = [str(('${PROJECT_SOURCE_DIR}/pymodule/pycatima.cpp').resolve())]
example_module = Pybind11Extension(
'pycatima',
SRC,
include_dirs=['${PROJECT_BINARY_DIR}/include','${PROJECT_SOURCE_DIR}/global'],
library_dirs=['${PROJECT_BINARY_DIR}/lib','${PROJECT_BINARY_DIR}','${PROJECT_BINARY_DIR}/Release'],
libraries=['catima']
)
setup(
name='pycatima',
version=1.71,
author='Andrej Prochazka',
author_email='hrocho@vodacionline.sk',
description='python interface to catima library',
url='https://github.com/hrosiak/catima',
ext_modules=[example_module],
cmdclass={"build_ext": build_ext},
)

10
spline.cpp
View File

@ -1,15 +1,7 @@
/*
* This is modification of Tino Kluge tk spline
* https://github.com/ttk592/spline/
*
* the modification is in LU caclulation,
* optimized for tridiagonal matrices
*/
#include "spline.h" #include "spline.h"
namespace catima{ namespace catima{
} // namespace tk }

8
spline.h
View File

@ -111,7 +111,6 @@ struct cspline_special{
cspline_special() = default; cspline_special() = default;
const T *table; const T *table;
const double *m_x;
const double *m_y; const double *m_y;
std::array<double,N> m_a,m_b,m_c; std::array<double,N> m_a,m_b,m_c;
double m_b0, m_c0; double m_b0, m_c0;
@ -120,7 +119,9 @@ struct cspline_special{
double operator()(double x)const{return evaluate(x);} double operator()(double x)const{return evaluate(x);}
double evaluate(double x) const double evaluate(double x) const
{ {
int idx=std::max( table->index(x), 0); const T& m_x = *table;
int idx=std::max( table->index(x), 0);
double h=x-m_x[idx]; double h=x-m_x[idx];
double interpol; double interpol;
if(x<m_x[0]) { if(x<m_x[0]) {
@ -139,6 +140,7 @@ struct cspline_special{
double deriv(double x) const double deriv(double x) const
{ {
const T& m_x = *table;
int idx=std::max( table->index(x), 0); int idx=std::max( table->index(x), 0);
double h=x-m_x[idx]; double h=x-m_x[idx];
@ -162,7 +164,7 @@ template<typename T>
cspline_special<T>::cspline_special(const T &x, cspline_special<T>::cspline_special(const T &x,
const std::vector<double>& y, const std::vector<double>& y,
bool boundary_second_deriv bool boundary_second_deriv
):table(&x),m_y(y.data()),m_x(x.values) ):table(&x),m_y(y.data())
{ {
static_assert (N>2, "N must be > 2"); static_assert (N>2, "N must be > 2");
tridiagonal_matrix<N> A{}; tridiagonal_matrix<N> A{};

4
storage.cpp
View File

@ -19,7 +19,11 @@
#include "catima/catima.h" #include "catima/catima.h"
namespace catima { namespace catima {
Data _storage; Data _storage;
#ifdef VETABLE
LogVArray<max_datapoints> energy_table(logEmin,logEmax);
#else
EnergyTable<max_datapoints> energy_table(logEmin,logEmax); EnergyTable<max_datapoints> energy_table(logEmin,logEmax);
#endif
bool operator==(const DataPoint &a, const DataPoint &b){ bool operator==(const DataPoint &a, const DataPoint &b){
if( (a.m == b.m) && (a.p == b.p) && (a.config == b.config)){ if( (a.m == b.m) && (a.p == b.p) && (a.config == b.config)){

108
storage.h
View File

@ -29,7 +29,7 @@
#include "catima/spline.h" #include "catima/spline.h"
//#define VETABLE
namespace catima{ namespace catima{
/** /**
@ -48,36 +48,31 @@ namespace catima{
static constexpr int size() {return N;}; static constexpr int size() {return N;};
double values[N]; double values[N];
double step; double step;
double* begin(){return values;} const double* begin()const{return values;}
double* end(){return &values[num];} const double* end()const{return &values[num];}
int index(double v)const noexcept{ int index(double v)const noexcept{
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;
#ifdef ET_CALCULATED_INDEX
double lxval = (std::log(v/values[0])/LN10);
int i = static_cast<int> (std::floor(lxval/step)); int i = static_cast<int> (std::floor(lxval/step));
if(v >= values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision if(v >= values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision
return i; return i;
#else
auto it=std::upper_bound(begin(),end(),v);
return int(it-begin())-1;
#endif
}; };
std::size_t num; std::size_t num;
}; };
extern EnergyTable<max_datapoints> energy_table; template<typename T>
double EnergyTable_interpolate(const T &table, double xval, double *y){
template<int N>
int EnergyTable_index(const EnergyTable<N> &table, double val){
if(val<table.values[0] || val>table.values[table.num-1])return -1;
double lxval = (log(val/table.values[0])/LN10);
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
return i;
}
template<int N>
double EnergyTable_interpolate(const EnergyTable<N> &table, double xval, double *y){
double r; double r;
if(xval<table.values[0] || xval>table.values[table.num-1])return 0.0; if(xval<table.values[0] || xval>table.values[table.size()-1])return 0.0;
if(xval==table.values[table.num-1])return y[table.num-1]; if(xval==table.values[table.size()-1])return y[table.size()-1];
int i = EnergyTable_index(table, xval); int i = table.index(xval);
double linstep = table.values[i+1] - table.values[i]; double linstep = table.values[i+1] - table.values[i];
if(linstep == 0.0)return table.values[i]; if(linstep == 0.0)return table.values[i];
double x = 1.0 - ((xval - table.values[i])/linstep); double x = 1.0 - ((xval - table.values[i])/linstep);
@ -85,6 +80,63 @@ namespace catima{
return r; return r;
} }
template <int N>
struct LogVArray{
LogVArray(double logmin, double logmax):logmin(logmin),logmax(logmax){
assert(logmax>logmin);
step = (logmax-logmin)/(N - 1.0);
}
double get_min()const noexcept{return logmin;}
double get_max()const noexcept{return logmax;}
constexpr static int size() noexcept{return N;}
constexpr double value(int i) const noexcept{return exp(LN10*(logmin + ((double)i)*step));}
double operator[](int i)const noexcept{return value(i);}
double operator()(int i)const noexcept{return value(i);}
double step_size()const noexcept{return step;}
int index(double v)const noexcept{
if(v<value(0) || step==0.0)return -1;
if(v>= (value(N-1)-numeric_epsilon))return N-1;
double lxval = (log(v/value(0))/LN10);
int i = static_cast<int> (std::floor(lxval/step));
if(v >= value(i+1)-numeric_epsilon)i++; // this is correction for floating point precision
return i;
}
double step=0.0;
private:
double logmin;
double logmax;
static_assert (N>2, "N must be more than 2");
};
template <int N>
struct LinearVArray{
LinearVArray(double min, double max):min(min),max(max){
if(max<=min)return;
step = (max-min)/(N-1);
}
double get_min()const noexcept{return min;}
double get_max()const noexcept{return max;}
constexpr static int size() noexcept{return N;}
double operator[](int i)const noexcept{return i*step + min;}
int index(double v)const noexcept{
if(v<min || step==0.0)return -1;
if(v>=max)return N-1;
assert(step>0.0);
return static_cast<int> (std::floor((v-min)/step));
}
private:
double step=0.0;
double min;
double max;
static_assert (N>2, "N must be more than 2");
};
#ifdef VETABLE
extern LogVArray<max_datapoints> energy_table;
#else
extern EnergyTable<max_datapoints> energy_table;
#endif
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
#ifdef GSL_INTERPOLATION #ifdef GSL_INTERPOLATION
/// Interpolation class, to store interpolated values /// Interpolation class, to store interpolated values
@ -107,12 +159,13 @@ namespace catima{
}; };
#endif #endif
template<typename xtype>
class InterpolatorCSpline{ class InterpolatorCSpline{
public: public:
using xtype = EnergyTable<max_datapoints>; //using xtype = EnergyTable<max_datapoints>;
InterpolatorCSpline()=default; InterpolatorCSpline()=default;
InterpolatorCSpline(const xtype &table, const std::vector<double> &y): InterpolatorCSpline(const xtype &table, const std::vector<double> &y):
min(table.values[0]), max(table.values[max_datapoints-1]), ss(table,y){} min(table[0]), max(table[max_datapoints-1]), ss(table,y){}
double operator()(double x)const{return eval(x);} double operator()(double x)const{return eval(x);}
double eval(double x)const{return ss.evaluate(x);} double eval(double x)const{return ss.evaluate(x);}
double derivative(double x)const{return ss.deriv(x);} double derivative(double x)const{return ss.deriv(x);}
@ -128,7 +181,14 @@ namespace catima{
#ifdef GSL_INTERPOLATION #ifdef GSL_INTERPOLATION
using Interpolator = InterpolatorGSL; using Interpolator = InterpolatorGSL;
#else #else
using Interpolator = InterpolatorCSpline; #ifdef VETABLE
//using Interpolator = InterpolatorSplineT<LogVArray<max_datapoints>>;
using Interpolator = InterpolatorCSpline<LogVArray<max_datapoints>>;
#else
//using Interpolator = InterpolatorSplineT<EnergyTable<max_datapoints>>;
using Interpolator = InterpolatorCSpline<EnergyTable<max_datapoints>>;
#endif
#endif #endif
#ifdef STORE_SPLINES #ifdef STORE_SPLINES

2
tests/test_calculations.cpp
View File

@ -425,7 +425,7 @@ using namespace std;
CHECK(res1.dEdxo == res2.dEdxo); CHECK(res1.dEdxo == res2.dEdxo);
CHECK(res1.tof == res2.tof); CHECK(res1.tof == res2.tof);
auto ra = catima::angular_straggling_from_E(p,res1.Ein,res1.Eout,graphite); auto ra = catima::angular_straggling_from_E(p(res1.Ein),res1.Eout,graphite);
CHECK(res1.sigma_a == approx(ra).R(1e-3)); CHECK(res1.sigma_a == approx(ra).R(1e-3));
auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite); auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite);

50
tests/test_storage.cpp
View File

@ -114,27 +114,45 @@ using catima::LN10;
} }
TEST_CASE("energy table"){ TEST_CASE("energy table"){
catima::LogVArray<catima::max_datapoints> etable(catima::logEmin,catima::logEmax);
catima::EnergyTable<catima::max_datapoints> energy_table(catima::logEmin,catima::logEmax);
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(energy_table.step==step);
CHECK(catima::energy_table.values[0]==approx(exp(LN10*(catima::logEmin))).R(1e-9)); CHECK(energy_table[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(energy_table[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(energy_table[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(energy_table[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(energy_table[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(energy_table[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)); CHECK(energy_table[catima::max_datapoints-1]==approx(exp(LN10*(catima::logEmax))).epsilon(1e-6));
CHECK(etable.step_size()==step);
CHECK(etable[0]==approx(exp(LN10*(catima::logEmin))).R(1e-9));
CHECK(etable[1]==approx(exp(LN10*(catima::logEmin+step))).R(1e-9));
CHECK(etable[2]==approx(exp(LN10*(catima::logEmin+2.0*step))).R(1e-9));
CHECK(etable[3]==approx(exp(LN10*(catima::logEmin+3.0*step))).R(1e-9));
CHECK(etable[4]==approx(exp(LN10*(catima::logEmin+4.0*step))).R(1e-9));
CHECK(etable[5]==approx(exp(LN10*(catima::logEmin+5.0*step))).R(1e-9));
CHECK(etable[catima::max_datapoints-1]==approx(exp(LN10*(catima::logEmax))).epsilon(1e-6));
} }
TEST_CASE("indexing"){ TEST_CASE("indexing"){
double val, dif; double val, dif;
catima::LogVArray<catima::max_datapoints> etable(catima::logEmin,catima::logEmax);
catima::EnergyTable<catima::max_datapoints> energy_table(catima::logEmin,catima::logEmax);
CHECK(EnergyTable_index(catima::energy_table, 0.0)==-1); CHECK(energy_table.index(0.0)==-1);
for(int i=0;i<catima::max_datapoints-1;i++){ for(int i=0;i<catima::max_datapoints-1;i++){
val = catima::energy_table.values[i]; val = energy_table[i];
dif = catima::energy_table.values[i+1] - val; dif = energy_table[i+1] - val;
CHECK(EnergyTable_index(catima::energy_table, val)==i); CHECK(energy_table.index(val)==i);
CHECK(EnergyTable_index(catima::energy_table, val+0.5*dif)==i); CHECK(energy_table.index(val+0.5*dif)==i);
CHECK(catima::energy_table.index(val)==i); CHECK(energy_table.index(val)==i);
CHECK(catima::energy_table.index(val+0.5*dif)==i); CHECK(energy_table.index(val+0.5*dif)==i);
CHECK(etable.index(val)==i);
CHECK(etable.index(val+0.5*dif)==i);
CHECK(etable.index(val+0.01*dif)==i);
CHECK(etable.index(val+0.99*dif)==i);
} }
} }