cwrapper update for mocadi

CMakeLists.txt

@@ -49,13 +49,15 @@ 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)
+list(APPEND EXTRA_LIBS fmt::fmt)
 
 configure_file(
                     "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
@@ -80,7 +82,7 @@ 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_include_directories(catima
                            PUBLIC  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
@@ -157,6 +159,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}

catima.cpp

@@ -137,14 +137,15 @@ 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 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);
+    spline_type range_spline = get_range_spline(data);    
-    double range = range_spline(T);
+    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,11 +182,11 @@ 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);    
+    t.thickness(th);
     return angular_straggling(p,t,c);
 }
 

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

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,51 +1,90 @@
 #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::Material mat;
-        catima::Projectile p(pa,pz);
         if(tz>200){
             mat = catima::get_material(tz);
         }
-        else{
+        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::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_angular_straggling_from_E(double pa, int pz, double Tin, double Tout,double ta, double tz){
+    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);          
-        catima::Material mat;
+        return energy_out(p(T), 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_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);
+
+        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);