reactions

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 

reactions.h

@@ -21,9 +21,37 @@
 
 #include "nurex/nurex.h"
 #include "catima/structures.h"
+#include "catima/config.h"
+#include "catima/integrator.h"
+#include <cmath>
 
 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<typename F>
+    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);
     
 }

structures.cpp

@@ -27,15 +27,19 @@ bool operator==(const Material &a, const Material&b){
 }
 
 
-Material::Material(std::initializer_list<std::array<double,3>>list,double _density, double _ipot):rho(_density),i_potential(_ipot){
+Material::Material(std::initializer_list<std::array<double,3>>list,double _density, double _ipot, double mass):rho(_density),i_potential(_ipot){
     std::initializer_list<std::array<double,3>>::iterator it;
     atoms.reserve(list.size());
     for ( it=list.begin(); it!=list.end(); ++it){
         add_element((*it)[0],(*it)[1],(*it)[2]);
     }
 
+    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){

structures.h

@@ -20,6 +20,7 @@
 #include <vector>
 #include <array>
 #include <initializer_list>
+#include "catima/constants.h"
 
 namespace catima{
 
@@ -89,7 +90,7 @@ namespace catima{
                 });
               * \endcode
               */
-            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0);
+            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0, double mass=0.0);
             
             /**
              * calculates internal variables if needed
@@ -133,6 +134,11 @@ namespace catima{
               */
             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);
     };

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} )

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);

tests/test_reaction.cpp

@@ -0,0 +1,34 @@
+#include "lest.hpp"
+#include <math.h>
+#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<<r<<"\n";
+
+        EXPECT(reaction_rate(0,10.0) == 0.0);
+        EXPECT(reaction_rate(1000,0.0) == 0.0);
+        
+        EXPECT(nonreaction_rate(0,10.0) == 1.0);
+        EXPECT(nonreaction_rate(1000,0.0) == 1.0);
+
+        auto fcc = [](double x){return 1000.0;};
+        r = reaction_rate(fcc, c.number_density_cm2());
+        std::cout<<r<<"\n";
+    }
+};
+
+int main( int argc, char * argv[] )
+{
+    return lest::run( specification, argc, argv );
+}
+

tests/test_structures.cpp

@@ -57,6 +57,15 @@ const lest::test specification[] =
             SECTION("default ionisation potential"){
               EXPECT(graphite.I()==0.0);
             }
+            SECTION("length"){
+                water.density(1.0);
+                water.thickness(1.0);
+                EXPECT(water.thickness()==approx(1.0,0.0001));
+                water.thickness_cm(1.0);
+                EXPECT(water.thickness()==approx(1.0,0.0001));
+                water.thickness_cm(2.0);
+                EXPECT(water.thickness()==approx(2.0,0.0001));
+            }
         }
     },
     CASE("Material automatic atomic weight"){
@@ -242,6 +251,8 @@ const lest::test specification[] =
                         {0, 8, 1}
                         });
 
+        auto air = catima::get_material(catima::material::Air);
+
         EXPECT(water1.weight_fraction(0)==0.111898);
         EXPECT(water2.weight_fraction(0)==approx(water1.weight_fraction(0)).R(1e-5));
         EXPECT(water1.weight_fraction(1)==0.888102);
@@ -263,6 +274,42 @@ const lest::test specification[] =
         EXPECT(mat.M()==approx(12.0,0.001));
         EXPECT(mat.weight_fraction(0)==approx(1.0).R(1e-6));
 
+        EXPECT(air.M() == approx(28.97,0.1));
+
+    },
+    CASE("number density"){
+        catima::Material c({12.0,6,1});
+        auto water = catima::get_material(catima::material::Water);
+        auto air = catima::get_material(catima::material::Air);
+        water.density(0.9982);
+        c.density(3.513);
+        air.density(1.2041e-3);
+
+        c.thickness_cm(1.0);
+        EXPECT(c.number_density()==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2()==approx(1.7662,0.01));
+        EXPECT(c.number_density(0)==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2(0)==approx(1.7662,0.01));
+        EXPECT(c.number_density(1)==0.0);
+        EXPECT(c.number_density_cm2(1)==0.0);
+        c.thickness_cm(2.0);
+        EXPECT(c.number_density()==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2()==approx(2.0*1.7662,0.01));
+        
+        water.thickness_cm(1.0);
+        EXPECT(water.number_density()==approx(0.3336,0.001));
+        EXPECT(water.number_density_cm2()==approx(0.3336,0.001));
+        EXPECT(water.number_density(0)==approx(2*0.3336,0.001));
+        EXPECT(water.number_density_cm2(0)==approx(2*0.3336,0.001));
+        EXPECT(water.number_density(1)==approx(0.3336,0.001));
+        EXPECT(water.number_density_cm2(1)==approx(0.3336,0.001));
+        water.thickness_cm(3.0);
+        EXPECT(water.number_density_cm2()==approx(3.0*0.3336,0.001));
+    
+        air.thickness_cm(1.0);
+        EXPECT(air.number_density(0)==approx(air.molar_fraction(0)*2*0.0002504,0.00001));
+        EXPECT(air.number_density(1)==approx(air.molar_fraction(1)*2*0.0002504,0.00001));
+        EXPECT(air.number_density(2)==approx(air.molar_fraction(2)*1*0.0002504,0.00001));
     }
 };