From 71d97e81c493b50a90feedd990843376678ef60e Mon Sep 17 00:00:00 2001 From: hrocho Date: Fri, 22 Apr 2022 22:52:10 +0200 Subject: [PATCH] energy table indexing --- CMakeLists.txt | 39 +++++++------ build_config.in | 1 + calculations.cpp | 43 +++++++------- calculations.h | 8 +++ catima.cpp | 10 ++-- cwrapper.cpp | 6 ++ pymodule/pycatima.cpp | 9 ++- pymodule/setup.py | 5 +- pymodule/setup.py.in | 22 ++++++++ spline.cpp | 10 +--- spline.h | 8 ++- storage.cpp | 4 ++ storage.h | 108 ++++++++++++++++++++++++++++-------- tests/test_calculations.cpp | 2 +- tests/test_storage.cpp | 50 +++++++++++------ 15 files changed, 224 insertions(+), 101 deletions(-) create mode 100644 pymodule/setup.py.in diff --git a/CMakeLists.txt b/CMakeLists.txt index f1d2611..1bdd288 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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_INTERPOLATION "use GSL inteRPOLATION" OFF) 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(PYTHON_WHEEL "make python wheel" OFF) ######## build type ############ if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug") 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") 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") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare") endif() @@ -56,18 +57,21 @@ endif() #endif(nurex_FOUND) find_package(fmt QUIET) -if(NOT fmt_FOUND) -message("fmt library not found, trying to dowload") - include(FetchContent) - FetchContent_Declare( - fmt - GIT_REPOSITORY https://github.com/fmtlib/fmt.git - GIT_TAG 8.1.1 - ) - set(CMAKE_POSITION_INDEPENDENT_CODE ON) - set(FMT_INSTALL ON) - FetchContent_MakeAvailable(fmt) -endif(NOT fmt_FOUND) +function(check_fmt) + if(NOT fmt_FOUND) + message("fmt library not found, trying to dowload") + include(FetchContent) + FetchContent_Declare( + fmt + GIT_REPOSITORY https://github.com/fmtlib/fmt.git + GIT_TAG 8.1.1 + ) + set(CMAKE_POSITION_INDEPENDENT_CODE ON) + set(FMT_INSTALL ON) + FetchContent_MakeAvailable(fmt) + endif(NOT fmt_FOUND) +endfunction(check_fmt) + configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in" "${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h") @@ -85,8 +89,9 @@ file(GLOB HEADERS *.h libs/*.h) add_library(catima ${SOURCES}) set_target_properties(catima PROPERTIES - POSITION_INDEPENDENT_CODE ON - LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib + POSITION_INDEPENDENT_CODE ON + LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib + ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib ) target_link_libraries(catima PUBLIC ${EXTRA_LIBS}) @@ -124,6 +129,7 @@ if(PYTHON_MODULE) FetchContent_MakeAvailable(pybind11) endif(NOT pybind11_FOUND) + check_fmt() #set(PYBIND11_CPP_STANDARD -std=c++14) pybind11_add_module(pycatima pymodule/pycatima.cpp) 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") if(PYTHON_WHEEL) + check_fmt() execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py bdist_wheel) endif(PYTHON_WHEEL) diff --git a/build_config.in b/build_config.in index 5d03161..001055d 100644 --- a/build_config.in +++ b/build_config.in @@ -8,5 +8,6 @@ #cmakedefine GLOBAL #cmakedefine REACTIONS #cmakedefine NUREX +#cmakedefine ET_CALCULATED_INDEX #endif diff --git a/calculations.cpp b/calculations.cpp index d2c3a54..cdbfda2 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -1,4 +1,4 @@ -#include +#include #include #include #include "catima/calculations.h" @@ -21,7 +21,7 @@ extern "C" namespace catima{ 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; 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 gamma = gamma_from_T(energy); 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; } @@ -485,7 +485,7 @@ double angular_scattering_power(const Projectile &p, const Target &t, double Es2 double lr = radiation_length(t.Z,t.A); //constexpr double Es2 = 198.81; //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){ @@ -496,7 +496,7 @@ double angular_scattering_power(const Projectile &p, const Material &mat, double double X0 = radiation_length(mat); //constexpr double Es2 = 198.81; //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){ @@ -507,10 +507,10 @@ double angular_scattering_power_xs(const Projectile &p, const Material &mat, dou double pv = _p*beta; double p1v1 = p1*beta1; 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 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 @@ -542,7 +542,7 @@ double radiation_length(int z, double m){ if(z==92){return 6.00;} 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 a2 = fine_structure*fine_structure*z2; 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 beta2 = beta*beta; 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) ){ - 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 = 24.89 * std::pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)* + log( 2.0*electron_mass*1e6*beta2/(33.05*std::pow(t.Z,1.6364))); cor = std::max(cor, 0.0 ); } //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){ - 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 B = 1.18-tz*(7.5e-03 - 4.53e-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){ @@ -705,7 +705,7 @@ double z_eff_Hubert(double pz, double E, double tz){ 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){ @@ -733,7 +733,7 @@ double z_eff_Winger(double pz, double beta, double tz){ if(tz==2){ tz = 2.8; } - x = beta /0.012 /pow(pz,0.45); + x = beta /0.012 /std::pow(pz,0.45); lnz =log(pz); lnzt=log(tz); @@ -752,13 +752,14 @@ double z_eff_global(double pz, double E, double tz){ else if(E<30.0 || pz<29){ return z_eff_Pierce_Blann(pz, E); } - else + else{ #ifdef GLOBAL return global_qmean(pz, tz, E); #else assert(false); return -1; #endif + } } 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 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); 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) ); } @@ -789,17 +790,17 @@ double z_eff_atima14(double pz, double T, double tz){ double c2 = 0.28; double c3 = 0.04; 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.; emin = 1000.; 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; } 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){ qwinger = z_eff_Winger(pz,beta,tz); qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin); diff --git a/calculations.h b/calculations.h index f61e3c8..72a91af 100644 --- a/calculations.h +++ b/calculations.h @@ -230,5 +230,13 @@ namespace catima{ inline double power(double x, double 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 diff --git a/catima.cpp b/catima.cpp index 9c4b0da..6550545 100644 --- a/catima.cpp +++ b/catima.cpp @@ -1,5 +1,5 @@ #include -#include +#include #include #include "catima/catima.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; double _p = p_from_T(p.T,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){ @@ -154,7 +154,7 @@ double angular_variance(Projectile p, const Material &t, const Config &c, int or auto fx0p = [&](double x)->double{ 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; if(c.scattering == scattering_types::dhighland){ 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{ 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); }; @@ -175,7 +175,7 @@ double angular_variance(Projectile p, const Material &t, const Config &c, int or if(c.scattering == scattering_types::gottschalk){ 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){ diff --git a/cwrapper.cpp b/cwrapper.cpp index e5d1e2e..8ed94a3 100644 --- a/cwrapper.cpp +++ b/cwrapper.cpp @@ -26,6 +26,7 @@ extern "C" { 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.scattering = 255; catima::Projectile p(pa,pz); catima::Material mat = make_material(ta,tz, thickness, density); 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){ catima::default_config.z_effective = catima_defaults.z_effective; + catima::default_config.scattering = 255; catima::Projectile p(pa,pz); catima::Material mat = make_material(ta,tz, thickness, density); 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){ catima::default_config.z_effective = catima_defaults.z_effective; + catima::default_config.scattering = 255; catima::Projectile p(pa,pz); catima::Material mat = make_material(ta,tz, 0, -1); return range(p, mat); @@ -60,12 +63,14 @@ extern "C" { 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.scattering = 255; catima::Projectile p(pa,pz); catima::Material mat = make_material(ta,tz, 0, -1); return range_straggling(p, T, mat); } 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::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){ + catima::default_config.scattering = 255; catima::Projectile p(pa,pz); catima::Material mat = make_material(ta,tz, 0, -1); return catima::energy_straggling_from_E(p,Tin,Tout,mat); diff --git a/pymodule/pycatima.cpp b/pymodule/pycatima.cpp index f2bfe7a..0cc231f 100644 --- a/pymodule/pycatima.cpp +++ b/pymodule/pycatima.cpp @@ -49,9 +49,14 @@ py::list storage_info(){ py::list get_energy_table(){ py::list r; - for(auto e : energy_table){ - r.append(e); + for (size_t i = 0; i < energy_table.size(); i++) + { + r.append(energy_table[i]); } + + //for(auto e : energy_table){ + //r.append(e); + //} return r; } diff --git a/pymodule/setup.py b/pymodule/setup.py index dee0b42..c8e67f2 100644 --- a/pymodule/setup.py +++ b/pymodule/setup.py @@ -5,10 +5,7 @@ from setuptools import setup DIR = Path(__file__).parents[0] -SRC = [str((DIR/'pycatima.cpp').resolve())]#+[str(fname.resolve()) for fname in DIR.glob('../*.cpp')] -#SRCG = [str(fname.resolve()) for fname in DIR.glob('../global/*.c')] -#SRC += SRCG -print (SRC) +SRC = [str((DIR/'pycatima.cpp').resolve())] example_module = Pybind11Extension( 'pycatima', SRC, diff --git a/pymodule/setup.py.in b/pymodule/setup.py.in new file mode 100644 index 0000000..afce630 --- /dev/null +++ b/pymodule/setup.py.in @@ -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}, +) diff --git a/spline.cpp b/spline.cpp index 76f4cf7..2c6871d 100644 --- a/spline.cpp +++ b/spline.cpp @@ -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" - namespace catima{ -} // namespace tk +} diff --git a/spline.h b/spline.h index 2abea56..153b1e6 100644 --- a/spline.h +++ b/spline.h @@ -111,7 +111,6 @@ struct cspline_special{ cspline_special() = default; const T *table; - const double *m_x; const double *m_y; std::array m_a,m_b,m_c; double m_b0, m_c0; @@ -120,7 +119,9 @@ struct cspline_special{ double operator()(double x)const{return evaluate(x);} 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 interpol; if(xindex(x), 0); double h=x-m_x[idx]; @@ -162,7 +164,7 @@ template cspline_special::cspline_special(const T &x, const std::vector& y, 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"); tridiagonal_matrix A{}; diff --git a/storage.cpp b/storage.cpp index 1372e7c..bcaf363 100644 --- a/storage.cpp +++ b/storage.cpp @@ -19,7 +19,11 @@ #include "catima/catima.h" namespace catima { Data _storage; + #ifdef VETABLE + LogVArray energy_table(logEmin,logEmax); + #else EnergyTable energy_table(logEmin,logEmax); + #endif bool operator==(const DataPoint &a, const DataPoint &b){ if( (a.m == b.m) && (a.p == b.p) && (a.config == b.config)){ diff --git a/storage.h b/storage.h index 3005536..93ce59e 100644 --- a/storage.h +++ b/storage.h @@ -29,7 +29,7 @@ #include "catima/spline.h" - +//#define VETABLE namespace catima{ /** @@ -48,36 +48,31 @@ namespace catima{ static constexpr int size() {return N;}; double values[N]; double step; - double* begin(){return values;} - double* end(){return &values[num];} - int index(double v)const noexcept{ - double lxval = (log(v/values[0])/LN10); + const double* begin()const{return values;} + const double* end()const{return &values[num];} + int index(double v)const noexcept{ 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 (std::floor(lxval/step)); if(v >= values[i+1]-numeric_epsilon)i++; // this is correction for floating point precision return i; + #else + auto it=std::upper_bound(begin(),end(),v); + return int(it-begin())-1; + #endif }; std::size_t num; }; - extern EnergyTable energy_table; - - template - int EnergyTable_index(const EnergyTable &table, double val){ - if(valtable.values[table.num-1])return -1; - 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; - } - - template - double EnergyTable_interpolate(const EnergyTable &table, double xval, double *y){ + template + double EnergyTable_interpolate(const T &table, double xval, double *y){ double r; - if(xvaltable.values[table.num-1])return 0.0; - if(xval==table.values[table.num-1])return y[table.num-1]; - int i = EnergyTable_index(table, xval); + if(xvaltable.values[table.size()-1])return 0.0; + if(xval==table.values[table.size()-1])return y[table.size()-1]; + int i = table.index(xval); double linstep = table.values[i+1] - table.values[i]; if(linstep == 0.0)return table.values[i]; double x = 1.0 - ((xval - table.values[i])/linstep); @@ -85,6 +80,63 @@ namespace catima{ return r; } + template + 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(N-1)-numeric_epsilon))return N-1; + double lxval = (log(v/value(0))/LN10); + int i = static_cast (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 + 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=max)return N-1; + assert(step>0.0); + return static_cast (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 energy_table; + #else + extern EnergyTable energy_table; + #endif + ////////////////////////////////////////////////////////////////////////////////////// #ifdef GSL_INTERPOLATION /// Interpolation class, to store interpolated values @@ -107,12 +159,13 @@ namespace catima{ }; #endif + template class InterpolatorCSpline{ public: - using xtype = EnergyTable; + //using xtype = EnergyTable; InterpolatorCSpline()=default; InterpolatorCSpline(const xtype &table, const std::vector &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 eval(double x)const{return ss.evaluate(x);} double derivative(double x)const{return ss.deriv(x);} @@ -128,7 +181,14 @@ namespace catima{ #ifdef GSL_INTERPOLATION using Interpolator = InterpolatorGSL; #else -using Interpolator = InterpolatorCSpline; +#ifdef VETABLE +//using Interpolator = InterpolatorSplineT>; +using Interpolator = InterpolatorCSpline>; +#else +//using Interpolator = InterpolatorSplineT>; +using Interpolator = InterpolatorCSpline>; +#endif + #endif #ifdef STORE_SPLINES diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 838d5e9..ee9cb30 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -425,7 +425,7 @@ using namespace std; CHECK(res1.dEdxo == res2.dEdxo); 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)); auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite); diff --git a/tests/test_storage.cpp b/tests/test_storage.cpp index 7446e96..3974c8f 100644 --- a/tests/test_storage.cpp +++ b/tests/test_storage.cpp @@ -114,27 +114,45 @@ using catima::LN10; } TEST_CASE("energy table"){ + catima::LogVArray etable(catima::logEmin,catima::logEmax); + catima::EnergyTable energy_table(catima::logEmin,catima::logEmax); double step = (catima::logEmax - catima::logEmin)/(catima::max_datapoints-1); - CHECK(catima::energy_table.step==step); - 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)); + CHECK(energy_table.step==step); + CHECK(energy_table[0]==approx(exp(LN10*(catima::logEmin))).R(1e-9)); + CHECK(energy_table[1]==approx(exp(LN10*(catima::logEmin+step))).R(1e-9)); + CHECK(energy_table[2]==approx(exp(LN10*(catima::logEmin+2.0*step))).R(1e-9)); + CHECK(energy_table[3]==approx(exp(LN10*(catima::logEmin+3.0*step))).R(1e-9)); + CHECK(energy_table[4]==approx(exp(LN10*(catima::logEmin+4.0*step))).R(1e-9)); + CHECK(energy_table[5]==approx(exp(LN10*(catima::logEmin+5.0*step))).R(1e-9)); + 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"){ double val, dif; + catima::LogVArray etable(catima::logEmin,catima::logEmax); + catima::EnergyTable 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