gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-23 02:38:51 -05:00
Code Issues Projects Releases Wiki Activity

update3

Browse Source
This commit is contained in:
hrocho 2018-10-21 22:08:16 +02:00
parent 3d240d2fd6
commit c971222984
20 changed files with 887 additions and 553 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
CMakeLists.txt
View File

@ -1,27 +1,30 @@
cmake_minimum_required(VERSION 3.2.0) cmake_minimum_required(VERSION 3.2.0)
project(catima) project(catima)
############ options ############# ############ options #############
#option(THREADS "Use multi-threading" ON) option(BUILD_SHARED_LIBS "build as shared library" ON)
option(PYTHON_MODULE "compile the Catima python module(requires numpy and cython installed)" OFF) option(PYTHON_MODULE "compile the Catima python module(requires numpy and cython installed)" OFF)
option(TESTS "build tests" OFF) option(TESTS "build tests" OFF)
option(EXAMPLES "build examples" ON) option(EXAMPLES "build examples" ON)
option(GSL_INTEGRATION "use GSL integration" ON) option(APPS "build catima applications" ON)
option(GENERATE_DATA "make data tables generator" OFF)
option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
option(DOCS "build documentation (requires doxygen)" OFF)
option(GLOBAL "build with global, sources are required" OFF) option(GLOBAL "build with global, sources are required" OFF)
option(REACTIONS "enable/disable nuclear reaction rate" ON) option(REACTIONS "enable/disable nuclear reaction rate" ON)
option(APPS "build catima applications" ON) option(STORE_SPLINES "store splines, if disables splines are always recreated" ON)
option(GSL_INTEGRATION "use GSL integration" OFF)
option(GSL_INTERPOLATION "use GSL inteRPOLATION" OFF)
option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
option(GENERATE_DATA "make data tables generator" OFF)
option(DOCS "build documentation (requires doxygen)" 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}")
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 "-Wall -Wextra -Wfatal-errors -g") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare")
endif() endif()
MESSAGE(STATUS "Build type Debug") MESSAGE(STATUS "Build type Debug")
endif() endif()
@ -41,18 +44,10 @@ else()
endif() endif()
############# Requirements ################## ############# Requirements ##################
find_package(GSL REQUIRED) if(GSL_INTEGRATION OR GSL_INTERPOLATION)
MESSAGE(STATUS "GSL include dirs: " ${GSL_INCLUDE_DIRS}) find_package(GSL REQUIRED)
#if(THREADS) MESSAGE(STATUS "GSL include dirs: " ${GSL_INCLUDE_DIRS})
# find_package(Threads REQUIRED) list(APPEND EXTRA_LIBS ${GSL_LIBRARIES} )
# set (EXTRA_LIBS ${EXTRA_LIBS} ${CMAKE_THREAD_LIBS_INIT})
# set (USE_THREADS ON)
# MESSAGE(STATUS "Nurex will use threads")
#endif(THREADS)
find_package(PythonInterp)
if(PYTHONINTERP_FOUND)
message("-- Python found: ${PYTHON_EXECUTABLE}")
endif() endif()
find_package(nurex QUIET) find_package(nurex QUIET)
@ -81,29 +76,19 @@ endif(GLOBAL)
file(GLOB HEADERS *.h) file(GLOB HEADERS *.h)
add_library(catima SHARED ${SOURCES}) add_library(catima SHARED ${SOURCES})
add_library(catima_static STATIC ${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
) )
set_target_properties(catima_static
PROPERTIES OUTPUT_NAME catima target_link_libraries(catima ${EXTRA_LIBS})
POSITION_INDEPENDENT_CODE ON
ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
)
target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES})
target_link_libraries(catima_static ${EXTRA_LIBS} ${GSL_LIBRARIES})
target_include_directories(catima target_include_directories(catima
PUBLIC $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> PUBLIC $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}> $<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
) )
target_include_directories(catima_static
PUBLIC $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
)
add_library(catima::catima ALIAS catima) add_library(catima::catima ALIAS catima)
add_library(catima::catima_static ALIAS catima_static)
FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima) FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima)
@ -111,10 +96,23 @@ FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima)
MESSAGE( STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER} ) MESSAGE( STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER} )
######## for python module ######## for python module
find_package(PythonInterp)
if(PYTHONINTERP_FOUND)
message("-- Python found: ${PYTHON_EXECUTABLE}")
endif()
if(PYTHON_MODULE) if(PYTHON_MODULE)
if(NOT PYTHONINTERP_FOUND) if(NOT PYTHONINTERP_FOUND)
MESSAGE(SEND_ERROR "Python is required to build nurex python modules") MESSAGE(SEND_ERROR "Python is required to build nurex python modules")
endif(NOT PYTHONINTERP_FOUND) endif(NOT PYTHONINTERP_FOUND)
#find_package(pybind11 REQUIRED)
#set(PYBIND11_CPP_STANDARD -std=c++14)
#pybind11_add_module(pycatima pymodule/pycatima)
#target_include_directories(pycatima PUBLIC
# $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include>
# $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/libs>
# $<INSTALL_INTERFACE:include>)
#target_link_libraries(pycatima PRIVATE catima)
find_program(CYTHON_EXECUTABLE find_program(CYTHON_EXECUTABLE
NAMES cython cython2 cython3 cython.bat NAMES cython cython2 cython3 cython.bat
DOC "path to the cython executable" DOC "path to the cython executable"
@ -182,7 +180,7 @@ endif(APPS)
####### install part ####### ####### install part #######
FILE(GLOB headers "*.h") FILE(GLOB headers "*.h")
include(GNUInstallDirs) include(GNUInstallDirs)
install (TARGETS catima catima_static install (TARGETS catima
EXPORT catimaConfig EXPORT catimaConfig
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}) ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
@ -194,7 +192,7 @@ install(EXPORT catimaConfig
DESTINATION lib/cmake/catima DESTINATION lib/cmake/catima
) )
export(TARGETS catima catima_static NAMESPACE catima:: FILE catimaConfig.cmake) export(TARGETS catima NAMESPACE catima:: FILE catimaConfig.cmake)
export(PACKAGE catima) export(PACKAGE catima)
###### packaging ####### ###### packaging #######

5
README.md
View File

@ -26,17 +26,18 @@ compile options, enable or disable with cmake:
> cmake ../ -D[OPTION] > cmake ../ -D[OPTION]
available options: available options:
* BUILD_SHARED_LIBS - if ON shared library is build, otherwise static
* PYTHON_MODULE - enable/disable building of the python bindigs, cython and numpy are required to build the catima python module, default OFF * PYTHON_MODULE - enable/disable building of the python bindigs, cython and numpy are required to build the catima python module, default OFF
* TESTS - build tests * TESTS - build tests
* EXAMPLES - build examples * EXAMPLES - build examples
* DOCS - prepare doxygen documentation (after cmake, __make docs__ needs to be executed) * DOCS - prepare doxygen documentation (after cmake, __make docs__ needs to be executed)
* GENERATE_DATA - makes program to re-generate precalculated tables (ie precalculated LS coefficients), default:OFF * GENERATE_DATA - makes program to re-generate precalculated tables (ie precalculated LS coefficients), default:OFF
* THIN_TARGET_APPROXIMATION - compile the library with thin target approximation, default: ON * THIN_TARGET_APPROXIMATION - compile the library with thin target approximation, default: ON
* GSL_INTEGRATION - use GSL integration functions, otherwise use built-in integrator, default: ON * GSL_INTEGRATION - use GSL integration functions, otherwise use built-in integrator, default: OFF
* GLOBAL - compile with GLOBAL code (source not included at the moment, needs to be manually added to __global__ directory, default:OFF) * GLOBAL - compile with GLOBAL code (source not included at the moment, needs to be manually added to __global__ directory, default:OFF)
ie: ie:
> cmake -DCATIMA_PYTHON=ON -DEXAMPLES=ON ../ > cmake -DPYTHON_MODULE=ON -DEXAMPLES=ON ../
after the compilation the libraries and headers must be either installed system-wide by make install or PATH and LD_LIBRARY_PATH must be adjusted to point to headers and library files. after the compilation the libraries and headers must be either installed system-wide by make install or PATH and LD_LIBRARY_PATH must be adjusted to point to headers and library files.

13
bin/catima_calculator.cpp
View File

@ -18,6 +18,19 @@ void help(){
std::cout<<"usage: catima_calculator config_file.json\n"; std::cout<<"usage: catima_calculator config_file.json\n";
} }
inline std::vector<double> linspace_vector(double a, double b, unsigned int num){
std::vector<double> res;
if(num>=2 && a<b){
res.resize(num);
double step = (b-a)/(num-1);
for(unsigned int i=0;i<(num-1);i++){
res[i]=a+(i*step);
}
res[num-1] = b;
}
return res;
}
json load_json(const char *fname); json load_json(const char *fname);
char* getCmdOption(char ** begin, char ** end, const std::string & option); char* getCmdOption(char ** begin, char ** end, const std::string & option);
Material json_material(json &j); Material json_material(json &j);

2
build_config.in
View File

@ -3,6 +3,8 @@
#cmakedefine THIN_TARGET_APPROXIMATION #cmakedefine THIN_TARGET_APPROXIMATION
#cmakedefine GSL_INTEGRATION #cmakedefine GSL_INTEGRATION
#cmakedefine GSL_INTERPOLATION
#cmakedefine STORE_SPLINES
#cmakedefine GLOBAL #cmakedefine GLOBAL
#cmakedefine REACTIONS #cmakedefine REACTIONS
#cmakedefine NUREX #cmakedefine NUREX

13
calculations.cpp
View File

@ -43,11 +43,17 @@ double dedx_n(const Projectile &p, const Target &t){
double asum = p.A + t.A; double asum = p.A + t.A;
double epsilon = 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV double epsilon = 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
double sn=0; double sn=0;
if(epsilon<=0){
if(epsilon<=30){ return 0.0;
}
else if(epsilon<=30){
assert(p.A>0);
assert(epsilon>0);
sn = log(1+(1.1383*epsilon))/ (2*(epsilon + 0.01321*pow(epsilon,0.21226) + 0.19593*pow(epsilon,0.5))); sn = log(1+(1.1383*epsilon))/ (2*(epsilon + 0.01321*pow(epsilon,0.21226) + 0.19593*pow(epsilon,0.5)));
} }
else{ else{
assert(p.A>0);
assert(epsilon>0);
sn = log(epsilon)/(2*epsilon); sn = log(epsilon)/(2*epsilon);
} }
sn = 100*8.4621*p.Z*t.Z*p.A*sn*Avogadro/(asum*zpowers*t.A); sn = 100*8.4621*p.Z*t.Z*p.A*sn*Avogadro/(asum*zpowers*t.A);
@ -798,6 +804,7 @@ double z_eff_global(double pz, double E, double tz){
#ifdef GLOBAL #ifdef GLOBAL
return global_qmean(pz, tz, E); return global_qmean(pz, tz, E);
#else #else
assert(false);
return -1; return -1;
#endif #endif
} }
@ -871,6 +878,8 @@ double z_eff_atima14(double pz, double T, double tz){
} }
} }
} }
#else
assert(false);
#endif #endif
return qmean; return qmean;
} }

6
calculations.h
View File

@ -39,9 +39,13 @@ namespace catima{
double reduced_energy_loss_unit(const Projectile &p, const Target &t); double reduced_energy_loss_unit(const Projectile &p, const Target &t);
/**
* @brief bethek_dedx_e - electronics stopping power
* @return stopping power
*/
double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0); double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config); double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config);
double bethek_barkas(double zp_eff,double eta, double zt); double bethek_barkas(double zp_eff,double eta, double zt);
double bethek_density_effect(double beta, int zt); double bethek_density_effect(double beta, int zt);

131
catima.cpp
View File

@ -17,7 +17,6 @@ namespace catima{
Config default_config; Config default_config;
bool operator==(const Config &a, const Config&b){ bool operator==(const Config &a, const Config&b){
if(std::memcmp(&a,&b,sizeof(Config)) == 0){ if(std::memcmp(&a,&b,sizeof(Config)) == 0){
return true; return true;
@ -30,11 +29,9 @@ bool operator==(const Config &a, const Config&b){
double dedx(Projectile &p, double T, const Material &mat, const Config &c){ double dedx(Projectile &p, double T, const Material &mat, const Config &c){
double sum = 0; double sum = 0;
if(T<=0)return 0.0; if(T<=0)return 0.0;
sum += dedx_n(p,mat);
double se=0;
p.T = T; p.T = T;
sum += dedx_n(p,mat);
double se=0;
if(p.T<=10){ if(p.T<=10){
se = sezi_dedx_e(p,mat); se = sezi_dedx_e(p,mat);
} }
@ -78,21 +75,23 @@ double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
double range(Projectile &p, double T, const Material &t, const Config &c){ double range(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return range_spline(T); return range_spline(T);
} }
double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c){ double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return p.A/range_spline.derivative(T); return p.A/range_spline.derivative(T);
} }
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
std::vector<double> dedx; std::vector<double> dedx;
dedx.reserve(T.size()); dedx.reserve(T.size());
for(auto e:T){ for(auto e:T){
@ -107,45 +106,52 @@ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T,
} }
double range_straggling(Projectile &p, double T, const Material &t, const Config &c){ double range_straggling(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
spline_type range_straggling_spline = get_range_straggling_spline(data);
return sqrt(range_straggling_spline(T)); return sqrt(range_straggling_spline(T));
} }
double range_variance(Projectile &p, double T, const Material &t, const Config &c){ double range_variance(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
spline_type range_straggling_spline = get_range_straggling_spline(data);
return range_straggling_spline(T); return range_straggling_spline(T);
} }
double domega2de(Projectile &p, double T, const Material &t, const Config &c){ double domega2de(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
spline_type range_straggling_spline = get_range_straggling_spline(data);
return range_straggling_spline.derivative(T); return range_straggling_spline.derivative(T);
} }
double da2de(Projectile &p, double T, const Material &t, const Config &c){ double da2de(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data);
return angular_variance_spline.derivative(T); return angular_variance_spline.derivative(T);
} }
double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){ double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
return sqrt(angular_straggling_spline(T) - angular_straggling_spline(Tout)); spline_type angular_variance_spline = get_angular_variance_spline(data);
return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
} }
double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){ double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
spline_type range_straggling_spline = get_range_straggling_spline(data);
double dEdxo = p.A/range_spline.derivative(Tout); double dEdxo = p.A/range_spline.derivative(Tout);
return dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(Tout))/p.A; return dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(Tout))/p.A;
} }
double energy_out(double T, double thickness, Interpolator &range_spline){ double energy_out(double T, double thickness, const Interpolator &range_spline){
constexpr double epsilon = 1E-5; constexpr double epsilon = 1E-5;
int counter = 0; int counter = 0;
double range; double range;
@ -172,14 +178,16 @@ double energy_out(double T, double thickness, Interpolator &range_spline){
} }
double energy_out(Projectile &p, double T, const Material &t, const Config &c){ double energy_out(Projectile &p, double T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return energy_out(T,t.thickness(),range_spline); return energy_out(T,t.thickness(),range_spline);
} }
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){ std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
std::vector<double> eout; std::vector<double> eout;
eout.reserve(T.size()); eout.reserve(T.size());
@ -199,15 +207,18 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
Result res; Result res;
double T = p.T; double T = p.T;
if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;} if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
auto data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
res.Ein = T; res.Ein = T;
res.range = range_spline(T); res.range = range_spline(T);
res.dEdxi = p.A/range_spline.derivative(T); res.dEdxi = p.A/range_spline.derivative(T);
res.Eout = energy_out(T,t.thickness(),range_spline); res.Eout = energy_out(T,t.thickness(),range_spline);
Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
spline_type range_straggling_spline = get_range_straggling_spline(data);
if(res.Eout<Ezero){ if(res.Eout<Ezero){
res.dEdxo = 0.0; res.dEdxo = 0.0;
@ -225,20 +236,23 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
double s2 = range_straggling_spline.derivative(res.Eout); double s2 = range_straggling_spline.derivative(res.Eout);
res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A; res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A;
Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data);
s1 = angular_variance_spline.derivative(T); s1 = angular_variance_spline.derivative(T);
s2 = angular_variance_spline.derivative(res.Eout); s2 = angular_variance_spline.derivative(res.Eout);
res.sigma_a = sqrt(0.5*(s1+s2)*edif); res.sigma_a = sqrt(0.5*(s1+s2)*edif);
} }
else{ else{
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A; res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data);
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout)); res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
} }
#else #else
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A; res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data);
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout)); res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
#endif #endif
if( !(c.skip&skip_tof) && t.thickness()>0){ if( !(c.skip&skip_tof) && t.thickness()>0){
@ -371,7 +385,7 @@ std::vector<double> calculate_tof(Projectile p, const Material &t, const Config
} }
return values; return values;
} }
/*
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
DataPoint dp(p,t,c); DataPoint dp(p,t,c);
dp.range.resize(max_datapoints); dp.range.resize(max_datapoints);
@ -411,6 +425,47 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
} }
return dp; return dp;
} }
*/
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
DataPoint dp(p,t,c);
dp.range.resize(max_datapoints);
dp.range_straggling.resize(max_datapoints);
dp.angular_variance.resize(max_datapoints);
auto fdedx = [&](double x)->double{
return 1.0/dedx(p,x,t,c);
};
auto fomega = [&](double x)->double{
//return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p,x,t,c),3);
};
double res=0.0;
//calculate 1st point to have i-1 element ready for loop
//res = integrator.integrate(fdedx,Ezero,energy_table(0));
//res = p.A*res;
//dp.range[0] = res;
dp.range[0] = 0.0;
dp.angular_variance[0] = 0.0;
//res = integrator.integrate(fomega,Ezero,energy_table(0));
//res = p.A*res;
dp.range_straggling[0]=0.0;
//p.T = energy_table(0);
for(int i=1;i<max_datapoints;i++){
res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
dp.range[i] = res + dp.range[i-1];
res = da2dx(p,energy_table(i),t)*res;
dp.angular_variance[i] = res + dp.angular_variance[i-1];
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
res = p.A*res;
dp.range_straggling[i] = res + dp.range_straggling[i-1];
}
return dp;
}
double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){
double res; double res;

3
catima.h
View File

@ -21,6 +21,7 @@
#include <vector> #include <vector>
// #define NDEBUG // #define NDEBUG
#include "catima/build_config.h"
#include "catima/config.h" #include "catima/config.h"
#include "catima/constants.h" #include "catima/constants.h"
#include "catima/structures.h" #include "catima/structures.h"
@ -141,7 +142,7 @@ namespace catima{
* @range_spline - precaclulated range spline for material * @range_spline - precaclulated range spline for material
* @return outcoming energy after the thickness in Mev/u * @return outcoming energy after the thickness in Mev/u
*/ */
double energy_out(double T, double thickness, Interpolator &range_spline); double energy_out(double T, double thickness, const Interpolator &range_spline);
/** /**
* calculates outcoming energy * calculates outcoming energy

2
catimac.pxd
View File

@ -141,4 +141,4 @@ cdef extern from "catima/storage.h" namespace "catima":
cdef EnergyTableType energy_table; cdef EnergyTableType energy_table;
cdef Data _storage; cdef Data _storage;
cdef DataPoint& get_data(const Projectile &p, const Material &t, Config c); cdef DataPoint& get_data(const Projectile &p, const Material &t, const Config c);

10
integrator.cpp
View File

@ -1,12 +1,13 @@
#include "integrator.h" #include "integrator.h"
//#ifdef GSL_INTEGRATION
#include "gsl/gsl_integration.h" #include "gsl/gsl_integration.h"
#include "gsl/gsl_errno.h" #include "gsl/gsl_errno.h"
//#endif
namespace catima{ namespace catima{
integrator_type integrator; integrator_type integrator;
IntegratorGSL integratorGSL(true); #ifdef GSL_INTEGRATION
double funcwrapper3(double x, void *_c){ double funcwrapper3(double x, void *_c){
std::function<double(double)> *f = (std::function<double(double)> *)_c; std::function<double(double)> *f = (std::function<double(double)> *)_c;
return (*f)(x); return (*f)(x);
@ -44,6 +45,5 @@ namespace catima{
} }
return result; return result;
}; };
#endif
} }

13
integrator.h
View File

@ -17,15 +17,20 @@
#ifndef INTEGRATOR_H #ifndef INTEGRATOR_H
#define INTEGRATOR_H #define INTEGRATOR_H
#include "catima/build_config.h" #include "catima/build_config.h"
#include "gsl/gsl_integration.h"
#include <functional> #include <functional>
#include <array> #include <array>
#ifdef USE_THREADS
#include <mutex> //#ifdef USE_THREADS
//#include <mutex>
//#endif
#ifdef GSL_INTEGRATION
#include "gsl/gsl_integration.h"
#endif #endif
namespace catima{ namespace catima{
#ifdef GSL_INTEGRATION
/// helper class to integrate functions using the GSL library /// helper class to integrate functions using the GSL library
class IntegratorGSL{ class IntegratorGSL{
public: public:
@ -45,6 +50,7 @@ class IntegratorGSL{
std::mutex integration_mutex; std::mutex integration_mutex;
#endif #endif
}; };
#endif
template<int order> template<int order>
struct GL_data{ struct GL_data{
@ -176,7 +182,6 @@ using integrator_type = GaussLegendreIntegration<8>;
#endif #endif
extern integrator_type integrator; extern integrator_type integrator;
extern IntegratorGSL integratorGSL;
} }
#endif #endif

4
reactions.cpp
View File

@ -23,8 +23,8 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
int ap = lround(projectile.A); int ap = lround(projectile.A);
int zp = lround(projectile.Z); int zp = lround(projectile.Z);
auto data = _storage.Get(projectile,target,c); auto& data = _storage.Get(projectile,target,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); spline_type range_spline = get_range_spline(data);
if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0; if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0;
auto sigma_r = [&](double th){ auto sigma_r = [&](double th){

241
spline.cpp
View File

@ -10,247 +10,6 @@
namespace catima{ namespace catima{
band_matrix::band_matrix(int dim)
{
resize(dim);
}
void band_matrix::resize(int dim)
{
assert(dim>0);
a.resize(dim);
d.resize(dim);
c.resize(dim);
}
int band_matrix::dim() const
{
return d.size();
}
// defines the new operator (), so that we can access the elements
// by A(i,j), index going from i=0,...,dim()-1
double & band_matrix::operator () (int i, int j)
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
assert(k<2 && k>-2);
if(k>0)return c[i];
else if(k==0) return d[i];
else return a[i];
}
double band_matrix::operator () (int i, int j) const
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
if(k>0)return c[i];
else if(k==0) return d[i];
else return a[i];
}
std::vector<double> band_matrix::trig_solve(const std::vector<double>& b) const
{
assert( this->dim()==(int)b.size() );
std::vector<double> x(this->dim());
std::vector<double> g(this->dim());
assert(d[0]!=0.0);
x[0] = b[0]/d[0];
double bet = d[0];
for(int j=1;j<this->dim();j++){
g[j] = c[j-1]/bet;
bet = d[j] - (a[j]*g[j]);
assert(bet != 0.0);
x[j] = (b[j]-a[j]*x[j-1])/bet;
}
for(int j=this->dim()-2;j>=0;j--){
x[j] -= g[j+1]*x[j+1];
}
return x;
}
// spline implementation
// -----------------------
void spline::set_boundary(spline::bd_type left, double left_value,
spline::bd_type right, double right_value,
bool force_linear_extrapolation)
{
assert(n==0); // set_points() must not have happened yet
m_left=left;
m_right=right;
m_left_value=left_value;
m_right_value=right_value;
m_force_linear_extrapolation=force_linear_extrapolation;
}
void spline::set_points(const double *x,
const double *y,
const size_t num
)
{
assert(num>2);
m_x=x;
m_y=y;
n=num;
// TODO: maybe sort x and y, rather than returning an error
for(int i=0; i<n-1; i++) {
assert(m_x[i]<m_x[i+1]);
}
// setting up the matrix and right hand side of the equation system
// for the parameters b[]
band_matrix A(n);
std::vector<double> rhs(n);
for(int i=1; i<n-1; i++) {
A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
}
// boundary conditions
if(m_left == spline::bd_type::second_deriv) {
// 2*b[0] = f''
A(0,0)=2.0;
A(0,1)=0.0;
rhs[0]=m_left_value;
} else{
// c[0] = f', needs to be re-expressed in terms of b:
// (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
A(0,0)=2.0*(x[1]-x[0]);
A(0,1)=1.0*(x[1]-x[0]);
rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value);
}
if(m_right == spline::bd_type::second_deriv) {
// 2*b[n-1] = f''
A(n-1,n-1)=2.0;
A(n-1,n-2)=0.0;
rhs[n-1]=m_right_value;
} else{
// c[n-1] = f', needs to be re-expressed in terms of b:
// (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
// = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
A(n-1,n-1)=2.0*(x[n-1]-x[n-2]);
A(n-1,n-2)=1.0*(x[n-1]-x[n-2]);
rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
// solve the equation system to obtain the parameters b[]
//m_b=A.lu_solve(rhs);
m_b=A.trig_solve(rhs);
// calculate parameters a[] and c[] based on b[]
m_a.resize(n);
m_c.resize(n);
for(int i=0; i<n-1; i++) {
m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
- 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
}
// for left extrapolation coefficients
m_b0 = (m_force_linear_extrapolation==false) ? m_b[0] : 0.0;
m_c0 = m_c[0];
// for the right extrapolation coefficients
// f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
double h=x[n-1]-x[n-2];
// m_b[n-1] is determined by the boundary condition
m_a[n-1]=0.0;
m_c[n-1]=3.0*m_a[n-2]*h*h+2.0*m_b[n-2]*h+m_c[n-2]; // = f'_{n-2}(x_{n-1})
if(m_force_linear_extrapolation==true)
m_b[n-1]=0.0;
}
double spline::operator() (double x) const
{
assert(n>2);
// find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
auto it=std::lower_bound(m_x,m_x+n,x);
//int idx=std::max( int(it-m_x)-1, 0);
if(it!=m_x)it--;
int idx = std::distance(m_x,it);
double mx = *it;
double h=x-mx;
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
interpol=(m_b0*h + m_c0)*h + m_y[0];
} else if(x>m_x[n-1]) {
// extrapolation to the right
interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1];
} else {
// interpolation
interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
}
return interpol;
}
double spline::deriv(int order, double x) const
{
assert(order>0);
// find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
auto it=std::lower_bound(m_x,m_x+n,x);
//int idx=std::max( int(it-m_x)-1, 0);
if(it!=m_x)it--;
int idx = std::distance(m_x,it);
double mx = *it;
double h=x-mx;
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
switch(order) {
case 1:
interpol=2.0*m_b0*h + m_c0;
break;
case 2:
interpol=2.0*m_b0*h;
break;
default:
interpol=0.0;
break;
}
} else if(x>m_x[n-1]) {
// extrapolation to the right
switch(order) {
case 1:
interpol=2.0*m_b[n-1]*h + m_c[n-1];
break;
case 2:
interpol=2.0*m_b[n-1];
break;
default:
interpol=0.0;
break;
}
} else {
// interpolation
switch(order) {
case 1:
interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
break;
case 2:
interpol=6.0*m_a[idx]*h + 2.0*m_b[idx];
break;
case 3:
interpol=6.0*m_a[idx];
break;
default:
interpol=0.0;
break;
}
}
return interpol;
}
} // namespace tk } // namespace tk

233
spline.h
View File

@ -17,77 +17,208 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef TK_SPLINE_H #ifndef CATIMA_SPLINE_H
#define TK_SPLINE_H #define CATIMA_SPLINE_H
#include <cstdio> #include <cstdio>
#include <cassert> #include <cassert>
#include <vector> #include <vector>
#include <algorithm> #include <algorithm>
#include <array>
#include "catima/constants.h"
#ifdef GSL_INTERPOLATION
#include <gsl/gsl_spline.h>
#endif
namespace catima namespace catima
{ {
// band matrix solver enum interpolation_t {cspline, linear};
class band_matrix
/**
* Tridiagonal matrix solver
*/
template<int N>
class tridiagonal_matrix
{ {
private: private:
std::vector<double> a; std::array<double,N> a;
std::vector<double> d; std::array<double,N> d;
std::vector<double> c; std::array<double,N> c;
std::vector<double> save;
public: public:
band_matrix() {}; // constructor tridiagonal_matrix() {}
band_matrix(int dim); // constructor
~band_matrix() {}; // destructor
void resize(int dim); // init with dim,n_u,n_l
int dim() const; // matrix dimension
// access operator // access operator
double & operator () (int i, int j); // write double & operator () (unsigned int i, unsigned int j); // write
double operator () (int i, int j) const; // read double operator () (unsigned int i, unsigned int j) const; // read
// we can store an additional diogonal (in m_lower) std::array<double, N> trig_solve(const std::array<double, N>& b) const;
std::vector<double> trig_solve(const std::vector<double>& b) const;
}; };
template<int N>
// spline interpolation double & tridiagonal_matrix<N>::operator () (unsigned int i, unsigned int j)
class spline
{ {
public: int k=j-i;
enum class bd_type { if(k == -1)return c[i];
first_deriv = 1, else if(k==0) return d[i];
second_deriv = 2 else return a[i];
}; }
private: template<int N>
const double *m_x, *m_y; // x,y coordinates of points double tridiagonal_matrix<N>::operator () (unsigned int i, unsigned int j) const
size_t n=0; {
// interpolation parameters int k=j-i;
// f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i if(k==-1)return c[i];
std::vector<double> m_a,m_b,m_c; // spline coefficients else if(k==0) return d[i];
double m_b0, m_c0; // for left extrapol else if(k==1)return a[i];
bd_type m_left = bd_type::second_deriv; else return 0.0;
bd_type m_right = bd_type::second_deriv; }
double m_left_value = 0.0;
double m_right_value = 0.0;
bool m_force_linear_extrapolation = false;
public: template<int N>
// set default boundary condition to be zero curvature at both ends std::array<double, N> tridiagonal_matrix<N>::trig_solve(const std::array<double, N>& b) const
spline(){} {
std::array<double, N> x;
if(d[0] == 0.0){return x;}
std::array<double, N> g;
x[0] = b[0]/d[0];
double bet = d[0];
for(std::size_t j=1, max=N;j<max;j++){
g[j] = c[j-1]/bet;
bet = d[j] - (a[j]*g[j]);
if(bet == 0.0){
x.fill(0.0);
return x;
}
x[j] = (b[j]-a[j]*x[j-1])/bet;
}
for(int j=N-2;j>=0;j--){
x[j] -= g[j+1]*x[j+1];
}
return x;
}
// optional, but if called it has to come be before set_points()
void set_boundary(bd_type left, double left_value, /**
bd_type right, double right_value, * Cubic Spline class, accepting EnergyTable type as x-variable
bool force_linear_extrapolation=false); */
void set_points(const double *x, template<typename T>
const double *y, struct cspline_special{
const size_t num); constexpr static int N = T::size();
double operator() (double x) const; cspline_special(const T& x,
double deriv(int order, double x) const; const std::vector<double>& y,
bool boundary_second_deriv = true);
cspline_special() = default;
const T *table;
const double *m_x;
const double *m_y;
std::array<double,N> m_a,m_b,m_c;
double m_b0, m_c0;
double operator()(double x)const{return evaluate(x);}
double evaluate(double x) const
{
int idx=std::max( table->index(x), 0);
double h=x-m_x[idx];
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
interpol=(m_b0*h + m_c0)*h + m_y[0];
} else if(x>m_x[N-1]) {
// extrapolation to the right
interpol=(m_b[N-1]*h + m_c[N-1])*h + m_y[N-1];
} else {
// interpolation
interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
}
return interpol;
}
double deriv(double x) const
{
int idx=std::max( table->index(x), 0);
double h=x-m_x[idx];
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
interpol=2.0*m_b0*h + m_c0;
} else if(x>m_x[N-1]) {
// extrapolation to the right
interpol=2.0*m_b[N-1]*h + m_c[N-1];
} else {
// interpolation
interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
}
return interpol;
}
static_assert (T::size()>2, "N must be > 2");
}; };
} // namespace tk template<typename T>
cspline_special<T>::cspline_special(const T &x,
const std::vector<double>& y,
bool boundary_second_deriv
):table(&x),m_y(y.data()),m_x(x.values)
{
static_assert (N>2, "N must be > 2");
tridiagonal_matrix<N> A{};
std::array<double, N> rhs;
for(std::size_t i=1; i<N-1; i++) {
A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
}
// boundary conditions
if(boundary_second_deriv) {
// 2*b[0] = f''
A(0,0)=2.0;
A(0,1)=0.0;
rhs[0]=0.0; // 0.0 is value of derivative
A(N-1,N-1)=2.0;
A(N-1,N-2)=0.0;
rhs[N-1]=0.0; // 0.0 is value of derivative
} else {
// c[0] = f', needs to be re-expressed in terms of b:
// (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
A(0,0)=2.0*(x[1]-x[0]);
A(0,1)=1.0*(x[1]-x[0]);
rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-0.0); // 0.0 is deriv value
#endif /* TK_SPLINE_H */ // c[n-1] = f', needs to be re-expressed in terms of b:
// (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
// = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
A(N-1,N-1)=2.0*(x[N-1]-x[N-2]);
A(N-1,N-2)=1.0*(x[N-1]-x[N-2]);
rhs[N-1]=3.0*(0.0-(y[N-1]-y[N-2])/(x[N-1]-x[N-2]));
}
m_b=A.trig_solve(rhs);
// calculate parameters a[] and c[] based on b[]
for(int i=0; i<N-1; i++) {
m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
- 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
}
// for left extrapolation coefficients
//s.m_b0 = (m_force_linear_extrapolation==false) ? s.m_b[0] : 0.0;
m_b0 = 0.0;
m_c0 = m_c[0];
double h=x[N-1]-x[N-2];
m_a[N-1]=0.0;
m_c[N-1]=3.0*m_a[N-2]*h*h+2.0*m_b[N-2]*h+m_c[N-2]; // = f'_{n-2}(x_{n-1})
m_b[N-1]=0.0;
}
} // namespace end
#endif

142
storage.cpp
View File

@ -30,10 +30,71 @@ namespace catima {
} }
} }
DataPoint::~DataPoint(){ #ifdef GSL_INTERPOLATION
//////////// Interpolator ////////////////////////////////
InterpolatorGSL::InterpolatorGSL(const EnergyTable<max_datapoints>& x, const std::vector<double>& y, interpolation_t type){
acc = gsl_interp_accel_alloc ();
const int num = y.size();
if(type==cspline)
spline = gsl_spline_alloc (gsl_interp_cspline, num);
else
spline = gsl_spline_alloc (gsl_interp_linear, num);
gsl_spline_init (spline, x.values, y.data(), num);
min= x[0];
max= x[num-1];
}
InterpolatorGSL::~InterpolatorGSL(){
gsl_interp_accel_free (acc);
gsl_spline_free (spline);
}
double InterpolatorGSL::eval(double x) const{
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval(spline, x, acc);
}
double InterpolatorGSL::derivative(double x)const{
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval_deriv (spline, x, acc);
}
#endif
#ifdef STORE_SPLINES
const Interpolator& get_range_spline(const DataPoint &data){
return data.range_spline;
} }
const Interpolator& get_range_straggling_spline(const DataPoint &data){
return data.range_straggling_spline;
}
const Interpolator& get_angular_variance_spline(const DataPoint &data){
return data.angular_variance_spline;
}
#else
Interpolator get_range_spline(const DataPoint &data){
//return Interpolator(energy_table.values,data.range);
//return data.range_spline;
return Interpolator(energy_table,data.range);
}
Interpolator get_range_straggling_spline(const DataPoint &data){
//return Interpolator(energy_table.values,data.range_straggling);
//return data.range_straggling_spline;
return Interpolator(energy_table,data.range_straggling);
}
Interpolator get_angular_variance_spline(const DataPoint &data){
//return Interpolator(energy_table.values,data.angular_variance);
//return data.angular_variance_spline;
return Interpolator(energy_table,data.angular_variance);
}
#endif
Data::Data(){ Data::Data(){
//storage.reserve(max_storage_data); // disabled because of "circular" storage //storage.reserve(max_storage_data); // disabled because of "circular" storage
storage.resize(max_storage_data); storage.resize(max_storage_data);
@ -57,12 +118,20 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){
index->angular_variance = calculate_angular_variance(p,t,c); index->angular_variance = calculate_angular_variance(p,t,c);
#else #else
*index = calculate_DataPoint(p,t,c); *index = calculate_DataPoint(p,t,c);
#endif #ifdef STORE_SPLINES
//index->range_spline = Interpolator(energy_table.values,index->range);
//index->range_straggling_spline = Interpolator(energy_table.values,index->range_straggling);
//index->angular_variance_spline = Interpolator(energy_table.values,index->angular_variance);
index->range_spline = Interpolator(energy_table, index->range);
index->range_straggling_spline = Interpolator(energy_table, index->range_straggling);
index->angular_variance_spline = Interpolator(energy_table, index->angular_variance);
#endif
#endif
index++; index++;
} }
DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
for(auto &e:storage){ for(auto &e:storage){
if( (e.p==p) && (e.m==t) && (e.config==c)){ if( (e.p==p) && (e.m==t) && (e.config==c)){
return e; return e;
@ -73,71 +142,4 @@ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
return *std::prev(index); return *std::prev(index);
} }
//////////// Interpolator ////////////////////////////////
InterpolatorGSL::InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type){
acc = gsl_interp_accel_alloc ();
if(type==cspline)
spline = gsl_spline_alloc (gsl_interp_cspline, num);
else
spline = gsl_spline_alloc (gsl_interp_linear, num);
gsl_spline_init (spline, x, y, num);
min= x[0];
max= x[num-1];
}
InterpolatorGSL::InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
//Interpolator(x.data(),y.data(),x.size());
acc = gsl_interp_accel_alloc ();
if(type==cspline)
spline = gsl_spline_alloc (gsl_interp_cspline, x.size());
else
spline = gsl_spline_alloc (gsl_interp_linear, x.size());
gsl_spline_init (spline, x.data(), y.data(), x.size());
min= x[0];
max= x[x.size()-1];
}
InterpolatorGSL::~InterpolatorGSL(){
gsl_interp_accel_free (acc);
gsl_spline_free (spline);
}
double InterpolatorGSL::eval(double x){
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval(spline, x, acc);
}
double InterpolatorGSL::derivative(double x){
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval_deriv (spline, x, acc);
}
//////////// Interpolator2 ////////////////////////////////
#ifdef BUILTIN_SPLINE
Interpolator2::Interpolator2(const double *x, const double *y, int num){
ss.set_points(x,y,num);
min= x[0];
max= x[num-1];
}
double Interpolator2::eval(double x){
if(x<min)x=min;
if(x>max)x=max;
return ss(x);
}
double Interpolator2::derivative(double x){
if(x<min)x=min;
if(x>max)x=max;
return ss.deriv(1,x);
}
#endif
} }

212
storage.h
View File

@ -18,24 +18,23 @@
#define STORAGE #define STORAGE
#include <vector> #include <vector>
#include <array>
#include <iterator> #include <iterator>
#include <cmath> #include <cmath>
//#include <unordered_set> //#include <unordered_set>
#include <gsl/gsl_spline.h>
#include "catima/build_config.h" #include "catima/build_config.h"
#include "catima/constants.h" #include "catima/constants.h"
#include "catima/structures.h" #include "catima/structures.h"
#include "catima/config.h" #include "catima/config.h"
#ifdef BUILTIN_SPLINE
#include "catima/spline.h" #include "catima/spline.h"
#endif
namespace catima{ namespace catima{
enum interpolation_t {cspline, linear}; /**
* Class to store energy points, log spaced from logmin to logmax.
*/
template<int N> template<int N>
struct EnergyTable{ struct EnergyTable{
EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){ EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){
@ -45,10 +44,18 @@ namespace catima{
} }
} }
double operator()(int i)const{return values[i];} double operator()(int i)const{return values[i];}
double operator[](int i)const{return values[i];}
static constexpr int size() {return N;};
double values[N]; double values[N];
double step; double step;
double* begin(){return values;} double* begin(){return values;}
double* end(){return &values[num-1];} double* end(){return &values[num];}
int index(double v)const noexcept{
double lxval = (log(v/values[0])/M_LN10);
if(v<values[0] || step==0.0)return -1;
if(v>=values[N-1])return N-1;
return static_cast<int> (std::floor(lxval/step));
};
std::size_t num; std::size_t num;
}; };
@ -75,75 +82,19 @@ namespace catima{
return r; return r;
} }
/* //////////////////////////////////////////////////////////////////////////////////////
template<int N> #ifdef GSL_INTERPOLATION
struct EnergyTableLinear{ /// Interpolation class, to store interpolated values
constexpr EnergyTableLinear():values(),num(N){
for(auto i=0;i<N;i++){
values[i]=exp(M_LN10*(logEmin + ((double)i)*(logEmax-logEmin)/(N - 1.0)));
}
}
double operator()(int i)const{return values[i];}
double values[N];
std::size_t num;
};
*/
// return vector with lineary spaced elements from a to b, num is number of elements
inline std::vector<double> linspace_vector(double a, double b, unsigned int num){
std::vector<double> res;
if(num>=2 && a<b){
res.resize(num);
double step = (b-a)/(num-1);
for(unsigned int i=0;i<(num-1);i++){
res[i]=a+(i*step);
}
res[num-1] = b;
}
return res;
}
class DataPoint{
public:
Projectile p;
Material m;
Config config;
std::vector<double> range;
std::vector<double> range_straggling;
std::vector<double> angular_variance;
DataPoint(){};
DataPoint(const Projectile _p, const Material _m,const Config &_c=default_config):p(_p),m(_m),config(_c){};
~DataPoint();
friend bool operator==(const DataPoint &a, const DataPoint &b);
};
class Data{
public:
Data();
~Data();
void Add(const Projectile &p, const Material &t, const Config &c=default_config);
int GetN() const {return storage.size();};
void Reset(){storage.clear();storage.resize(max_storage_data);index=storage.begin();};
DataPoint& Get(const Projectile &p, const Material &t, const Config &c=default_config);
DataPoint& Get(unsigned int i){return storage[i];};
int get_index() {return std::distance(storage.begin(),index);}
private:
std::vector<DataPoint> storage;
std::vector<DataPoint>::iterator index;
};
/// Interpolation class, to store interpolated values
class InterpolatorGSL{ class InterpolatorGSL{
public: public:
InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type=cspline); InterpolatorGSL(){};
InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline); InterpolatorGSL(const EnergyTable<max_datapoints>& x, const std::vector<double>& y, interpolation_t type=cspline);
~InterpolatorGSL(); ~InterpolatorGSL();
double operator()(double x){return eval(x);}; double operator()(double x)const{return eval(x);};
double eval(double x); double eval(double x) const;
double derivative(double x); double derivative(double x) const;
double get_min(){return min;}; double get_min()const{return min;};
double get_max(){return max;}; double get_max()const{return max;};
private: private:
double min=0; double min=0;
@ -151,33 +102,126 @@ inline std::vector<double> linspace_vector(double a, double b, unsigned int num)
gsl_interp_accel *acc; gsl_interp_accel *acc;
gsl_spline *spline; gsl_spline *spline;
}; };
#endif
#ifdef BUILTIN_SPLINE class InterpolatorCSpline{
class Interpolator2{
public: public:
Interpolator2(const double *x, const double *y, int num); using xtype = EnergyTable<max_datapoints>;
double operator()(double x){return eval(x);}; InterpolatorCSpline()=default;
double eval(double x); InterpolatorCSpline(const xtype &table, const std::vector<double> &y):
double derivative(double x); min(table.values[0]), max(table.values[max_datapoints-1]), ss(table,y){}
double get_min(){return min;}; double operator()(double x)const{return eval(x);}
double get_max(){return max;}; double eval(double x)const{return ss.evaluate(x);}
double derivative(double x)const{return ss.deriv(x);}
double get_min()const{return min;}
double get_max()const{return max;}
private: private:
double min=0; double min=0;
double max=0; double max=0;
spline ss; cspline_special<xtype> ss;
}; };
#ifdef GSL_INTERPOLATION
using Interpolator = InterpolatorGSL;
#else
using Interpolator = InterpolatorCSpline;
#endif #endif
#ifdef STORE_SPLINES
using spline_type = const Interpolator&;
#else
using spline_type = Interpolator;
#endif
// return vector with lineary spaced elements from a to b, num is number of elements
/**
* @brief structure to store calculated data points and optionally also splines
*/
class DataPoint{
public:
Projectile p;
Material m;
Config config;
std::vector<double> range;
std::vector<double> range_straggling;
std::vector<double> angular_variance;
#ifdef STORE_SPLINES
Interpolator range_spline;
Interpolator range_straggling_spline;
Interpolator angular_variance_spline;
#endif
DataPoint()=default;
DataPoint(const Projectile _p, const Material _m,const Config &_c=default_config):p(_p),m(_m),config(_c){}
DataPoint(const DataPoint&)=delete;
DataPoint(DataPoint&&)=default;
DataPoint& operator=(const DataPoint&)=default;
DataPoint& operator=(DataPoint&&)=default;
friend bool operator==(const DataPoint &a, const DataPoint &b);
};
#ifdef STORE_SPLINES
const Interpolator& get_range_spline(const DataPoint &data);
const Interpolator& get_range_straggling_spline(const DataPoint &data);
const Interpolator& get_angular_variance_spline(const DataPoint &data);
#else
Interpolator get_range_spline(const DataPoint &data);
Interpolator get_range_straggling_spline(const DataPoint &data);
Interpolator get_angular_variance_spline(const DataPoint &data);
#endif
/**
* @brief The Data class to store DataPoints
*/
class Data{
public:
Data();
~Data();
/**
* @brief Add new DataPoint
* @param p - Projectile
* @param t - Material
* @param c - Config
*/
void Add(const Projectile &p, const Material &t, const Config &c=default_config);
int GetN() const {return storage.size();};
void Reset(){storage.clear();storage.resize(max_storage_data);index=storage.begin();};
/**
* @brief Get DataPoint reference for projectile-target-config combination
* @param p - Projectile
* @param t - Material
* @param c - Config
* @return reference to DataPoint
*/
DataPoint& Get(const Projectile &p, const Material &t, const Config &c=default_config);
DataPoint& Get(unsigned int i){return storage[i];};
int get_index() {return std::distance(storage.begin(),index);}
private:
std::vector<DataPoint> storage;
std::vector<DataPoint>::iterator index;
};
extern Data _storage; extern Data _storage;
inline DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){ /**
* @brief get_data - Get DataPoint from the global storage class
* @param p - Projectile
* @param t - Material
* @param c - Config
* @return const reference to DataPoint
*/
inline const DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
return _storage.Get(p, t, c); return _storage.Get(p, t, c);
} }
bool operator==(const DataPoint &a, const DataPoint &b); bool operator==(const DataPoint &a, const DataPoint &b);
using InterpolatorLinear = InterpolatorGSL;
using Interpolator = InterpolatorGSL;
} }
#endif #endif

1
tests/test.py
View File

@ -217,6 +217,7 @@ class TestStructures(unittest.TestCase):
energies = [100,500,1000] energies = [100,500,1000]
res2 = catima.dedx_from_range(p,graphite,energy=energies) res2 = catima.dedx_from_range(p,graphite,energy=energies)
self.assertEqual(len(res2),len(energies)) self.assertEqual(len(res2),len(energies))
self.assertEqual(len(res2),3)
for i,e in enumerate(energies): for i,e in enumerate(energies):
r = catima.dedx_from_range(p, graphite, energy=e) r = catima.dedx_from_range(p, graphite, energy=e)
self.assertAlmostEqual(res2[i], r, 0.1) self.assertAlmostEqual(res2[i], r, 0.1)

309
tests/test2.py Normal file
View File

@ -0,0 +1,309 @@
import sys
sys.path.insert(0,"../build")
import unittest
import pycatima as catima
import math
class TestStructures(unittest.TestCase):
def test_Projectile(self):
print(catima.storage_info())
p = catima.Projectile(238,92)
self.assertEqual(p.A(),238)
self.assertEqual(p.Z(),92)
self.assertEqual(p.Q(),92)
p = catima.Projectile(238,92,90)
self.assertEqual(p.A(),238)
self.assertEqual(p.Z(),92)
self.assertEqual(p.Q(),90)
p.T(1000)
self.assertEqual(p.T(),1000)
p(500)
self.assertEqual(p.T(),500)
p = catima.Projectile(238,92,90, T=100)
self.assertEqual(p.T(),100)
def test_Material(self):
mat = catima.Material()
mat.add_element(12,6,1)
self.assertEqual(mat.ncomponents(),1)
mat.add_element(1,1,2)
self.assertEqual(mat.ncomponents(),2)
mat2 = catima.Material([[12.01,6,1]])
self.assertEqual(mat2.ncomponents(),1)
self.assertEqual(mat2.molar_mass(),12.01)
mat3 = catima.Material([[12,6,1]])
self.assertEqual(mat3.ncomponents(),1)
self.assertEqual(mat3.molar_mass(),12)
Water = catima.Material([[1,1,2],[16,8,1]])
self.assertEqual(Water.molar_mass(),18)
mat2 = catima.Material([[0,6,1]])
self.assertEqual(mat2.ncomponents(),1)
self.assertAlmostEqual(mat2.molar_mass(),12,1)
mat5 = catima.Material([[0,6,1]],density=1.9, thickness=0.5)
self.assertEqual(mat5.ncomponents(),1)
self.assertEqual(mat5.thickness(),0.5)
self.assertEqual(mat5.density(),1.9)
mat6 = catima.Material([[0,6,1]],density=1.9, thickness=0.5,i_potential=80.0)
self.assertEqual(mat6.ncomponents(),1)
self.assertEqual(mat6.thickness(),0.5)
self.assertEqual(mat6.density(),1.9)
self.assertEqual(mat6.I(),80.0)
# copy
mat3.density(1.8)
matc = mat3.copy()
self.assertEqual(matc.ncomponents(),1)
self.assertEqual(matc.molar_mass(),12)
self.assertEqual(matc.density(),1.8)
mat3.density(2.0)
self.assertEqual(matc.density(),1.8)
self.assertEqual(mat3.density(),2.0)
def test_default_material(self):
m1 = catima.get_material(6);
self.assertAlmostEqual(m1.molar_mass(),12,1)
self.assertEqual(m1.ncomponents(),1)
self.assertAlmostEqual(m1.density(),2.0,1)
m2 = catima.get_material(catima.material.Water)
self.assertEqual(m2.ncomponents(),2)
self.assertAlmostEqual(m2.molar_mass(),18,1)
self.assertAlmostEqual(m2.density(),1.0,1)
m3 = catima.get_material(3001)
self.assertEqual(m3.ncomponents(),0)
self.assertAlmostEqual(m3.molar_mass(),0,1)
self.assertAlmostEqual(m3.density(),0.0,1)
def test_layers(self):
graphite = catima.get_material(6)
graphite.thickness(0.5)
p10 = catima.get_material(catima.material.P10)
p10.thickness(0.01)
n2 = catima.get_material(7)
n2.thickness(0.02)
mat= catima.Layers()
self.assertEqual(mat.num(),0)
mat.add(graphite)
self.assertEqual(mat.num(),1)
self.assertAlmostEqual(mat[0].molar_mass(),12,1)
self.assertAlmostEqual(mat[0].thickness(),0.5,1)
self.assertAlmostEqual(mat[0].density(),2.0,1)
mat.add(p10)
self.assertEqual(mat.num(),2)
graphite.thickness(1.0)
graphite.density(1.8)
mat.add(graphite)
self.assertEqual(mat.num(),3)
self.assertAlmostEqual(mat[2].molar_mass(),12,1)
self.assertAlmostEqual(mat[0].thickness(),0.5,1)
self.assertAlmostEqual(mat[0].density(),2.0,1)
self.assertAlmostEqual(mat[2].thickness(),1.0,1)
self.assertAlmostEqual(mat[2].density(),1.8,1)
mat[2].thickness(1.2)
mat[2].density(1.9)
self.assertAlmostEqual(mat.materials[2].thickness(),1.2,1)
self.assertAlmostEqual(mat.materials[2].density(),1.9,1)
#self.assertAlmostEqual(mat.materials[0].thickness(),0.5,1)
#self.assertAlmostEqual(mat.materials[0].density(),2.0,1)
self.assertEqual(mat[3],None)
self.assertEqual(mat["a"],None)
mat2 = catima.Layers()
mat2.add(n2)
self.assertEqual(mat2.num(),1)
mats = mat2 + mat
self.assertEqual(mats.num(),4)
self.assertAlmostEqual(mats[0].molar_mass(),14,1)
self.assertEqual(mats[0].thickness(),0.02)
self.assertAlmostEqual(mats[1].molar_mass(),12,1)
self.assertAlmostEqual(mats[3].molar_mass(),12,1)
n2.thickness(0.5)
mats = mats + n2
self.assertEqual(mats.num(),5)
self.assertAlmostEqual(mats[0].molar_mass(),14,1)
self.assertEqual(mats[0].thickness(),0.02)
self.assertAlmostEqual(mats[4].molar_mass(),14,1)
self.assertEqual(mats[4].thickness(),0.5)
def test_material_calculation(self):
Water = catima.get_material(catima.material.Water)
p = catima.Projectile(1,1)
p(1000)
res = catima.calculate(p,Water)
res2 = catima.dedx(p,Water)
self.assertAlmostEqual(res.dEdxi,2.23,1)
self.assertAlmostEqual(res["dEdxi"],2.23,1)
self.assertAlmostEqual(res.dEdxi,res2,3)
res = catima.calculate(p(500),Water)
res2 = catima.dedx(p,Water)
self.assertAlmostEqual(res.dEdxi,2.76,1)
self.assertAlmostEqual(res.dEdxi,res2,3)
res = catima.calculate(p(9),Water)
res2 = catima.dedx(p,Water)
self.assertAlmostEqual(res.dEdxi,51.17,1)
self.assertAlmostEqual(res.dEdxi,res2,3)
res = catima.calculate(p(9),Water)
res = catima.calculate(p(9),Water)
self.assertAlmostEqual(res.dEdxi,51.17,1)
p(900000)
res = catima.calculate(p,Water)
res2 = catima.dedx_from_range(p,Water)
self.assertAlmostEqual(res.dEdxi,res2,3)
def test_config(self):
Water = catima.get_material(catima.material.Water)
Water.density(1.0)
Water.thickness(1.0)
p = catima.Projectile(1,1)
conf = catima.Config()
conf.dedx_straggling = catima.omega_type.bohr
conf2 = catima.Config()
conf2.dedx_straggling = catima.omega_type.atima
p(1000)
res = catima.calculate(p,Water,config=conf)
res2 = catima.calculate(p,Water,config=conf2)
self.assertAlmostEqual(res.dEdxi,res2.dEdxi,delta=1e-6)
self.assertNotAlmostEqual(res.sigma_E,res2.sigma_E,delta=1e-4)
self.assertNotAlmostEqual(res.sigma_r,res2.sigma_r,delta=1e-4)
def test_eout(self):
graphite = catima.get_material(6)
graphite.thickness(0.5)
p = catima.Projectile(12,6)
res = catima.calculate(p(1000),graphite)
res2 = catima.energy_out(p(1000),graphite)
self.assertAlmostEqual(res.Eout,997.077,1)
self.assertAlmostEqual(res["Eout"],997.077,1)
self.assertAlmostEqual(res.Eout,res2,3)
def test_eout_list(self):
graphite = catima.get_material(6)
graphite.thickness(0.5)
p = catima.Projectile(12,6)
energies = [100,500,1000]
res = catima.calculate(p(1000),graphite)
self.assertAlmostEqual(res.Eout,997.077,1)
res2 = catima.energy_out(p,graphite,energy=energies)
self.assertEqual(len(res2),len(energies))
self.assertAlmostEqual(res2[2], 997.077,1)
self.assertAlmostEqual(res2[0], catima.calculate(p(energies[0]),graphite).Eout ,1)
self.assertAlmostEqual(res2[1], catima.calculate(p(energies[1]),graphite).Eout ,1)
def test_dedx_from_range_list(self):
graphite = catima.get_material(6)
graphite.thickness(0.5)
p = catima.Projectile(12,6)
energies = [100,500,1000]
res2 = catima.dedx_from_range(p,graphite,energy=energies)
self.assertEqual(len(res2),len(energies))
self.assertEqual(len(res2),3)
for i,e in enumerate(energies):
r = catima.dedx_from_range(p, graphite, energy=e)
print(r)
print(res2)
self.assertAlmostEqual(res2[i], r, 0.1)
def test_layer_calculation(self):
p = catima.Projectile(12,6)
Water = catima.get_material(catima.material.Water)
Water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
graphite.density(2.0)
mat = catima.Layers()
mat.add(Water)
mat.add(graphite)
res = catima.calculate_layers(p(1000),mat)
self.assertEqual(len(res.results),2)
self.assertAlmostEqual(res.total_result.Eout,926.3,1)
self.assertAlmostEqual(res.total_result.sigma_a,0.00269,1)
self.assertAlmostEqual(res["Eout"],926.3,1)
self.assertAlmostEqual(res["sigma_a"],0.00269,4)
self.assertAlmostEqual(res["tof"],0.402,2)
self.assertAlmostEqual(res["Eloss"],884,0)
self.assertAlmostEqual(res[0]["Eout"],932.24,0)
self.assertAlmostEqual(res[1]["Eout"],926.3,0)
self.assertAlmostEqual(res[0]["sigma_a"],0.00258,4)
self.assertAlmostEqual(res[1]["sigma_a"],0.000774,4)
self.assertAlmostEqual(res[0]["range"],107.1,0)
self.assertAlmostEqual(res[1]["range"],111.3,0)
def test_energy_table(self):
table = catima.get_energy_table()
self.assertEqual(table[0],math.log(catima.logEmin))
#self.assertEqual(table[10],catima.energy_table(10))
self.assertEqual(len(table),catima.max_datapoints)
def test_storage(self):
p = catima.Projectile(12,6)
Water = catima.get_material(catima.material.Water)
Water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
data = catima.get_data(p, Water)
print(data[0])
et = catima.get_energy_table()
self.assertEqual(len(data),3)
self.assertEqual(len(data[0]),len(et))
self.assertEqual(len(data[0]),catima.max_datapoints)
res = catima.calculate(p(et[10]),Water)
self.assertAlmostEqual(res.range,data[0][10],6)
self.assertAlmostEqual(catima.projectile_range(p,Water),data[0][10],6)
#self.assertAlmostEqual(catima.domega2de(p,Water),data[1][10],6)
res = catima.calculate(p(et[100]),Water)
self.assertAlmostEqual(res.range,data[0][100],6)
self.assertAlmostEqual(catima.projectile_range(p,Water),data[0][100],6)
#self.assertAlmostEqual(catima.domega2de(p,Water),data[1][100],6)
res = catima.calculate(p(et[200]),Water)
self.assertAlmostEqual(res.range,data[0][200],6)
self.assertAlmostEqual(catima.projectile_range(p,Water),data[0][200],6)
#self.assertAlmostEqual(catima.domega2de(p,Water),data[1][200],6)
res = catima.calculate(p(et[401]),Water)
self.assertAlmostEqual(res.range,data[0][401],6)
self.assertAlmostEqual(catima.projectile_range(p,Water),data[0][401],6)
#self.assertAlmostEqual(catima.domega2de(p,Water),data[1][401],6)
def test_python_storage_access(self):
p = catima.Projectile(12,6)
Water = catima.get_material(catima.material.Water)
Water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
data = catima.get_data(p, Water)
self.assertEqual(catima.max_storage_data,100) # assuming 50, this has to be changed manually
r = catima.storage_info()
#self.assertAlmostEqual(catima.da2de(p,Water,et[100]),data[2][100],6)
#self.assertAlmostEqual(catima.da2de(p,Water,et[400]),data[2][400],6)
if __name__ == "__main__":
unittest.main()

6
tests/test_storage.cpp
View File

@ -57,18 +57,18 @@ const lest::test specification[] =
EXPECT(catima::_storage.get_index()==0); EXPECT(catima::_storage.get_index()==0);
catima::_storage.Add(p,water); catima::_storage.Add(p,water);
auto dp = catima::_storage.Get(0); auto& dp = catima::_storage.Get(0);
EXPECT(catima::_storage.get_index()==1); EXPECT(catima::_storage.get_index()==1);
EXPECT(dp.p.A==12); EXPECT(dp.p.A==12);
EXPECT(dp.m.ncomponents()==2); EXPECT(dp.m.ncomponents()==2);
catima::_storage.Add(p,water); catima::_storage.Add(p,water);
auto dp2 = catima::_storage.Get(1); auto& dp2 = catima::_storage.Get(1);
EXPECT(catima::_storage.get_index()==1); EXPECT(catima::_storage.get_index()==1);
EXPECT(dp2.p.A==0); EXPECT(dp2.p.A==0);
EXPECT(dp2.m.ncomponents()==0); EXPECT(dp2.m.ncomponents()==0);
catima::_storage.Add(p,graphite); catima::_storage.Add(p,graphite);
auto dp3 = catima::_storage.Get(1); auto& dp3 = catima::_storage.Get(1);
EXPECT(catima::_storage.get_index()==2); EXPECT(catima::_storage.get_index()==2);
EXPECT(dp3.p.A==12); EXPECT(dp3.p.A==12);
EXPECT(dp3.m.ncomponents()==1); EXPECT(dp3.m.ncomponents()==1);