Fix typo in gwm_integrators

energy table indexing

.gitignore

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

CMakeLists.txt

@@ -1,27 +1,28 @@
-cmake_minimum_required(VERSION 3.2.0)
+cmake_minimum_required(VERSION 3.14)
 project(catima)
 
 ############ options #############
-option(BUILD_SHARED_LIBS "build as shared library" ON)
+option(BUILD_SHARED_LIBS "build as shared library" OFF)
 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" ON)
+option(APPS "build catima applications" OFF)
 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}")
+    #set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -ffast-math")
     MESSAGE(STATUS "Build type Release")
 else()
-    set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")
+    #set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")
     if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang" OR ${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
         set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare")
     endif()
@@ -31,8 +32,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})
 
@@ -49,39 +50,30 @@ 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 *.cpp)
+file(GLOB SOURCES catima/*.cpp)
 if(GLOBAL)
 file(GLOB GLOBAL_SOURCES global/*.c)
 LIST (APPEND SOURCES ${GLOBAL_SOURCES})
 endif(GLOBAL)
-file(GLOB HEADERS *.h libs/*.h)
+file(GLOB HEADERS catima/*.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
                       )
 
-target_link_libraries(catima ${EXTRA_LIBS})
-
+target_link_libraries(catima PUBLIC ${EXTRA_LIBS})
+target_compile_features(catima PRIVATE cxx_std_17)
 target_include_directories(catima
                            PUBLIC  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
                                    $<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
@@ -94,26 +86,43 @@ 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(PythonInterp)
-if(PYTHONINTERP_FOUND)
-    message(STATUS "Python found: ${PYTHON_EXECUTABLE}")
+######## for python module #########
+find_package(Python COMPONENTS Interpreter Development)
+if(Python_FOUND)
+    message(STATUS "Python found: ${Python_EXECUTABLE}")
 endif()
 if(PYTHON_MODULE)
-    if(NOT PYTHONINTERP_FOUND)
+    if(NOT Python_FOUND)
         MESSAGE(SEND_ERROR "Python is required to build nurex python modules")
-    endif(NOT PYTHONINTERP_FOUND)
-    find_package(pybind11 REQUIRED)
-    set(PYBIND11_CPP_STANDARD -std=c++14)
-    pybind11_add_module(pycatima pymodule/pycatima)
+    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)
     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)
-
+    target_link_libraries(pycatima PRIVATE catima fmt::fmt)    
 endif(PYTHON_MODULE )
 
+configure_file("${PROJECT_SOURCE_DIR}/pymodule/setup.py.in" "${PROJECT_BINARY_DIR}/setup.py")
+if(PYTHON_WHEEL)
+    check_fmt()
+	execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py bdist_wheel)
+endif(PYTHON_WHEEL)
+
 ########## Sub Directories ###########
 if(EXAMPLES)
     file(GLOB EXAMPLES examples/*.cpp)
@@ -142,8 +151,12 @@ add_subdirectory("bin")
 endif(APPS)
 
 ####### install part #######
-FILE(GLOB headers "*.h")
+FILE(GLOB headers "catima/*.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}

README.md

@@ -1,5 +1,6 @@
 [![Codacy Badge](https://api.codacy.com/project/badge/Grade/dc251db65f7a4c06ae07380544ea08fc)](https://www.codacy.com/manual/hrosiak/catima?utm_source=github.com&utm_medium=referral&utm_content=hrosiak/catima&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
@@ -32,6 +33,13 @@ 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:

appveyor.yml

@@ -0,0 +1,42 @@
+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\*

build_config.in

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

catima/calculations.cpp

@@ -1,4 +1,4 @@
-#include <math.h>
+#include <cmath>
 #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 = pow(p.Z,0.23)+pow(t.Z,0.23);
+    double zpowers = std::pow(p.Z,0.23)+std::pow(t.Z,0.23);
     double asum = p.A + t.A;
     return 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
 }
@@ -473,17 +473,44 @@ double p_from_T(double T, double M){
 double energy_straggling_firsov(double z1,double energy, double z2, double m2){
     double gamma = gamma_from_T(energy);
     double beta2=1.0-1.0/(gamma*gamma);
-    double factor=4.8184E-3*pow(z1+z2,8.0/3.0)/m2;
+    double factor=4.8184E-3*std::pow(z1+z2,8.0/3.0)/m2;
     return factor*beta2/fine_structure/fine_structure;
     }
 
-double angular_scattering_variance(const Projectile &p, const Target &t){
+double angular_scattering_power(const Projectile &p, const Target &t, 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 lr = radiation_length(t.Z,t.A);
-    return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
+    //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));
 }
 
 /// radiation lengths are taken from Particle Data Group 2014
@@ -515,7 +542,7 @@ double radiation_length(int z, double m){
     if(z==92){return 6.00;}
 
     double z2 = z*z;
-    double z_13 = 1.0/pow(z,1./3.);
+    double z_13 = 1.0/std::pow(z,1./3.);
     double z_23 = z_13*z_13;
     double a2 = fine_structure*fine_structure*z2;
     double a4 = a2*a2;
@@ -524,6 +551,31 @@ 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;
@@ -567,11 +619,11 @@ double dedx_variance(const Projectile &p, const Target &t, const Config &c){
     double beta = beta_from_T(p.T);
     double beta2 = beta*beta;
     double zp_eff = z_effective(p,t,c);
-    double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
+    double f = domega2dx_constant*ipow(zp_eff,2)*t.Z/t.A;
 
     if( (c.calculation == omega_types::atima) ){
-        cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
-			log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364)));
+        cor = 24.89 * std::pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
+			log( 2.0*electron_mass*1e6*beta2/(33.05*std::pow(t.Z,1.6364)));
 	    cor = std::max(cor, 0.0 );
     }
 	//double X = bethek_lindhard_X(p);
@@ -613,18 +665,19 @@ 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/pow(z,2.0/3.0)));
+    return z*(1.0-exp(-0.95*fine_structure_inverted*beta/std::pow(z,2.0/3.0)));
 }
 
 double z_eff_Anthony_Landford(double pz, double beta, double tz){
     double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz);
     double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz);
-    return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0))));
+    return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/std::pow(pz,2.0/3.0))));
 }
 
 double z_eff_Hubert(double pz, double E, double tz){
@@ -652,7 +705,7 @@ double z_eff_Hubert(double pz, double E, double tz){
         x4 = 0.5218 + 0.02521*lntz;
     }
 
-    return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4)));
+    return pz*(1-x1*exp(-x2*std::pow(E,x3)*std::pow(pz,-x4)));
 }
 
 double z_eff_Winger(double pz, double beta, double tz){
@@ -680,7 +733,7 @@ double z_eff_Winger(double pz, double beta, double tz){
     if(tz==2){
         tz = 2.8;
     }
-	x = beta /0.012 /pow(pz,0.45);
+	x = beta /0.012 /std::pow(pz,0.45);
 	lnz =log(pz);
 	lnzt=log(tz);
 
@@ -696,13 +749,17 @@ 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
+    else if(E<30.0 || pz<29){        
+        return z_eff_Pierce_Blann(pz, E);
+    }
+    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){
@@ -710,10 +767,10 @@ double z_eff_Schiwietz(double pz, double beta, double tz){
     double c1, c2;
     double x;
 
-    scaled_velocity = catima::power(pz,-0.543)*beta/bohr_velocity;
+    scaled_velocity = std::pow(pz,-0.543)*beta/bohr_velocity;
     c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9);
     c2 = 1+0.030*scaled_velocity*log(tz);
-    x = c1*catima::power(scaled_velocity/c2/1.54,1+(1.83/pz));
+    x = c1*std::pow(scaled_velocity/c2/1.54,1+(1.83/pz));
     return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) );
 
 }
@@ -733,17 +790,17 @@ double z_eff_atima14(double pz, double T, double tz){
         double c2 = 0.28;
         double c3 = 0.04;
         qglobal = z_eff_global(pz,T,tz);
-        qglobal = (qglobal - qpb)*c1/catima::power(tz+1,c2)*(1-exp(-c3*T)) + qpb;
+        qglobal = (qglobal - qpb)*c1/std::pow(tz+1,c2)*(1-exp(-c3*T)) + qpb;
     }
 
     emax = 1500.;
     emin = 1000.;
     if(T>emax){
-        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
+        qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
         qmean = qhigh;
     }
     else if(T>=emin && T<=emax){
-        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
+        qhigh = pz * (1.0-exp(-180.0*beta*std::pow(gamma,0.18)*std::pow(pz,-0.82)*std::pow(tz,0.1)));
         if(pz<=28){
             qwinger = z_eff_Winger(pz,beta,tz);
             qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin);

catima/calculations.h

@@ -100,9 +100,25 @@ 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;
     
-    double angular_scattering_variance(const Projectile &p, const Target &t);
+    /**
+     * 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);
     /**
       * returns radiation length of the (M,Z) material
       * for certain z the radiation length is tabulated, otherwise calculated
@@ -112,7 +128,27 @@ namespace catima{
       */
     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
       */ 
@@ -194,5 +230,13 @@ namespace catima{
     inline double power(double x, double y){
         return exp(log(x)*y);
     }
+
+    inline double ipow(double x, int y){
+      if(y==0)return 1.0;
+      else if(y==1) return x;      
+      else if(y==2) return x*x;      
+      else return std::pow(x,y);
+    }
+
 }
 #endif

catima/catima.cpp

@@ -1,5 +1,5 @@
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include <algorithm>
 #include "catima/catima.h"
 #include "catima/constants.h"
@@ -57,18 +57,6 @@ 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);
@@ -120,26 +108,90 @@ 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 angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
+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);    
-    //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));
+    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){
+    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);
 }
 
 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);    
-
-    //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);
@@ -159,7 +211,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_epsilon)return e;
+        if(fabs(r)<Eout_th_epsilon)return e;
         double step = -r*dedx;
         e = e-step;
         if(e<Ezero)return 0.0;
@@ -220,17 +272,29 @@ 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.Eout = energy_out(T,t.thickness(),range_spline);
+    res.sigma_r = sqrt(range_straggling_spline(T));        
+
+    if(t.thickness()==0){
+        res.dEdxo = res.dEdxi;
+	res.Eout = res.Ein;
+        return res;
+    }
+    
+    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;
@@ -238,6 +302,7 @@ Result calculate(Projectile p, const Material &t, const Config &c){
         res.sigma_E = 0.0;
     }    
     else{        
+        spline_type angular_variance_spline = get_angular_variance_spline(data);
         res.dEdxo = p.A/range_spline.derivative(res.Eout);        
         #ifdef THIN_TARGET_APPROXIMATION
         if(thin_target_limit*res.Ein<res.Eout){
@@ -245,59 +310,56 @@ Result calculate(Projectile p, const Material &t, const Config &c){
             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);
             }        
+        }
         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));
             }
             
+        }        
         #else
         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);
-        spline_type angular_variance_spline = get_angular_variance_spline(data);
+        if(use_angular_spline){
             res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
+        }        
         #endif
-        if( t.thickness()>0){
-            //auto tofdata = calculate_tof(p,t,c);
+        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);
-            }
-    }
-    res.sigma_r = sqrt(range_straggling_spline(T));    
-    res.Eloss = (res.Ein - res.Eout)*p.A;
             
+    } //end of else for non stopped case
+    
+    // position straggling in material    
     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 = angular_variance(p(T),t,c,2);
     res.sigma_x = sqrt(res.sigma_x);
 
-    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();
-            };
+    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);
 
-    //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 Layers &layers, const Config &c){
+MultiResult calculate(const Projectile &p, const Phasespace &ps, 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);
@@ -310,7 +372,6 @@ MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c
         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
@@ -321,13 +382,13 @@ MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c
     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(res.total_result.sigma_x);
+        res.total_result.sigma_x = sqrt(std::abs(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(res.total_result.sigma_x);
+        res.total_result.sigma_x = sqrt(std::abs(res.total_result.sigma_x));
         }
     return res;
 }
@@ -349,6 +410,9 @@ 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
@@ -366,8 +430,10 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     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];
+        //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 = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
         res = p.A*res;

catima/catima.h

@@ -113,6 +113,23 @@ 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
@@ -121,7 +138,7 @@ namespace catima{
       * @param mat - Material
       * @return angular straggling
       */
-    double angular_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
+    double angular_straggling_from_E(const Projectile &p, double Tout,Material t, const Config &c=default_config);
 
     /**
       * calculates Energy straggling in the material from difference of incoming a nd outgoing energies
@@ -186,7 +203,16 @@ namespace catima{
       * @return results stored in MultiResult structure
       *
       */
-    MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
+    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);
+    };
     inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
         return calculate(p(T), layers, c);
     }

catima/config.h

@@ -45,6 +45,16 @@ 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
       */
@@ -57,7 +67,8 @@ namespace catima{
 
         unsigned char corrections = 0;
         unsigned char calculation = 1;
-        unsigned char low_energy = 0;
+        unsigned char low_energy = low_energy_types::srim_85;
+        unsigned char scattering = scattering_types::atima_scattering;        
     };
 
 

catima/constants.h

@@ -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_epsilon = 1e-5;  //
+constexpr double Eout_th_epsilon = 1e-5;  //
 
 constexpr double thin_target_limit = 1 - 1e-3;
 

catima/convert.h

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

catima/cwrapper.cpp

@@ -0,0 +1,96 @@
+#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);
+    }
+
+}

catima/cwrapper.h

@@ -34,14 +34,18 @@ struct CatimaConfig {
         char z_effective;
 };
 
-struct CatimaConfig catima_defaults = {none};
+extern struct CatimaConfig catima_defaults;
 
 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
 }

catima/gwm_integrators.cpp

@@ -0,0 +1,97 @@
+/*
+ * 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;
+            }
+        }
+    }
+}

catima/gwm_integrators.h

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

catima/integrator.h

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

catima/material_database.cpp

@@ -158,6 +158,7 @@ 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();

catima/material_database.h

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

catima/spline.h

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

catima/storage.cpp

@@ -19,7 +19,11 @@
 #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)){

catima/storage.h

@@ -29,7 +29,7 @@
 
 #include "catima/spline.h"
 
-
+//#define VETABLE
 namespace catima{
 
     /**
@@ -48,36 +48,31 @@ namespace catima{
     static constexpr int size() {return N;};
 	double values[N];
 	double step;
-	double* begin(){return values;}
-    double* end(){return &values[num];}
+	const double* begin()const{return values;}
+    const double* end()const{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;
     };
 
-	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){
+	template<typename T>
+	double EnergyTable_interpolate(const T &table, double xval, double *y){
 	    double r;
-	    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);
+	    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);
 	    double linstep = table.values[i+1] - table.values[i];
         if(linstep == 0.0)return table.values[i];
 	    double x = 1.0 - ((xval - table.values[i])/linstep);
@@ -85,6 +80,63 @@ namespace catima{
 	    return r;
 	}
 
+    template <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
@@ -107,12 +159,13 @@ 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.values[0]), max(table.values[max_datapoints-1]), ss(table,y){}
+            min(table[0]), max(table[max_datapoints-1]), ss(table,y){}
         double operator()(double x)const{return eval(x);}
         double eval(double x)const{return ss.evaluate(x);}
         double derivative(double x)const{return ss.deriv(x);}
@@ -128,7 +181,14 @@ namespace catima{
 #ifdef GSL_INTERPOLATION
 using Interpolator = InterpolatorGSL;
 #else
-using Interpolator = InterpolatorCSpline;
+#ifdef VETABLE
+//using Interpolator = InterpolatorSplineT<LogVArray<max_datapoints>>;
+using Interpolator = InterpolatorCSpline<LogVArray<max_datapoints>>;
+#else 
+//using Interpolator = InterpolatorSplineT<EnergyTable<max_datapoints>>;
+using Interpolator = InterpolatorCSpline<EnergyTable<max_datapoints>>;
+#endif
+
 #endif
 
 #ifdef STORE_SPLINES

catima/structures.h

@@ -243,6 +243,12 @@ namespace catima{
         #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
       */

cwrapper.cpp

@@ -1,51 +0,0 @@
-#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);
-    }
-
-}

docs/README.md

@@ -18,6 +18,12 @@ 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
 -------------

docs/pycatima.md

@@ -1,5 +1,13 @@
 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
 ---------
@@ -158,7 +166,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 MeV/u
+  * sigma_a - Angular straggling in rad
   * sigma_r - range straggling in g/cm^2
   * sigma_x - position straggling in cm
   * sp - non-reaction probability

gwm_test/Makefile

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

gwm_test/gwm_test.cpp

@@ -0,0 +1,22 @@
+/*
+ * 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;
+
+}

pymodule/pycatima.cpp

@@ -6,16 +6,15 @@
 #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;
 
-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 catima_info(){
+    return "CATIMA version = 1.7\n";
 }
 
 std::string  material_to_string(const Material &r){
@@ -50,9 +49,14 @@ py::list storage_info(){
 
 py::list get_energy_table(){
     py::list r;
-    for(auto e : energy_table){
-        r.append(e);
+    for (size_t i = 0; i < energy_table.size(); i++)
+    {
+        r.append(energy_table[i]);
     }
+    
+    //for(auto e : energy_table){
+        //r.append(e);
+    //}
     return r;
 }
 
@@ -101,6 +105,7 @@ 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;
@@ -127,6 +132,12 @@ 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")
@@ -176,9 +187,13 @@ 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("get_dict",&get_result_dict)
+             .def("__repr__",[](const Result &self){
+                 return py::str(get_result_dict(self)); 
+             });
 
      py::class_<MultiResult>(m,"MultiResult")
              .def(py::init<>(),"constructor")
@@ -214,6 +229,16 @@ 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")
@@ -240,6 +265,13 @@ PYBIND11_MODULE(pycatima,m){
             .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)
@@ -262,7 +294,133 @@ PYBIND11_MODULE(pycatima,m){
 		 .value("Suprasil", material::Suprasil)
 		 .value("HAVAR", material::HAVAR)
 		 .value("Steel", material::Steel)
-            .value("CO2", material::CO2);
+		 .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);
+
 
 
     py::class_<Config>(m,"Config")
@@ -271,12 +429,14 @@ 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){
@@ -285,12 +445,17 @@ 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);
@@ -309,6 +474,7 @@ 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);

pymodule/setup.py

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

pymodule/setup.py.in

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

spline.cpp

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

tests/test_calculations.cpp

@@ -183,12 +183,104 @@ 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));
@@ -259,7 +351,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.01));
+        CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.025));
         }
       
 
@@ -333,8 +425,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 == ra);
+        auto ra = catima::angular_straggling_from_E(p(res1.Ein),res1.Eout,graphite);
+        CHECK(res1.sigma_a == approx(ra).R(1e-3));
 
         auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
         CHECK(res1.sigma_E == re);
@@ -404,6 +496,7 @@ 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;
@@ -481,3 +574,22 @@ using namespace std;
         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));
+    }
+    

tests/test_storage.cpp

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