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

Merge pull request #90 from hrosiak/u2

Browse Source 
energy table indexing
This commit is contained in:
Andrej Prochazka 2022-04-22 22:56:20 +02:00 committed by GitHub
parent 910bf5a7d6 71d97e81c4
commit b88f5149b8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
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);
} }
} }