gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2025-04-13 03:18:50 -04:00
Code Issues Projects Releases Wiki Activity

Compare commits

base: gwm17:master
Branches Tags
gwm17:master
gwm17:v1.6
gwm17:v1.53.1
gwm17:v1.53
gwm17:v1.5
gwm17:v1.4
..
compare: gwm17:v1.53.1
Branches Tags
gwm17:master
gwm17:v1.6
gwm17:v1.53.1
gwm17:v1.53
gwm17:v1.5
gwm17:v1.4

No commits in common. "master" and "v1.53.1" have entirely different histories.
master ... v1.53.1

48 changed files with 282 additions and 8076 deletions
Show Stats Expand all files Collapse all files

1
.gitignore vendored
View File

@ -37,4 +37,3 @@ benchmark/*
build/*
global/*
.vscode/*
gwm_test/gwm_test

89
CMakeLists.txt
View File

@ -1,28 +1,27 @@
cmake_minimum_required(VERSION 3.14)
cmake_minimum_required(VERSION 3.2.0)
project(catima)
############ options #############
option(BUILD_SHARED_LIBS "build as shared library" OFF)
option(BUILD_SHARED_LIBS "build as shared library" ON)
option(PYTHON_MODULE "compile the Catima python module(requires numpy and cython installed)" OFF)
option(TESTS "build tests" OFF)
option(EXAMPLES "build examples" OFF)
option(APPS "build catima applications" OFF)
option(APPS "build catima applications" ON)
option(GLOBAL "build with global, sources are required" OFF)
option(REACTIONS "enable/disable nuclear reaction rate" 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(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} -ffast-math")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}")
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()
@ -32,8 +31,8 @@ MESSAGE(STATUS "Build type: " ${CMAKE_BUILD_TYPE})
################################
######### compiler flags ###########
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_EXTENSIONS OFF)
set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
MESSAGE(STATUS "install prefix: " ${CMAKE_INSTALL_PREFIX})
@ -50,30 +49,39 @@ if(GSL_INTEGRATION OR GSL_INTERPOLATION)
list(APPEND EXTRA_LIBS ${GSL_LIBRARIES} )
endif()
find_package(nurex QUIET)
if(nurex_FOUND)
message(STATUS "nurex library found")
set(NUREX ON)
list(APPEND EXTRA_LIBS nurex::nurex)
endif(nurex_FOUND)
configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
"${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h")
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
"${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h"
)
configure_file("${PROJECT_SOURCE_DIR}/init.sh.in"
"${PROJECT_BINARY_DIR}/init.sh")
"${PROJECT_BINARY_DIR}/init.sh"
)
############### main build ###########################
file(GLOB SOURCES catima/*.cpp)
file(GLOB SOURCES *.cpp)
if(GLOBAL)
file(GLOB GLOBAL_SOURCES global/*.c)
LIST (APPEND SOURCES ${GLOBAL_SOURCES})
endif(GLOBAL)
file(GLOB HEADERS catima/*.h libs/*.h)
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
ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
POSITION_INDEPENDENT_CODE ON
LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
)
target_link_libraries(catima PUBLIC ${EXTRA_LIBS})
target_compile_features(catima PRIVATE cxx_std_17)
target_link_libraries(catima ${EXTRA_LIBS})
target_include_directories(catima
PUBLIC $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
@ -86,42 +94,25 @@ FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima)
# the compiler used for C++ files
MESSAGE( STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER} )
######## for python module #########
find_package(Python COMPONENTS Interpreter Development)
if(Python_FOUND)
message(STATUS "Python found: ${Python_EXECUTABLE}")
######## for python module
find_package(PythonInterp)
if(PYTHONINTERP_FOUND)
message(STATUS "Python found: ${PYTHON_EXECUTABLE}")
endif()
if(PYTHON_MODULE)
if(NOT Python_FOUND)
if(NOT PYTHONINTERP_FOUND)
MESSAGE(SEND_ERROR "Python is required to build nurex python modules")
endif(NOT Python_FOUND)
find_package(pybind11 QUIET)
if(NOT pybind11_FOUND)
message("pybind11 not found, trying to dowload")
include(FetchContent)
FetchContent_Declare(
pybind11
GIT_REPOSITORY https://github.com/pybind/pybind11.git
GIT_TAG v2.6.2
)
FetchContent_MakeAvailable(pybind11)
endif(NOT pybind11_FOUND)
check_fmt()
#set(PYBIND11_CPP_STANDARD -std=c++14)
pybind11_add_module(pycatima pymodule/pycatima.cpp)
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 fmt::fmt)
endif(PYTHON_MODULE )
target_link_libraries(pycatima PRIVATE catima)
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)
endif(PYTHON_MODULE )
########## Sub Directories ###########
if(EXAMPLES)
@ -151,12 +142,8 @@ add_subdirectory("bin")
endif(APPS)
####### install part #######
FILE(GLOB headers "catima/*.h")
FILE(GLOB headers "*.h")
include(GNUInstallDirs)
include(CMakePackageConfigHelpers)
write_basic_package_version_file(catimaConfigVersion.cmake VERSION 1.7 COMPATIBILITY AnyNewerVersion)
install (TARGETS catima
EXPORT catimaConfig
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}

8
README.md
View File

@ -1,6 +1,5 @@
[![Codacy Badge](https://api.codacy.com/project/badge/Grade/dc251db65f7a4c06ae07380544ea08fc)](https://www.codacy.com/manual/hrosiak/catima?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=hrosiak/catima&amp;utm_campaign=Badge_Grade)
[![Language grade: C/C++](https://img.shields.io/lgtm/grade/cpp/g/hrosiak/catima.svg?logo=lgtm&logoWidth=18)](https://lgtm.com/projects/g/hrosiak/catima/context:cpp)
[![Build status](https://ci.appveyor.com/api/projects/status/39gva190iyyfrtym/branch/master?svg=true)](https://ci.appveyor.com/project/hrosiak/catima-09rdl/branch/master)
[![Documentation Status](https://readthedocs.org/projects/catima/badge/?version=latest)](https://catima.readthedocs.io/en/latest/?badge=latest)
CATima
@ -33,13 +32,6 @@ This can be done sourcing the init.sh file, which is generated in the build dire
> source init.sh
```
Python Module
-------------
Python module can be installed for Linux and Windows and python versions 3.7-3.10 using pip:
```
pip install pycatima
```
cmake options
-------------
compile options, enable or disable with cmake:

0
catima/abundance_database.cpp → abundance_database.cpp
View File

0
catima/abundance_database.h → abundance_database.h
View File

42
appveyor.yml
View File

@ -1,42 +0,0 @@
version: '{build}'
platform:
- x64
environment:
matrix:
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
CONFIG: Release
PYTHON: "C:\\Python37-x64"
PYTHON_ARCH: 64
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
CONFIG: Release
PYTHON: "C:\\Python38-x64"
PYTHON_ARCH: 64
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
CONFIG: Release
PYTHON: "C:\\Python39-x64"
PYTHON_ARCH: 64
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
CONFIG: Release
PYTHON: "C:\\Python310-x64"
PYTHON_ARCH: 64
install:
- set PATH=%PYTHON%;%PYTHON%\Scripts;%PATH%
- pip --version
- pip install pybind11
- pip install wheel
build_script:
- echo %PATH%
# - git submodule add -b stable ../../pybind/pybind11 extern/pybind11
# - git submodule update --init
- mkdir build && cd build
- cmake -DBUILD_SHARED_LIBS=OFF -DAPPS=OFF -DPYTHON_WHEEL=OFF -G "Visual Studio 16 2019" -A%PLATFORM% ../
- cmake --build ./ --config "%CONFIG%"
- cmake -DBUILD_SHARED_LIBS=OFF -DAPPS=OFF -DPYTHON_WHEEL=ON -G "Visual Studio 16 2019" -A%PLATFORM% ../
# - python ../pymodule/setup.py bdist_wheel
artifacts:
- path: build\dist\*

1
build_config.in
View File

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

95
catima/calculations.cpp → calculations.cpp
View File

@ -1,4 +1,4 @@
#include <cmath>
#include <math.h>
#include <algorithm>
#include <cassert>
#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 = std::pow(p.Z,0.23)+std::pow(t.Z,0.23);
double zpowers = pow(p.Z,0.23)+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,44 +473,17 @@ 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*std::pow(z1+z2,8.0/3.0)/m2;
double factor=4.8184E-3*pow(z1+z2,8.0/3.0)/m2;
return factor*beta2/fine_structure/fine_structure;
}
double angular_scattering_power(const Projectile &p, const Target &t, double Es2){
double angular_scattering_variance(const Projectile &p, const Target &t){
if(p.T<=0)return 0.0;
double e=p.T;
double _p = p_from_T(e,p.A);
double beta = _p/((e+atomic_mass_unit)*p.A);
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 * ipow(p.Z,2)/(lr*ipow(_p*beta,2));
}
double angular_scattering_power(const Projectile &p, const Material &mat, double Es2){
if(p.T<=0)return 0.0;
double e=p.T;
double _p = p_from_T(e,p.A);
double beta = _p/((e+atomic_mass_unit)*p.A);
double X0 = radiation_length(mat);
//constexpr double Es2 = 198.81;
//constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass;
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){
if(p.T<=0)return 0.0;
double e=p.T;
double _p = p_from_T(e,p.A);
double beta = _p/((e+atomic_mass_unit)*p.A);
double pv = _p*beta;
double p1v1 = p1*beta1;
double Xs = scattering_length(mat);
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 * ipow(p.Z,2)/(Xs*ipow(_p*beta,2));
return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
}
/// radiation lengths are taken from Particle Data Group 2014
@ -542,7 +515,7 @@ double radiation_length(int z, double m){
if(z==92){return 6.00;}
double z2 = z*z;
double z_13 = 1.0/std::pow(z,1./3.);
double z_13 = 1.0/pow(z,1./3.);
double z_23 = z_13*z_13;
double a2 = fine_structure*fine_structure*z2;
double a4 = a2*a2;
@ -551,31 +524,6 @@ double radiation_length(int z, double m){
return lr;
}
double radiation_length(const Material &material){
double sum = 0;
for(int i=0;i<material.ncomponents();i++){
auto t = material.get_element(i);
double w = material.weight_fraction(i);
sum += w/radiation_length(t.Z,t.A);
}
return 1/sum;
}
double scattering_length(int a, int z){
double c = 0.00034896*z*z*(2*std::log(33219*std::pow(a*z,-1.0/3.0))-1)/a;
return 1.0/c;
}
double scattering_length(const Material& material){
double sum = 0;
for(int i=0;i<material.ncomponents();i++){
auto t = material.get_element(i);
double w = material.weight_fraction(i);
sum += w/scattering_length(t.A,t.Z);
}
return 1/sum;
}
double precalculated_lindhard(const Projectile &p){
double T = p.T;
@ -619,11 +567,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*ipow(zp_eff,2)*t.Z/t.A;
double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
if( (c.calculation == omega_types::atima) ){
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 = 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 = std::max(cor, 0.0 );
}
//double X = bethek_lindhard_X(p);
@ -665,19 +613,18 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){
return z_eff_Schiwietz(p.Z, beta, t.Z);
}
else{
assert(false);
return 0.0;
}
}
double z_eff_Pierce_Blann(double z, double beta){
return z*(1.0-exp(-0.95*fine_structure_inverted*beta/std::pow(z,2.0/3.0)));
return z*(1.0-exp(-0.95*fine_structure_inverted*beta/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/std::pow(pz,2.0/3.0))));
return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0))));
}
double z_eff_Hubert(double pz, double E, double tz){
@ -705,7 +652,7 @@ double z_eff_Hubert(double pz, double E, double tz){
x4 = 0.5218 + 0.02521*lntz;
}
return pz*(1-x1*exp(-x2*std::pow(E,x3)*std::pow(pz,-x4)));
return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4)));
}
double z_eff_Winger(double pz, double beta, double tz){
@ -733,7 +680,7 @@ double z_eff_Winger(double pz, double beta, double tz){
if(tz==2){
tz = 2.8;
}
x = beta /0.012 /std::pow(pz,0.45);
x = beta /0.012 /pow(pz,0.45);
lnz =log(pz);
lnzt=log(tz);
@ -749,17 +696,13 @@ double z_eff_Winger(double pz, double beta, double tz){
double z_eff_global(double pz, double E, double tz){
if(E>2000)
return pz;
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){
@ -767,10 +710,10 @@ double z_eff_Schiwietz(double pz, double beta, double tz){
double c1, c2;
double x;
scaled_velocity = std::pow(pz,-0.543)*beta/bohr_velocity;
scaled_velocity = catima::power(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*std::pow(scaled_velocity/c2/1.54,1+(1.83/pz));
x = c1*catima::power(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) );
}
@ -790,17 +733,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/std::pow(tz+1,c2)*(1-exp(-c3*T)) + qpb;
qglobal = (qglobal - qpb)*c1/catima::power(tz+1,c2)*(1-exp(-c3*T)) + qpb;
}
emax = 1500.;
emin = 1000.;
if(T>emax){
qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
qmean = qhigh;
}
else if(T>=emin && T<=emax){
qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
if(pz<=28){
qwinger = z_eff_Winger(pz,beta,tz);
qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin);

50
catima/calculations.h → calculations.h
View File

@ -100,25 +100,9 @@ namespace catima{
*/
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
//constexpr double Es2_FR =2*PI/fine_structure* electron_mass * electron_mass;
constexpr double Es2_FR = 198.81;
/**
* angular scattering power in form of da^2/dx in units rad^2/ g/cm^2
* @param p - Projectile class
* @param t - Target class
* @Es2 - energy constant squared, default is 14.1^2 = 198.81
*/
double angular_scattering_power(const Projectile &p, const Target &t, double Es2=Es2_FR);
/**
* angular scattering power in form of da^2/dx in units rad^2/ g/cm^2
* @param p - Projectile class
* @param t - Material class
* @Es2 - energy constant squared, default is 14.1^2 = 198.81
*/
double angular_scattering_power(const Projectile &p, const Material &material, double Es2=Es2_FR);
double angular_scattering_power_xs(const Projectile &p, const Material &mat, double p1, double beta1, double Es2=225.0);
double angular_scattering_variance(const Projectile &p, const Target &t);
/**
* returns radiation length of the (M,Z) material
* for certain z the radiation length is tabulated, otherwise calculated
@ -127,28 +111,8 @@ namespace catima{
* @return radiation length in g/cm^2
*/
double radiation_length(int z, double m);
/**
* returns radiation length of the Material class
* radiation length if calculated if not specified in Material
* or return specified radiation length
* @param Material - Material class
* @return radiation length in g/cm^2
*/
double radiation_length(const Material &material);
/**
* returns radiation length of the Material class
* radiation length if calculated if not specified in Material
* or return specified radiation length
* @param Material - Material class
* @return radiation length in g/cm^2
*/
double scattering_length(const Material &material);
/** returns effective Z of the projectile
* @param p - Projectile class
* @param t - Target class
* @param c - Configuration, the z effective will be calculated according to c.z_effective value
* @return - z effective
*/
@ -230,13 +194,5 @@ 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

208
catima/catima.cpp → catima.cpp
View File

@ -1,5 +1,5 @@
#include <iostream>
#include <cmath>
#include <math.h>
#include <algorithm>
#include "catima/catima.h"
#include "catima/constants.h"
@ -57,6 +57,18 @@ double domega2dx(const Projectile &p, const Material &mat, const Config &c){
return sum;
}
double da2dx(const Projectile &p, const Material &mat, const Config &c){
double sum = 0;
for(int i=0;i<mat.ncomponents();i++){
auto t = mat.get_element(i);
double w = mat.weight_fraction(i);
sum += w*angular_scattering_variance(p,t);
}
return sum;
}
double range(const Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
@ -108,90 +120,26 @@ double domega2de(const Projectile &p, double T, const Material &t, const Config
spline_type range_straggling_spline = get_range_straggling_spline(data);
return range_straggling_spline.derivative(T);
}
double da2dx(const Projectile &p, const Material &mat, const Config &c){
double Es2 = 198.81;
double f = 1.0;
if(c.scattering == scattering_types::dhighland)Es2 = 15*15;
if(c.scattering == scattering_types::fermi_rossi)Es2 = 15*15;
return f*angular_scattering_power(p,mat, Es2);
}
/*
double da2de(const Projectile &p, double T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
//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);
}
*/
double da2de(const Projectile &p, const Material &mat, const Config &c){
return da2dx(p,mat,c)/dedx(p,mat,c);
}
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*ipow(_p*beta,2));
}
double angular_variance(Projectile p, const Material &t, const Config &c, int order){
const double T = p.T;
const double p1 = p_from_T(T,p.A);
const double beta1 = p1/((T+atomic_mass_unit)*p.A);
assert(T>0.0);
assert(t.density()>0.0);
assert(t.thickness()>0.0);
auto& data = _storage.Get(p,t,c);
spline_type range_spline = get_range_spline(data);
double range = range_spline(T);
double rrange = std::min(range/t.density(), t.thickness_cm()); // residual range, in case of stopping inside material
double X0 = radiation_length(t);
double Es2 = 198.81;
if(c.scattering == scattering_types::dhighland)Es2 = 15*15;
if(c.scattering == scattering_types::fermi_rossi)Es2 = 15*15;
auto fx0p = [&](double x)->double{
double e =energy_out(T,x*t.density(),range_spline);
double d = ipow((rrange-x),order);
double ff = 1;
if(c.scattering == scattering_types::dhighland){
double l = x*t.density()/X0;
double lnl = log(l);
ff = 0.97*(1+lnl/20.7)*(1+lnl/22.7);
}
//return d*ff*da2dx(p(e), t, c);
return d*ff*Tfr(p(e),1, 1.0);
};
auto fx0p_2 = [&](double x)->double{
double e =energy_out(T,x*t.density(),range_spline);
double d = ipow((rrange-x),order);
return d*angular_scattering_power_xs(p(e),t,p1,beta1);
};
// corrections
if(c.scattering == scattering_types::gottschalk){
return integrator.integrate(fx0p_2,0, rrange)*t.density();
}
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){
return sqrt(angular_variance(p,t,c));
}
double angular_straggling_from_E(const Projectile &p, double Tout, Material t, const Config &c){
double angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
spline_type range_spline = get_range_spline(data);
double th = range_spline(p.T)-range_spline(Tout);
t.thickness(th);
return angular_straggling(p,t,c);
//Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
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(const 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_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);
@ -211,7 +159,7 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){
e = T - (thickness*dedx);
while(1){
r = range - range_spline(e) - thickness;
if(fabs(r)<Eout_th_epsilon)return e;
if(fabs(r)<Eout_epsilon)return e;
double step = -r*dedx;
e = e-step;
if(e<Ezero)return 0.0;
@ -272,94 +220,84 @@ Result calculate(Projectile p, const Material &t, const Config &c){
if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
auto& data = _storage.Get(p,t,c);
bool use_angular_spline = false;
if(c.scattering == scattering_types::atima_scattering){
use_angular_spline = true;
}
//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);
res.Ein = T;
res.range = range_spline(T);
res.dEdxi = p.A/range_spline.derivative(T);
res.sigma_r = sqrt(range_straggling_spline(T));
if(t.thickness()==0){
res.dEdxo = res.dEdxi;
res.Eout = res.Ein;
return res;
}
res.dEdxi = p.A/range_spline.derivative(T);
res.Eout = energy_out(T,t.thickness(),range_spline);
res.Eloss = (res.Ein - res.Eout)*p.A;
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
spline_type range_straggling_spline = get_range_straggling_spline(data);
spline_type angular_variance_spline = get_angular_variance_spline(data);
if(res.Eout<Ezero){
res.dEdxo = 0.0;
res.sigma_a = 0.0;
res.tof = 0.0;
res.sigma_E = 0.0;
}
else{
spline_type angular_variance_spline = get_angular_variance_spline(data);
}
else{
res.dEdxo = p.A/range_spline.derivative(res.Eout);
#ifdef THIN_TARGET_APPROXIMATION
if(thin_target_limit*res.Ein<res.Eout){
double edif = (res.Ein-res.Eout);
double edif = (res.Ein-res.Eout);
double s1 = range_straggling_spline.derivative(T);
double s2 = range_straggling_spline.derivative(res.Eout);
res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A;
if(use_angular_spline){
s1 = angular_variance_spline.derivative(T);
s2 = angular_variance_spline.derivative(res.Eout);
res.sigma_a = sqrt(0.5*(s1+s2)*edif);
}
s1 = angular_variance_spline.derivative(T);
s2 = angular_variance_spline.derivative(res.Eout);
res.sigma_a = sqrt(0.5*(s1+s2)*edif);
}
else{
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
if(use_angular_spline){
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
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
if(use_angular_spline){
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
}
//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));
#endif
if( (!use_angular_spline) && res.range>t.thickness()){ // do not calculate angle scattering when stopped inside material
res.sigma_a = angular_straggling(p(T),t,c);
}
//Interpolator tof_spline(energy_table.values, tofdata.data(), energy_table.num,interpolation_t::linear);
//res.tof = tof_spline(res.Ein) - tof_spline(res.Eout);
res.tof = calculate_tof_from_E(p,res.Eout,t);
} //end of else for non stopped case
if( t.thickness()>0){
//auto tofdata = calculate_tof(p,t,c);
//Interpolator tof_spline(energy_table.values, tofdata.data(), energy_table.num,interpolation_t::linear);
//res.tof = tof_spline(res.Ein) - tof_spline(res.Eout);
res.tof = calculate_tof_from_E(p,res.Eout,t);
}
}
res.sigma_r = sqrt(range_straggling_spline(T));
res.Eloss = (res.Ein - res.Eout)*p.A;
// position straggling in material
double rrange = std::min(res.range/t.density(), t.thickness_cm());
res.sigma_x = angular_variance(p(T),t,c,2);
double rrange = std::min(res.range/t.density(), t.thickness_cm());
auto fx2p = [&](double x)->double{
double range = range_spline(T);
double e =energy_out(T,x*t.density(),range_spline);
return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
};
//res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3, 1e-3,1);
res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange);
res.sigma_x = sqrt(res.sigma_x);
rrange = std::min(res.range/t.density(), t.thickness_cm());
// position vs angle covariance, needed later for final position straggling
res.cov = angular_variance(p(T),t,c,1);
auto fx1p = [&](double x)->double{
double e =energy_out(T,x*t.density(),range_spline);
return (rrange-x)*da2dx(p(e), t, c)*t.density();
};
//res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,1);
res.cov = integrator.integrate(fx1p,0, rrange);
#ifdef REACTIONS
res.sp = nonreaction_rate(p,t,c);
#endif
return res;
}
MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &layers, const Config &c){
MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c){
MultiResult res;
double e = p.T;
res.total_result.Ein = e;
res.total_result.sigma_a = ps.sigma_a*ps.sigma_a;
res.total_result.sigma_x = ps.sigma_x*ps.sigma_x;
res.total_result.cov = ps.cov_x;
res.results.reserve(layers.num());
for(auto&m:layers.get_materials()){
Result r = calculate(p,m,e,c);
@ -372,6 +310,7 @@ MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &l
res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov)
+ (a2*m.thickness_cm()*m.thickness_cm())
+ r.sigma_x*r.sigma_x;
//res.total_result.sigma_x += (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x;
res.total_result.cov += a2*m.thickness_cm() + r.cov;
res.total_result.sigma_a += r.sigma_a*r.sigma_a;
#ifdef REACTIONS
@ -382,13 +321,13 @@ MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &l
if(e>Ezero){
res.total_result.sigma_a = sqrt(res.total_result.sigma_a);
res.total_result.sigma_E = sqrt(res.total_result.sigma_E);
res.total_result.sigma_x = sqrt(std::abs(res.total_result.sigma_x));
res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
}
else{
res.total_result.sigma_a = 0.0;
res.total_result.sigma_E = 0.0;
res.total_result.sigma_x = sqrt(std::abs(res.total_result.sigma_x));
res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
}
return res;
}
@ -410,9 +349,6 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
auto fomega = [&](double x)->double{
return domega2dx(p(x),t,c)/catima::power(dedx(p(x),t,c),3);
};
auto ftheta = [&](double x)->double{
return da2de(p(x),t,c);
};
//double res=0.0;
//calculate 1st point to have i-1 element ready for loop
@ -426,14 +362,12 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
//res = integrator.integrate(fomega,Ezero,energy_table(0));
//res = p.A*res;
dp.range_straggling[0]=0.0;
//p.T = energy_table(0);
//p.T = energy_table(0);
for(int i=1;i<max_datapoints;i++){
double 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];
dp.angular_variance[i] = p.A*integrator.integrate(ftheta,energy_table(i-1),energy_table(i))
+ dp.angular_variance[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;

32
catima/catima.h → catima.h
View File

@ -113,23 +113,6 @@ namespace catima{
*/
double da2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
/**
* returns the planar RMS angular straggling in rad
* @param p - Projectile
* @param t - Material class
* @param c - Config class
* @return angular RMS straggling in rad
*/
double angular_straggling(Projectile p, const Material &t, const Config &c=default_config);
/**
* returns the planar RMS angular variance in rad
* @param p - Projectile
* @param t - Material class
* @param c - Config class
* @return angular RMS variance in rad
*/
double angular_variance(Projectile p, const Material &t, const Config &c=default_config, int order = 0);
/**
* calculates angular scattering in the material from difference of incoming a nd outgoing energies
* @param p - Projectile
@ -138,7 +121,7 @@ namespace catima{
* @param mat - Material
* @return angular straggling
*/
double angular_straggling_from_E(const Projectile &p, double Tout,Material t, const Config &c=default_config);
double angular_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
/**
* calculates Energy straggling in the material from difference of incoming a nd outgoing energies
@ -149,7 +132,7 @@ namespace catima{
* @return angular straggling
*/
double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
/**
* calculates outcoming energy from range spline
* @param T - incoming energy
@ -203,16 +186,7 @@ namespace catima{
* @return results stored in MultiResult structure
*
*/
MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &layers, const Config &c=default_config);
/**
* calculate observables for multiple layers of material defined by Layers
* @return results stored in MultiResult structure
*
*/
inline MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config){
return calculate(p, {}, layers, c);
};
MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
return calculate(p(T), layers, c);
}

81
catima/convert.h
View File

@ -1,81 +0,0 @@
/*
* Author: Andrej Prochazka
* Copyright(C) 2017
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Modifed by GWM to remove unecessary fmt dependance (only place used in entire library)
*/
#ifndef CONVERT_H
#define CONVERT_H
#include "catima/structures.h"
#include "catima/material_database.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
using namespace catima;
bool mocadi_material_match(const Material &a, const Material&b){
if(std::fabs(a.density() - b.density())> 1e-6)return false;
if(a.ncomponents() != b.ncomponents())return false;
for(int i=0;i<a.ncomponents();i++){
if(a.get_element(i).stn != b.get_element(i).stn)return false;
if(a.get_element(i).A != b.get_element(i).A)return false;
if(a.get_element(i).Z != b.get_element(i).Z)return false;
}
if(a.M() != b.M())return false;
return true;
}
int materialdb_id(const Material &m){
for(int i=201;i<=347;i++){
if(mocadi_material_match(get_material(i),m)){
return i;
}
}
return 0;
}
bool save_mocadi(const char* filename, const Projectile p, const Layers &layers, const Phasespace &psx={}, const Phasespace &psy={}){
std::ofstream fw;
fw.open(filename, std::ios::out);
if (!fw.is_open()) { return false;}
fw<<"epax 2\natima-1.0\noption listmode root\n";
std::string beam = "BEAM\n100000\n" + std::to_string(p.T) + ", 0, " + std::to_string(p.A) + ", " + std::to_string(p.Z) + "\n2\n" +
std::to_string(psx.sigma_x) + ", " + std::to_string(1000*psx.sigma_a) + ", 0, 0, 0\n2\n" +
std::to_string(psy.sigma_x) + ", " + std::to_string(1000*psy.sigma_a) + ", 0, 0, 0\n1\n0, 0, 0, 0, 0";
fw<<beam;
int c = 0;
for (auto& m: layers.get_materials()){
int z = 0;
double a = 0.0;
if(m.ncomponents()==1){
z = m.get_element(0).Z;
a = m.get_element(0).A;
}
else{
z = materialdb_id(m);
}
std::string mstr = "*******\nMATTER\n" + std::to_string(a) + ", " + std::to_string(z) + ", " + std::to_string(m.density()*1000.0) +
"\n" + "2,"+std::to_string(m.thickness_cm()) + "\n0.\n0.,0.,0.\n1,1,0\n0,0,0,0\n";
fw<<mstr;
c++;
}
fw<<"ERWARTUNGSWERTE\nSAVE\nEND";
fw.close();
return true;
}
#endif

96
catima/cwrapper.cpp
View File

@ -1,96 +0,0 @@
#include "catima/cwrapper.h"
#include "catima/catima.h"
#include "catima/material_database.h"
#include "catima/nucdata.h"
#include "catima/reactions.h"
#include <cstring>
extern "C" {
struct CatimaConfig catima_defaults = {1};
catima::Material make_material(double ta, double tz, double thickness, double density){
catima::Material mat;
if(tz>200){
mat = catima::get_material(tz);
}
else{
mat.add_element(ta,tz,1.0);
}
mat.density(density).thickness(thickness);
if(density<=0.0){
catima::Material m0 = catima::get_material(tz);
mat.density(m0.density());
}
return mat;
}
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);
CatimaResult res;
res.Ein = r.Ein;
res.Eout = r.Eout;
res.Eloss = r.Eloss;
res.range = r.range;
res.dEdxi = r.dEdxi;
res.dEdxo = r.dEdxo;
res.sigma_E = r.sigma_E;
res.sigma_a = r.sigma_a;
res.sigma_r = r.sigma_r;
res.tof = r.tof;
// printf("%d\n",catima::_storage.get_index());
return res;
}
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);
}
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);
}
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);
return catima::angular_straggling_from_E(p(Tin),Tout,mat);
}
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);
}
double atomic_weight(int i){return catima::element_atomic_weight(i);}
double catima_nonreaction_rate(double pa, int pz, double T, double ta, double tz, double thickness){
catima::Projectile p(pa,pz);
p.T = T;
catima::Material mat = make_material(ta,tz, thickness, -1);
return catima::nonreaction_rate(p,mat);
}
}

97
catima/gwm_integrators.cpp
View File

@ -1,97 +0,0 @@
/*
* Simple energy loss integrators for forward and reverse cases.
*
* Dec 2022, Gordon McCann
*/
#include "catima/gwm_integrators.h"
namespace catima {
double integrate_energyloss(Projectile& proj, const Material& mat, const Config& c)
{
static double s_estep_max = 0.001;
static int s_depth_max = 100;
int depth = 0;
double e_in = proj.T; // MeV/u
double e_final = e_in;
double x_step = 0.25*mat.thickness(); //g/cm^2
double x_traversed = 0.0;
double e_step = dedx(proj, mat, c)*x_step;
double A_recip = 1.0/proj.A;
while(true)
{
if(e_step/e_final > s_estep_max && depth < s_depth_max)
{
++depth;
x_step *= 0.5;
e_step = dedx(proj, mat, c)*x_step*A_recip;
}
else if(x_step + x_traversed >= mat.thickness())
{
x_step = mat.thickness() - x_traversed;
e_step = dedx(proj, mat, c)*x_step*A_recip;
e_final -= e_step;
proj.T = e_final;
return (e_in - e_final)*proj.A;
}
else if(depth == s_depth_max)
{
return e_in*proj.A;
}
else if (e_final < 0.0)
{
return e_in * proj.A; //In case an integration step takes us below 0
}
else
{
e_step = dedx(proj, mat, c)*x_step*A_recip;
e_final -= e_step;
proj.T = e_final;
x_traversed += x_step;
}
}
}
double reverse_integrate_energyloss(Projectile& proj, const Material& mat, const Config& c)
{
static double s_estep_max = 0.001;
static int s_depth_max = 100;
int depth = 0;
double e_out = proj.T; //MeV/u
double e_initial = e_out;
double x_step = 0.25*mat.thickness(); //g/cm^2
double x_traversed = 0.0;
double e_step = dedx(proj, mat, c)*x_step;
double A_recip = 1.0/proj.A;
while(true)
{
if(e_step/e_initial > s_estep_max && depth < s_depth_max)
{
++depth;
x_step *= 0.5;
e_step = dedx(proj, mat, c)*x_step*A_recip;
}
else if(x_step + x_traversed >= mat.thickness())
{
x_step = mat.thickness() - x_traversed;
e_step = dedx(proj, mat, c)*x_step;
e_initial += e_step*A_recip;
proj.T = e_initial;
return (e_initial - e_out)*proj.A;
}
else if(depth == s_depth_max)
{
return e_out*proj.A;
}
else
{
e_step = dedx(proj, mat, c)*x_step;
e_initial += e_step*A_recip;
proj.T = e_initial;
x_traversed += x_step;
}
}
}
}

18
catima/gwm_integrators.h
View File

@ -1,18 +0,0 @@
/*
* Simple energy loss integrators for forward and reverse cases.
*
* Dec 2022, Gordon McCann
*/
#ifndef GWM_INTEGRATORS_H
#define GWM_INTEGRATORS_H
#include "catima/catima.h"
namespace catima {
double integrate_energyloss(Projectile& proj, const Material& mat, const Config& c=default_config);
double reverse_integrate_energyloss(Projectile& proj, const Material& mat, const Config& c=default_config);
}
#endif

13
catima/config.h → config.h
View File

@ -45,16 +45,6 @@ namespace catima{
srim_95 = 1,
};
/**
* enum to select angular scattering
*/
enum scattering_types:unsigned char{
fermi_rossi = 0,
dhighland = 1,
gottschalk = 2,
atima_scattering = 255,
};
/**
* structure to store calculation configuration
*/
@ -67,8 +57,7 @@ namespace catima{
unsigned char corrections = 0;
unsigned char calculation = 1;
unsigned char low_energy = low_energy_types::srim_85;
unsigned char scattering = scattering_types::atima_scattering;
unsigned char low_energy = 0;
};

2
catima/constants.h → constants.h
View File

@ -11,7 +11,7 @@ constexpr double logEmax = 7.0; // log of max energy
constexpr int max_datapoints = 600; // how many datapoints between logEmin and logEmax
constexpr int max_storage_data = 60; // number of datapoints which can be stored in cache
constexpr double numeric_epsilon = 10*std::numeric_limits<double>::epsilon();
constexpr double Eout_th_epsilon = 1e-5; //
constexpr double Eout_epsilon = 1e-5; //
constexpr double thin_target_limit = 1 - 1e-3;

51
cwrapper.cpp Normal file
View File

@ -0,0 +1,51 @@
#include "catima/cwrapper.h"
#include "catima/catima.h"
#include "catima/material_database.h"
#include <cstring>
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::Material mat;
catima::Projectile p(pa,pz);
if(tz>200){
mat = catima::get_material(tz);
}
else{
mat.add_element(ta,tz,1.0);
}
mat.density(density).thickness(thickness);
catima::Result r = catima::calculate(p(T),mat);
CatimaResult res;
std::memcpy(&res,&r,sizeof(res));
return res;
}
double catima_angular_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){
catima::Projectile p(pa,pz);
catima::Material mat;
if(tz>200){
mat = catima::get_material(tz);
}
else{
mat.add_element(ta,tz,1.0);
}
return catima::angular_straggling_from_E(p,Tin,Tout,mat);
}
double catima_energy_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){
catima::Projectile p(pa,pz);
catima::Material mat;
if(tz>200){
mat = catima::get_material(tz);
}
else{
mat.add_element(ta,tz,1.0);
}
return catima::energy_straggling_from_E(p,Tin,Tout,mat);
}
}

8
catima/cwrapper.h → cwrapper.h
View File

@ -34,18 +34,14 @@ struct CatimaConfig {
char z_effective;
};
extern struct CatimaConfig catima_defaults;
struct CatimaConfig catima_defaults = {none};
typedef struct CatimaResult CatimaResult;
CatimaResult catima_calculate(double pa, int pz, double T, double ta, double tz, double thickness, double density);
double catima_Eout(double pa, int pz, double T, double ta, double tz, double thickness, double density);
double catima_range(double pa, int pz, double T, double ta, double tz);
double catima_range_straggling(double pa, int pz, double T, double ta, double tz);
double catima_angular_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz);
double catima_energy_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz);
double atomic_weight(int i);
double catima_nonreaction_rate(double pa, int pz, double T, double ta, double tz, double thickness);
#ifdef __cplusplus
}

0
catima/data_atima.h → data_atima.h
View File

0
catima/data_ionisation_potential.h → data_ionisation_potential.h
View File

6
docs/README.md
View File

@ -18,12 +18,6 @@ This can be done sourcing the init.sh file, which is generated in the build dire
source init.sh
```
Python Module
-------------
Python module can be installed also using pip:
```
pip install pycatima
```
cmake options
-------------

10
docs/pycatima.md
View File

@ -1,13 +1,5 @@
Usage
=====
Installation
------------
Easiest way is to install pycatima on Linux and Windows is using pip:
```
pip install pycatima
```
note: python 3.7-3.9 is required
Pojectile
---------
@ -166,7 +158,7 @@ __Results__ class stores results for 1 layer of __Material__. It has following v
* dEdxo - Stopping power at exit
* range - range in the material in g/cm^2
* sigma_E - Energy straggling in MeV/u
* sigma_a - Angular straggling in rad
* sigma_a - Angular straggling in MeV/u
* sigma_r - range straggling in g/cm^2
* sigma_x - position straggling in cm
* sp - non-reaction probability

0
catima/generated_LS_coeff.h → generated_LS_coeff.h
View File

18
gwm_test/Makefile
View File

@ -1,18 +0,0 @@
CC=g++
EXE=gwm_test
CXX_FLAGS= -std=c++17 -g -Wall
CATIMA_PATH=../build/
LIB_PATH=$(CATIMA_PATH)lib/
INCLUDE_PATH=$(CATIMA_PATH)include/
CPP=gwm_test.cpp
.PHONY: all clean
all: $(EXE)
$(EXE): $(CPP)
$(CC) $(CXX_FLAGS) -I$(INCLUDE_PATH) $^ $(LIB_PATH)libcatima.a -o $@
clean:
$(RM) $(EXE)

22
gwm_test/gwm_test.cpp
View File

@ -1,22 +0,0 @@
/*
* Testing for gwm_integrators
*
* Dec 2022, Gordon McCann
*/
#include "catima/gwm_integrators.h"
#include "catima/nucdata.h"
#include <iostream>
int main()
{
std::cout<<"-------Testing GWM Energy Loss Integration-------"<<std::endl;
catima::Projectile p1(catima::element_atomic_weight(1), 1.0, 0, 3.0);
catima::Material mat1(catima::get_material(6));
mat1.density(2.23).thickness(500.0*1e-6);
double result = catima::integrate_energyloss(p1, mat1);
std::cout<<"Energy loss (MeV): "<<result<<" Final energy: "<<p1.T<<std::endl;
result = catima::reverse_integrate_energyloss(p1, mat1);
std::cout<<"Reverse Energy loss (MeV): "<<result<<" Initial energy: "<<p1.T<<std::endl;
}

0
catima/integrator.cpp → integrator.cpp
View File

2
catima/integrator.h → integrator.h
View File

@ -62,7 +62,7 @@ using integrator_type = IntegratorGSL;
#else
using integrator_type = GaussLegendreIntegration<8>;
#endif
using integrator_adaptive_type = GaussKronrodIntegration<21>;
using integrator_adaptive_type = GaussKronrodIntegration<15>;
extern integrator_type integrator;
extern integrator_adaptive_type integrator_adaptive;

1
catima/material_database.cpp → material_database.cpp
View File

@ -158,7 +158,6 @@ namespace catima{
case material::C21H24O4: return Material({{0,6,21},{0,1,24},{0,8,4}},1.18);
case material::CoRe_Alloy: return Material({{0,27,17},{0,75,23},{0,24,1}},11.5);
case material::LLZO_electrolyte: return Material({{0,3,7},{0,57,3},{0,40,2},{0,8,12}},5.1);
case material::Nylon: return Material({{0,6,6},{0,1,11},{0,7,1},{0,8,1}},1.14);
default:break;
}
return Material();

3
catima/material_database.h → material_database.h
View File

@ -150,8 +150,7 @@ namespace catima{
Iodonaphthalene = 343,
C21H24O4 = 344,
CoRe_Alloy = 345,
LLZO_electrolyte = 346,
Nylon = 347
LLZO_electrolyte = 346
};
Material get_compound(material m);

0
catima/nucdata.cpp → nucdata.cpp
View File

0
catima/nucdata.h → nucdata.h
View File

228
pymodule/pycatima.cpp
View File

@ -6,15 +6,16 @@
#include "catima/catima.h"
#include "catima/srim.h"
#include "catima/nucdata.h"
#include "catima/convert.h"
#include <iostream>
#include <string>
namespace py = pybind11;
using namespace catima;
std::string catima_info(){
return "CATIMA version = 1.7\n";
void catima_info(){
printf("CATIMA version = 1.5\n");
printf("number of energy points = %d\n",max_datapoints);
printf("min energy point = 10^%lf MeV/u\n",logEmin);
printf("max energy point = 10^%lf MeV/u\n",logEmax);
}
std::string material_to_string(const Material &r){
@ -49,14 +50,9 @@ py::list storage_info(){
py::list get_energy_table(){
py::list r;
for (size_t i = 0; i < energy_table.size(); i++)
{
r.append(energy_table[i]);
for(auto e : energy_table){
r.append(e);
}
//for(auto e : energy_table){
//r.append(e);
//}
return r;
}
@ -105,7 +101,6 @@ py::dict get_result_dict(const Result& r){
d["sigma_r"] = r.sigma_r;
d["sigma_a"] = r.sigma_a;
d["sigma_x"] = r.sigma_x;
d["cov"] = r.cov;
d["tof"] = r.tof;
d["sp"] = r.sp;
return d;
@ -132,12 +127,6 @@ PYBIND11_MODULE(pycatima,m){
.def_readwrite("Z",&Target::Z)
.def_readwrite("stn",&Target::stn);
py::class_<Phasespace>(m, "Phasespace")
.def(py::init<>(),"constructor")
.def_readwrite("sigma_x", &Phasespace::sigma_x)
.def_readwrite("sigma_a", &Phasespace::sigma_a)
.def_readwrite("cov_x", &Phasespace::cov_x);
py::class_<Material>(m,"Material")
.def(py::init<>(),"constructor")
@ -187,13 +176,9 @@ PYBIND11_MODULE(pycatima,m){
.def_readwrite("sigma_a", &Result::sigma_a)
.def_readwrite("sigma_r", &Result::sigma_r)
.def_readwrite("sigma_x", &Result::sigma_x)
.def_readwrite("cov", &Result::cov)
.def_readwrite("tof", &Result::tof)
.def_readwrite("sp", &Result::sp)
.def("get_dict",&get_result_dict)
.def("__repr__",[](const Result &self){
return py::str(get_result_dict(self));
});
.def("get_dict",&get_result_dict);
py::class_<MultiResult>(m,"MultiResult")
.def(py::init<>(),"constructor")
@ -229,17 +214,7 @@ PYBIND11_MODULE(pycatima,m){
}
d["partial"] = p;
return d;
})
.def("__repr__",[](const MultiResult &r){
py::dict d;
py::list p;
d["total_result"] = get_result_dict(r.total_result);
for(auto& entry:r.results){
p.append(get_result_dict(entry));
}
d["results"] = p;
return py::str(d);
});
});
py::enum_<z_eff_type>(m,"z_eff_type")
.value("none", z_eff_type::none)
@ -264,163 +239,30 @@ PYBIND11_MODULE(pycatima,m){
py::enum_<low_energy_types>(m,"low_energy_types")
.value("srim_85", low_energy_types::srim_85)
.value("srim_95", low_energy_types::srim_95);
py::enum_<scattering_types>(m,"scattering_types")
.value("fermi_rossi", scattering_types::fermi_rossi)
.value("dhighland", scattering_types::dhighland)
.value("gottschalk", scattering_types::gottschalk)
.value("atima_scattering", scattering_types::atima_scattering);
py::enum_<material>(m,"material")
.value("Plastics", material::Plastics)
.value("Air", material::Air)
.value("CH2", material::CH2)
.value("lH2", material::lH2)
.value("lD2", material::lD2)
.value("Water", material::Water)
.value("Diamond", material::Diamond)
.value("Glass", material::Glass)
.value("ALMG3", material::ALMG3)
.value("ArCO2_30", material::ArCO2_30)
.value("CF4", material::CF4)
.value("Isobutane", material::Isobutane)
.value("Kapton", material::Kapton)
.value("Mylar", material::Mylar)
.value("NaF", material::NaF)
.value("P10", material::P10)
.value("Polyolefin", material::Polyolefin)
.value("CmO2", material::CmO2)
.value("Suprasil", material::Suprasil)
.value("HAVAR", material::HAVAR)
.value("Steel", material::Steel)
.value("CO2", material::CO2)
.value("Methane", material::Methane)
.value("Methanol", material::Methanol)
.value("Acetone", material::Acetone)
.value("Acetylene", material::Acetylene)
.value("Adenine", material::Adenine)
.value("Adipose_Tissue", material::Adipose_Tissue)
.value("Alanine", material::Alanine)
.value("Bakelite", material::Bakelite)
.value("AgBr", material::AgBr)
.value("AgCl", material::AgCl)
.value("AgI", material::AgI)
.value("Al2O3", material::Al2O3)
.value("Amber", material::Amber)
.value("Ammonia", material::Ammonia)
.value("Aniline", material::Aniline)
.value("Anthracene", material::Anthracene)
.value("A_150", material::A_150)
.value("B_100", material::B_100)
.value("BaF2", material::BaF2)
.value("BaSO4", material::BaSO4)
.value("Benzene", material::Benzene)
.value("BeO", material::BeO)
.value("BGO", material::BGO)
.value("Blood_ICRP", material::Blood_ICRP)
.value("Bone_Compact", material::Bone_Compact)
.value("Bone_Cortical", material::Bone_Cortical)
.value("Brain_ICRP", material::Brain_ICRP)
.value("B4C", material::B4C)
.value("BC_400", material::BC_400)
.value("nButanol", material::nButanol)
.value("C_552", material::C_552)
.value("CdTe", material::CdTe)
.value("CdWO4", material::CdWO4)
.value("CaCO3", material::CaCO3)
.value("CaF2", material::CaF2)
.value("CaO", material::CaO)
.value("CaWO4", material::CaWO4)
.value("CsF", material::CsF)
.value("CsI", material::CsI)
.value("CCl4", material::CCl4)
.value("Tetrachloroethylene", material::Tetrachloroethylene)
.value("Cellophane", material::Cellophane)
.value("Chlorobenzene", material::Chlorobenzene)
.value("Chloroform", material::Chloroform)
.value("Cyclohexane", material::Cyclohexane)
.value("Concrete", material::Concrete)
.value("Diethyl_Ether", material::Diethyl_Ether)
.value("Ethane", material::Ethane)
.value("Ethanol", material::Ethanol)
.value("Ethylene", material::Ethylene)
.value("Eye_lens", material::Eye_lens)
.value("Fe2O3", material::Fe2O3)
.value("FeO", material::FeO)
.value("Freon_12", material::Freon_12)
.value("Freon_12B2", material::Freon_12B2)
.value("Freon_13", material::Freon_13)
.value("Freon_13B1", material::Freon_13B1)
.value("Freon_13I1", material::Freon_13I1)
.value("Gd2O2S", material::Gd2O2S)
.value("GaAs", material::GaAs)
.value("Gel_Photo_Emulsion", material::Gel_Photo_Emulsion)
.value("Glass_Pyrex", material::Glass_Pyrex)
.value("Glass_Lead", material::Glass_Lead)
.value("Glucose", material::Glucose)
.value("Glutamine", material::Glutamine)
.value("Glycerol", material::Glycerol)
.value("Guanine", material::Guanine)
.value("Gypsum", material::Gypsum)
.value("nHeptane", material::nHeptane)
.value("nHexane", material::nHexane)
.value("KI", material::KI)
.value("K2O", material::K2O)
.value("LaBr3", material::LaBr3)
.value("LaOBr", material::LaOBr)
.value("La2O2S", material::La2O2S)
.value("Lung", material::Lung)
.value("MgCO3", material::MgCO3)
.value("MgF2", material::MgF2)
.value("MgO", material::MgO)
.value("MS20_Tissue", material::MS20_Tissue)
.value("Muscle_skeletal", material::Muscle_skeletal)
.value("Muscle_strained", material::Muscle_strained)
.value("Muscle_sucrose", material::Muscle_sucrose)
.value("Muscle_no_sucrose", material::Muscle_no_sucrose)
.value("Na2CO3", material::Na2CO3)
.value("NaI", material::NaI)
.value("NaCl", material::NaCl)
.value("Na2O", material::Na2O)
.value("NaNO3", material::NaNO3)
.value("Naphthalene", material::Naphthalene)
.value("Nitrobenzene", material::Nitrobenzene)
.value("N2O", material::N2O)
.value("Octane", material::Octane)
.value("Paraffin", material::Paraffin)
.value("nPentane", material::nPentane)
.value("PhotoEmulsion", material::PhotoEmulsion)
.value("PuO2", material::PuO2)
.value("Polyacrylonitrile", material::Polyacrylonitrile)
.value("Polycarbonate", material::Polycarbonate)
.value("PMMA", material::PMMA)
.value("POM", material::POM)
.value("Polypropylene", material::Polypropylene)
.value("Polystyrene", material::Polystyrene)
.value("Propane", material::Propane)
.value("nPropanol", material::nPropanol)
.value("PVC", material::PVC)
.value("Pyridine", material::Pyridine)
.value("SiO2", material::SiO2)
.value("Skin", material::Skin)
.value("Sucrose", material::Sucrose)
.value("Teflon", material::Teflon)
.value("TlCl", material::TlCl)
.value("Toluene", material::Toluene)
.value("Trichloroethylene", material::Trichloroethylene)
.value("WF6", material::WF6)
.value("UC2", material::UC2)
.value("UC", material::UC)
.value("UO2", material::UO2)
.value("Urea", material::Urea)
.value("Valine", material::Valine)
.value("Iodonaphthalene", material::Iodonaphthalene)
.value("C21H24O4", material::C21H24O4)
.value("CoRe_Alloy", material::CoRe_Alloy)
.value("LLZO_electrolyte", material::LLZO_electrolyte)
.value("Nylon", material::Nylon);
.value("Plastics", material::Plastics)
.value("Air", material::Air)
.value("CH2", material::CH2)
.value("lH2", material::lH2)
.value("lD2", material::lD2)
.value("Water", material::Water)
.value("Diamond", material::Diamond)
.value("Glass", material::Glass)
.value("ALMG3", material::ALMG3)
.value("ArCO2_30", material::ArCO2_30)
.value("CF4", material::CF4)
.value("Isobutane", material::Isobutane)
.value("Kapton", material::Kapton)
.value("Mylar", material::Mylar)
.value("NaF", material::NaF)
.value("P10", material::P10)
.value("Polyolefin", material::Polyolefin)
.value("CmO2", material::CmO2)
.value("Suprasil", material::Suprasil)
.value("HAVAR", material::HAVAR)
.value("Steel", material::Steel)
.value("CO2", material::CO2);
py::class_<Config>(m,"Config")
@ -429,14 +271,12 @@ PYBIND11_MODULE(pycatima,m){
.def_readwrite("corrections", &Config::corrections)
.def_readwrite("calculation", &Config::calculation)
.def_readwrite("low_energy", &Config::low_energy)
.def_readwrite("scattering", &Config::scattering)
.def("get",[](const Config &r){
py::dict d;
d["z_effective"] = r.z_effective;
d["corrections"] = r.corrections;
d["calculation"] = r.calculation;
d["low_energy"] = r.low_energy;
d["scattering"] = r.scattering;
return d;
})
.def("__str__",[](const Config &r){
@ -445,17 +285,12 @@ PYBIND11_MODULE(pycatima,m){
s += ", corrections = "+std::to_string(r.corrections);
s += ", calculation = "+std::to_string(r.calculation);
s += ", low_energy = "+std::to_string(r.low_energy);
s += ", scattering = "+std::to_string(r.scattering);
return s;
});
m.def("angular_scattering_power",py::overload_cast<const Projectile&, const Material&, double>(&angular_scattering_power),"angular scattering power in rad^2/g/cm^2",py::arg("projectile"),py::arg("material"),py::arg("Es2")=Es2_FR);
m.def("radiation_length",py::overload_cast<const Material&>(radiation_length));
m.def("srim_dedx_e",&srim_dedx_e);
m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e", py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("calculate",py::overload_cast<Projectile, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("calculate",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("layers"), py::arg("config")=default_config);
m.def("calculate",py::overload_cast<const Projectile&, const Phasespace&, const Layers&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("phasespace"),py::arg("layers"), py::arg("config")=default_config);
m.def("calculate_layers",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("dedx_from_range",py::overload_cast<const Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
m.def("dedx_from_range",py::overload_cast<const Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
@ -474,7 +309,6 @@ PYBIND11_MODULE(pycatima,m){
l.append(r.second);
return l;
});
m.def("save_mocadi", &save_mocadi,py::arg("filename"),py::arg("projectile"),py::arg("layers"),py::arg("psx")=Phasespace(), py::arg("psy")=Phasespace());
m.def("catima_info",&catima_info);
m.def("storage_info",&storage_info);
m.def("get_energy_table",&get_energy_table);

26
pymodule/setup.py
View File

@ -1,26 +0,0 @@
from pathlib import Path
from pybind11.setup_helpers import Pybind11Extension, build_ext
from setuptools import setup
DIR = Path(__file__).parents[0]
SRC = [str((DIR/'pycatima.cpp').resolve())]
example_module = Pybind11Extension(
'pycatima',
SRC,
include_dirs=['../build/include','../global'],
library_dirs=['../build/lib','../build','../build/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},
)

22
pymodule/setup.py.in
View File

@ -1,22 +0,0 @@
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},
)

0
catima/reactions.cpp → reactions.cpp
View File

0
catima/reactions.h → reactions.h
View File

10
spline.cpp
View File

@ -1,7 +1,15 @@
/*
* 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

8
catima/spline.h → spline.h
View File

@ -111,6 +111,7 @@ struct cspline_special{
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;
@ -119,9 +120,7 @@ struct cspline_special{
double operator()(double x)const{return evaluate(x);}
double evaluate(double x) const
{
const T& m_x = *table;
int idx=std::max( table->index(x), 0);
int idx=std::max( table->index(x), 0);
double h=x-m_x[idx];
double interpol;
if(x<m_x[0]) {
@ -140,7 +139,6 @@ struct cspline_special{
double deriv(double x) const
{
const T& m_x = *table;
int idx=std::max( table->index(x), 0);
double h=x-m_x[idx];
@ -164,7 +162,7 @@ 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())
):table(&x),m_y(y.data()),m_x(x.values)
{
static_assert (N>2, "N must be > 2");
tridiagonal_matrix<N> A{};

0
catima/srim.cpp → srim.cpp
View File

0
catima/srim.h → srim.h
View File

4
catima/storage.cpp → storage.cpp
View File

@ -19,11 +19,7 @@
#include "catima/catima.h"
namespace catima {
Data _storage;
#ifdef VETABLE
LogVArray<max_datapoints> energy_table(logEmin,logEmax);
#else
EnergyTable<max_datapoints> 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)){

108
catima/storage.h → storage.h
View File

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

0
catima/structures.cpp → structures.cpp
View File

10
catima/structures.h → structures.h
View File

@ -215,7 +215,7 @@ namespace catima{
double number_density_cm2(int i)const{
if(i>=atoms.size())return 0.0;
return number_density_cm2()*molar_fraction(i);
}
}
friend bool operator==(const Material &a, const Material&b);
};
@ -236,19 +236,13 @@ namespace catima{
double sigma_a=0.0;
double sigma_r=0.0;
double sigma_x=0.0;
double cov = 0.0;
double cov = 0.0;
double tof=0.0;
#ifdef REACTIONS
double sp = 1.0;
#endif
};
struct Phasespace{
double sigma_x=0.0;
double sigma_a=0.0;
double cov_x=0.0;
};
/**
* structure to store results for calculation for multiple layers of materials, ie in catima::Layers
*/

6816
tests/doctest.h
View File

File diff suppressed because it is too large Load Diff

122
tests/test_calculations.cpp
View File

@ -183,104 +183,12 @@ using namespace std;
auto res = catima::calculate(p,water);
dif = res.tof - 0.038;
CHECK( fabs(dif) < 0.01);
}
TEST_CASE("scattering length"){
catima::Material be = catima::get_material(4);
double x0 = catima::radiation_length(be);
double xs = catima::scattering_length(be);
CHECK(xs/x0 == approx(1.4,0.1));
catima::Material pb = catima::get_material(82);
x0 = catima::radiation_length(pb);
xs = catima::scattering_length(pb);
CHECK(xs/x0 == approx(1.04,0.04));
}
TEST_CASE("scattering_kanematsu"){
catima::Config conf;
catima::Projectile p{1,1,1,158.6};
{ //p->Be
catima::Material be = catima::get_material(4);
double r0 = catima::range(p,be);
conf.scattering = catima::scattering_types::fermi_rossi;
be.thickness(r0*0.01);
auto r1 = catima::calculate(p,be,conf);
be.thickness(r0*0.1);
auto r10 = catima::calculate(p,be,conf);
CHECK(r1.sigma_a*1000 == approx(2.93,0.1));
CHECK(r10.sigma_a*1000 == approx(9.49,0.4));
conf.scattering = catima::scattering_types::dhighland;
be.thickness(r0*0.01);
r1 = catima::calculate(p,be,conf);
be.thickness(r0*0.1);
r10 = catima::calculate(p,be,conf);
CHECK(r1.sigma_a*1000 == approx(2.04,0.1));
CHECK(r10.sigma_a*1000 == approx(7.61,0.3));
}
{ //p->Cu
catima::Material cu = catima::get_material(29);
double r0 = catima::range(p,cu);
conf.scattering = catima::scattering_types::fermi_rossi;
cu.thickness(r0*0.01);
auto r1 = catima::calculate(p,cu,conf);
cu.thickness(r0*0.1);
auto r10 = catima::calculate(p,cu,conf);
CHECK(r1.sigma_a*1000 == approx(7.23,0.3));
CHECK(r10.sigma_a*1000 == approx(23.5,0.8));
conf.scattering = catima::scattering_types::dhighland;
cu.thickness(r0*0.01);
r1 = catima::calculate(p,cu,conf);
cu.thickness(r0*0.1);
r10 = catima::calculate(p,cu,conf);
CHECK(r1.sigma_a*1000 == approx(5.63,0.25));
CHECK(r10.sigma_a*1000 == approx(20.7,0.8));
}
{ //p->Pb
catima::Material pb = catima::get_material(82);
pb.density(11.35);
pb.I(823);
double r0 = catima::range(p,pb);
conf.scattering = catima::scattering_types::fermi_rossi;
pb.thickness(r0*0.01);
auto r1 = catima::calculate(p,pb,conf);
pb.thickness(r0*0.1);
auto r10 = catima::calculate(p,pb,conf);
CHECK(r1.sigma_a*1000 == approx(11.88,0.5));
CHECK(r10.sigma_a*1000 == approx(38.6,1.2));
conf.scattering = catima::scattering_types::dhighland;
pb.thickness(r0*0.01);
r1 = catima::calculate(p,pb,conf);
pb.thickness(r0*0.1);
r10 = catima::calculate(p,pb,conf);
CHECK(r1.sigma_a*1000 == approx(9.75,0.25));
CHECK(r10.sigma_a*1000 == approx(35.8,1.2));
}
}
TEST_CASE("angular scattering"){
catima::Config conf;
catima::Projectile p{1,1,1,158.6};
catima::Material cu = catima::get_material(29);
catima::Material water = catima::get_material(catima::material::Water);
p.T = 158.6;
/*
for(double th = 0.01;th<3;th+=0.02){
cu.thickness(th);
conf.scattering = catima::scattering_types::fermi_rossi;
auto r1 = catima::calculate(p,cu,conf);
conf.scattering = catima::scattering_types::atima_scattering;
auto r2= catima::calculate(p,cu, conf);
CHECK( r2.sigma_a == approx(r1.sigma_a).R(1e-3));
}
*/
cu.thickness_cm(0.02963);
auto res = catima::calculate(p,cu);
CHECK( 1000*res.sigma_a == approx(7.2,0.5));
@ -351,7 +259,7 @@ using namespace std;
water.thickness_cm(29.4/n);
for(int i=0;i<n;i++)lll.add(water);
res2 = catima::calculate(p(215),lll);
CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.025));
CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
}
@ -425,8 +333,8 @@ 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);
CHECK(res1.sigma_a == approx(ra).R(1e-3));
auto ra = catima::angular_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
CHECK(res1.sigma_a == ra);
auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
CHECK(res1.sigma_E == re);
@ -496,7 +404,6 @@ using namespace std;
auto water = catima::get_material(catima::material::Water);
auto res2 = catima::calculate(p(600),water,600);
CHECK(res2.dEdxi == approx(92.5).epsilon(2));
CHECK(catima::radiation_length(water)==approx(36.1,0.2));
}
TEST_CASE("z_eff"){
using namespace catima;
@ -573,23 +480,4 @@ using namespace std;
CHECK(0.1*hbar*c_light/atomic_mass_unit == approx(0.21183,0.0001));
CHECK(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*PI) == approx(5.21721169334564e-7).R(1e-3));
}
TEST_CASE("phasespace"){
using namespace catima;
catima::Projectile p{1,1,1,250};
catima::Material graphite;
graphite.add_element(12,6,1);
graphite.density(2.0);
graphite.thickness_cm(1.0);
Phasespace ps;
ps.sigma_a = 0.01;
Layers l;
l.add(graphite);
auto res = calculate(p, ps, l);
CHECK(res.total_result.sigma_a == approx(0.012,0.002));
CHECK(res.total_result.cov == approx(1.23e-4,1e-5));
}

50
tests/test_storage.cpp
View File

@ -114,45 +114,27 @@ using catima::LN10;
}
TEST_CASE("energy table"){
catima::LogVArray<catima::max_datapoints> etable(catima::logEmin,catima::logEmax);
catima::EnergyTable<catima::max_datapoints> energy_table(catima::logEmin,catima::logEmax);
double step = (catima::logEmax - catima::logEmin)/(catima::max_datapoints-1);
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));
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));
}
TEST_CASE("indexing"){
double val, dif;
catima::LogVArray<catima::max_datapoints> etable(catima::logEmin,catima::logEmax);
catima::EnergyTable<catima::max_datapoints> energy_table(catima::logEmin,catima::logEmax);
CHECK(energy_table.index(0.0)==-1);
for(int i=0;i<catima::max_datapoints-1;i++){
val = energy_table[i];
dif = energy_table[i+1] - val;
CHECK(energy_table.index(val)==i);
CHECK(energy_table.index(val+0.5*dif)==i);
CHECK(energy_table.index(val)==i);
CHECK(energy_table.index(val+0.5*dif)==i);
CHECK(EnergyTable_index(catima::energy_table, 0.0)==-1);
CHECK(etable.index(val)==i);
CHECK(etable.index(val+0.5*dif)==i);
CHECK(etable.index(val+0.01*dif)==i);
CHECK(etable.index(val+0.99*dif)==i);
for(int i=0;i<catima::max_datapoints-1;i++){
val = catima::energy_table.values[i];
dif = catima::energy_table.values[i+1] - val;
CHECK(EnergyTable_index(catima::energy_table, val)==i);
CHECK(EnergyTable_index(catima::energy_table, val+0.5*dif)==i);
CHECK(catima::energy_table.index(val)==i);
CHECK(catima::energy_table.index(val+0.5*dif)==i);
}
}