Merge pull request #89 from hrosiak/cwrap

Cwrap

CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.12.0)
+cmake_minimum_required(VERSION 3.14)
 project(catima)
 
 ############ options #############
@@ -32,7 +32,6 @@ MESSAGE(STATUS "Build type: " ${CMAKE_BUILD_TYPE})
 
 ######### compiler flags ###########
 set(CMAKE_CXX_EXTENSIONS OFF)
-set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 MESSAGE(STATUS "install prefix: " ${CMAKE_INSTALL_PREFIX})
 
@@ -49,22 +48,32 @@ if(GSL_INTEGRATION OR GSL_INTERPOLATION)
     list(APPEND EXTRA_LIBS ${GSL_LIBRARIES} )
 endif()
 
-find_package(nurex QUIET)
+#find_package(nurex QUIET)
-if(nurex_FOUND)
+#if(nurex_FOUND)
-message(STATUS "nurex library found")
+#message(STATUS "nurex library found")
-set(NUREX ON)
+#set(NUREX ON)
-list(APPEND EXTRA_LIBS nurex::nurex)
+#list(APPEND EXTRA_LIBS nurex::nurex)
-endif(nurex_FOUND)
+#endif(nurex_FOUND)
 
-
+find_package(fmt QUIET)
-configure_file(
+if(NOT fmt_FOUND)
-                    "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
+message("fmt library not found, trying to dowload")
-                    "${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h"
+	    include(FetchContent)
+	    FetchContent_Declare(
+		    fmt
+		    GIT_REPOSITORY https://github.com/fmtlib/fmt.git
+		    GIT_TAG 8.1.1
 		    )
+      set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+      set(FMT_INSTALL ON)
+	    FetchContent_MakeAvailable(fmt)
+endif(NOT fmt_FOUND)
+
+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)
@@ -80,8 +89,8 @@ set_target_properties(catima PROPERTIES
                       LIBRARY_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,7 +103,7 @@ 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
+######## for python module #########
 find_package(Python COMPONENTS Interpreter Development)
 if(Python_FOUND)
     message(STATUS "Python found: ${Python_EXECUTABLE}")
@@ -105,7 +114,7 @@ if(PYTHON_MODULE)
     endif(NOT Python_FOUND)
     find_package(pybind11 QUIET)
     if(NOT pybind11_FOUND)
-	    message(INFO "pybind11 not found, trying to dowload")
+	    message("pybind11 not found, trying to dowload")
 	    include(FetchContent)
 	    FetchContent_Declare(
 		    pybind11
@@ -121,10 +130,12 @@ if(PYTHON_MODULE)
                                     $<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)
-	execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_SOURCE_DIR}/pymodule/setup.py bdist_wheel)
+	execute_process(COMMAND ${Python_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py bdist_wheel)
 endif(PYTHON_WHEEL)
 
 ########## Sub Directories ###########
@@ -157,6 +168,10 @@ endif(APPS)
 ####### install part #######
 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}

appveyor.yml

@@ -33,10 +33,10 @@ build_script:
 #  - git submodule add -b stable ../../pybind/pybind11 extern/pybind11
 #  - git submodule update --init
   - mkdir build && cd build
-#  - cmake -DBUILD_SHARED_LIBS=OFF -DAPPS=OFF -G "Visual Studio 16 2019" -A%PLATFORM% -DPYTHON_EXECUTABLE=$(python -c "import sys; print(sys.executable)") ../
+  - cmake -DBUILD_SHARED_LIBS=OFF -DAPPS=OFF -DPYTHON_WHEEL=OFF -G "Visual Studio 16 2019" -A%PLATFORM% ../
-  - cmake -DBUILD_SHARED_LIBS=OFF -DAPPS=OFF -G "Visual Studio 16 2019" -A%PLATFORM% ../
   - cmake --build ./ --config "%CONFIG%"
-  - python ../pymodule/setup.py bdist_wheel
+  - 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\*

catima.cpp

@@ -138,13 +138,14 @@ double Tfr(const Projectile &p, double X, double Es2){
 
 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;
@@ -181,10 +182,10 @@ 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 T, double Tout, Material t, const Config &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(T)-range_spline(Tout);
+    double th = range_spline(p.T)-range_spline(Tout);    
     t.thickness(th);
     return angular_straggling(p,t,c);
 }
@@ -210,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;

catima.h

@@ -138,7 +138,7 @@ namespace catima{
       * @param mat - Material
       * @return angular straggling
       */
-    double angular_straggling_from_E(const Projectile &p, double T, double Tout,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

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;
 

convert.h

@@ -0,0 +1,77 @@
+/*
+ *  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/>.
+ */
+
+#ifndef CONVERT_H
+#define CONVERT_H
+#include "catima/structures.h"
+#include "catima/material_database.h"
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include "fmt/format.h"
+
+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 = fmt::format("BEAM\n100000\n{}, 0, {}, {}\n2\n{}, {}, 0, 0, 0\n2\n{}, {}, 0, 0, 0\n1\n0,  0, 0, 0, 0\n",p.T,p.A,p.Z,psx.sigma_x,1000*psx.sigma_a, psy.sigma_x,1000*psy.sigma_a);
+    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 = fmt::format("*******\nMATTER\n{}, {}, {}\n2,{}\n0.\n0.,0.,0.\n1,1,0\n0,0,0,0\n",a,z,m.density()*1000,m.thickness_cm());
+        fw<<mstr;        
+        c++;
+    }
+    fw<<"ERWARTUNGSWERTE\nSAVE\nEND";
+    fw.close();
+    return true;
+}
+
+#endif

cwrapper.cpp

@@ -1,14 +1,15 @@
 #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};
 
-    CatimaResult catima_calculate(double pa, int pz, double T, double ta, double tz, double thickness, double density){
+    catima::Material make_material(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);
         }
@@ -16,36 +17,74 @@ extern "C" {
             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::Projectile p(pa,pz);
+        catima::Material mat = make_material(ta,tz, thickness, density);
         catima::Result r =  catima::calculate(p(T),mat);    
         CatimaResult res;
-        std::memcpy(&res,&r,sizeof(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::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::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::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::Projectile p(pa,pz);        
+        catima::Material mat = make_material(ta,tz, 0, -1);
 
-        catima::Material mat;
+        return catima::angular_straggling_from_E(p(Tin),Tout,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;
+        catima::Material mat = make_material(ta,tz, 0, -1);
-        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);
     }
 
+    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);
+    }
+
 }

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
 }

pymodule/pycatima.cpp

@@ -6,8 +6,10 @@
 #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;
 
@@ -467,6 +469,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);