From eb0db58dba6d81aad7d5da9afde977bfc85c0a89 Mon Sep 17 00:00:00 2001 From: hrocho Date: Tue, 29 Mar 2022 11:19:50 +0200 Subject: [PATCH] cwrapper update for mocadi --- CMakeLists.txt | 20 ++++++---- catima.cpp | 17 +++++---- catima.h | 2 +- convert.h | 77 +++++++++++++++++++++++++++++++++++++++ cwrapper.cpp | 85 +++++++++++++++++++++++++++++++------------ cwrapper.h | 8 +++- pymodule/pycatima.cpp | 3 ++ 7 files changed, 171 insertions(+), 41 deletions(-) create mode 100644 convert.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 27964c5..ecbe0c1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,13 +49,15 @@ 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) +#find_package(nurex QUIET) +#if(nurex_FOUND) +#message(STATUS "nurex library found") +#set(NUREX ON) +#list(APPEND EXTRA_LIBS nurex::nurex) +#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 $ @@ -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} diff --git a/catima.cpp b/catima.cpp index a2f7ecb..2b9809a 100644 --- a/catima.cpp +++ b/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); - double range = range_spline(T); - + 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,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); - t.thickness(th); + double th = range_spline(p.T)-range_spline(Tout); + t.thickness(th); return angular_straggling(p,t,c); } diff --git a/catima.h b/catima.h index 505a2c5..9667b13 100644 --- a/catima.h +++ b/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 diff --git a/convert.h b/convert.h new file mode 100644 index 0000000..874039a --- /dev/null +++ b/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 . + */ + +#ifndef CONVERT_H +#define CONVERT_H +#include "catima/structures.h" +#include "catima/material_database.h" +#include +#include +#include +#include +#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 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::default_config.z_effective = catima_defaults.z_effective; - catima::Material mat; - catima::Projectile p(pa,pz); + catima::Material make_material(double ta, double tz, double thickness, double density){ + catima::Material mat; 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; - if(tz>200){ - mat = catima::get_material(tz); - } - else{ - mat.add_element(ta,tz,1.0); - } + catima::Material mat = make_material(ta,tz, thickness, density); + return energy_out(p(T), mat); + } - 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; - if(tz>200){ - mat = catima::get_material(tz); - } - else{ - mat.add_element(ta,tz,1.0); - } + 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); + } + } \ No newline at end of file diff --git a/cwrapper.h b/cwrapper.h index 90ae782..ed17974 100644 --- a/cwrapper.h +++ b/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 } diff --git a/pymodule/pycatima.cpp b/pymodule/pycatima.cpp index 50024ca..f2bfe7a 100644 --- a/pymodule/pycatima.cpp +++ b/pymodule/pycatima.cpp @@ -6,8 +6,10 @@ #include "catima/catima.h" #include "catima/srim.h" #include "catima/nucdata.h" +#include "catima/convert.h" #include #include + 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);