diff --git a/CMakeLists.txt b/CMakeLists.txt index 8897d27..203e292 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -86,7 +86,7 @@ set_target_properties(catima_static POSITION_INDEPENDENT_CODE ON ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib ) -target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES} ) +target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES} nurex::nurex) target_link_libraries(catima_static ${EXTRA_LIBS} ${GSL_LIBRARIES}) target_include_directories(catima diff --git a/reactions.h b/reactions.h index b856539..8473993 100644 --- a/reactions.h +++ b/reactions.h @@ -21,9 +21,37 @@ #include "nurex/nurex.h" #include "catima/structures.h" +#include "catima/config.h" +#include "catima/integrator.h" +#include namespace catima{ + /** + * return reaction probability + * @param sigma - cross section in mb + * @param t - number of targets per cm2 in 10^23 unit + */ + inline double reaction_rate(double sigma, double t){ + return 1.0 - std::exp(-sigma*t*0.0001); + } + + /** + * return nonreaction rate + * @param sigma - cross section in mb + * @param t - number of targets per cm2 in 10^23 unit + */ + inline double nonreaction_rate(double sigma, double t){ + return std::exp(-sigma*t*0.0001); + } + + template + double reaction_rate(F& f, double t){ + GaussLegendreIntegration<8> ii; + double i = ii.integrate(f,0,t); + return 1.0 - std::exp(-i*0.0001); + } + double reaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config); } diff --git a/structures.cpp b/structures.cpp index e10b067..f46690d 100644 --- a/structures.cpp +++ b/structures.cpp @@ -27,15 +27,19 @@ bool operator==(const Material &a, const Material&b){ } -Material::Material(std::initializer_list>list,double _density, double _ipot):rho(_density),i_potential(_ipot){ +Material::Material(std::initializer_list>list,double _density, double _ipot, double mass):rho(_density),i_potential(_ipot){ std::initializer_list>::iterator it; atoms.reserve(list.size()); for ( it=list.begin(); it!=list.end(); ++it){ add_element((*it)[0],(*it)[1],(*it)[2]); } - calculate(); // calculate if needed, ie average molar mass - + if(mass!=0.0){ + molar_mass=mass; + } + else{ + calculate(); // calculate if needed, ie average molar mass + } } Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th){ diff --git a/structures.h b/structures.h index b6053a4..58d824a 100644 --- a/structures.h +++ b/structures.h @@ -20,6 +20,7 @@ #include #include #include +#include "catima/constants.h" namespace catima{ @@ -89,7 +90,7 @@ namespace catima{ }); * \endcode */ - Material(std::initializer_list>list,double _density=0.0, double ipot = 0.0); + Material(std::initializer_list>list,double _density=0.0, double ipot = 0.0, double mass=0.0); /** * calculates internal variables if needed @@ -132,6 +133,11 @@ namespace catima{ * @return returns Molar Mass of the Material */ double M() const {return molar_mass;} + + /** + * sets molar mass of the Material + */ + Material& M(double mass){molar_mass=mass; return *this;} /** * @return returns density in g/cm^3 @@ -153,6 +159,11 @@ namespace catima{ */ Material& thickness(double val){th = val;return *this;}; + /** + * set length in cm, density should be set before + */ + Material& thickness_cm(double l){th = rho*l; return *this;} + /** * set the mean ionization potential, if non elemental I should be used */ @@ -164,6 +175,38 @@ namespace catima{ */ double I() const {return i_potential;}; + /** + * return number density of atoms/molecules per cm3 in 10^23 units + */ + double number_density()const{ + return Avogadro*rho/molar_mass; + } + + /** + * return number density of atoms of i-th element in 10^23 units + */ + double number_density(int i)const{ + if(i>=atoms.size())return 0.0; + return number_density()*molar_fraction(i); + } + + /** + * return number density of atoms/molecules per cm2 in 10^23 units + */ + double number_density_cm2()const{ + return Avogadro*th/molar_mass; + } + + /** + * return number density of atoms per cm2 of i-th element in 10^23 units + */ + double number_density_cm2(int i)const{ + if(i>=atoms.size())return 0.0; + return number_density_cm2()*molar_fraction(i); + } + + + friend bool operator==(const Material &a, const Material&b); }; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4e490c7..0c33416 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -10,6 +10,11 @@ foreach(entry ${CATIMA_TESTS}) add_test(${entry} ${PROJECT_BINARY_DIR}/tests/${entry}) endforeach(entry in ${CATIMA_TESTS}) +if(NUREX) + add_executable(test_reaction test_reaction.cpp) + target_link_libraries(test_reaction nurex::nurex catima) + add_test(test_reaction ${PROJECT_BINARY_DIR}/tests/test_reaction) +endif(NUREX) set(C_TESTS test_c) MESSAGE( STATUS "CMAKE_C_COMPILER: " ${CMAKE_C_COMPILER} ) diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index cbaa505..baa713d 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -313,6 +313,13 @@ const lest::test specification[] = EXPECT(res.Eout == approx(345.6).epsilon(1.0)); EXPECT(res.sigma_a == approx(0.0013).epsilon(1e-4)); EXPECT(res.sigma_E == approx(0.12).epsilon(1e-3)); + EXPECT(res.dEdxi == approx(103.5).epsilon(1e-1)); + + res = catima::calculate(p(150),air); + EXPECT(res.dEdxi == approx(173.6).epsilon(1e0)); + res = catima::calculate(p(1000),air); + EXPECT(res.dEdxi == approx(70.69).epsilon(1e-0)); + auto water = catima::get_material(catima::material::Water); auto res2 = catima::calculate(p(600),water,600); diff --git a/tests/test_reaction.cpp b/tests/test_reaction.cpp new file mode 100644 index 0000000..3a6b0c4 --- /dev/null +++ b/tests/test_reaction.cpp @@ -0,0 +1,34 @@ +#include "lest.hpp" +#include +#include "catima/catima.h" +#include "catima/reactions.h" +#include "testutils.h" +using namespace std; +using catima::approx; +using catima::reaction_rate; +using catima::nonreaction_rate; +const lest::test specification[] = +{ + CASE("reaction"){ + auto c = catima::get_material(6); + c.thickness(4.0); + double r = reaction_rate(1000,c.number_density_cm2()); + std::cout<