Merge pull request #41 from hrosiak/reactions3

Reactions3

CMakeLists.txt

@@ -3,7 +3,7 @@ project(catima)
 
 ############ options #############
 #option(THREADS "Use multi-threading" ON)
-option(CATIMA_PYTHON "compile the Catima python module(requires numpy and cython installed)" OFF)
+option(PYTHON_MODULE "compile the Catima python module(requires numpy and cython installed)" OFF)
 option(TESTS "build tests" OFF)
 option(EXAMPLES "build examples" ON)
 option(GSL_INTEGRATION "use GSL integration" ON)
@@ -110,7 +110,7 @@ MESSAGE( STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER} )
 MESSAGE( STATUS "CMAKE_CXX_FLAGS: " ${CMAKE_CXX_FLAGS_RELEASE} )
 
 ######## for python module
-if(CATIMA_PYTHON)
+if(PYTHON_MODULE)
     if(NOT PYTHONINTERP_FOUND)
         MESSAGE(SEND_ERROR "Python is required to build nurex python modules") 
     endif(NOT PYTHONINTERP_FOUND)
@@ -139,7 +139,7 @@ if(CATIMA_PYTHON)
     add_custom_command(TARGET target_python
                         COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py  build_ext ${CYTHON_DEFINES} -i                         
                         )
-endif(CATIMA_PYTHON )
+endif(PYTHON_MODULE )
 
 ########## Sub Directories ###########
 if(EXAMPLES)

README.md

@@ -26,7 +26,7 @@ compile options, enable or disable with cmake:
 > cmake ../ -D[OPTION]
 
 available options:
-  * CATIMA_PYTHON - enable/disable building of the python bindigs, cython and numpy are required to build the catima python module, default OFF
+  * PYTHON_MODULE - enable/disable building of the python bindigs, cython and numpy are required to build the catima python module, default OFF
   * TESTS - build tests
   * EXAMPLES - build examples
   * DOCS - prepare doxygen documentation (after cmake, __make docs__ needs to be executed)

examples/makefile

@@ -1,8 +1,8 @@
-PROGRAMS=simple dedx example2 materials ls_coefficients
+PROGRAMS=simple dedx example2 materials ls_coefficients reactions
 
 GCC=g++ -Wall -std=c++14
 INCDIR=-I$(CATIMAPATH)/include
-LIBDIR=-L$(CATIMAPATH)
+LIBDIR=-L$(CATIMAPATH)/lib
 LIBS=-lcatima
 
 

examples/reactions.cpp

@@ -0,0 +1,27 @@
+#include "catima/catima.h"
+#include "catima/reactions.h"
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+
+int main(){
+    catima::Material target = catima::get_material(4);
+    target.thickness(1.0);      // thickness in g/cm2
+    catima::Projectile p(12,6); // define projectile, ie 12C
+    
+    double cs = 45;
+    double rcsi = 870;
+    double rcso = 860;
+    
+    cout<<"C->Be\n";
+    cout<<"t(g/cm2)\t rate"<<endl;
+    for(double t=0.25; t<=5;t+=0.25){
+        target.thickness(t);
+        double r = production_rate(45,rcsi, rcso, target);
+        cout<<t<<"\t"<<r<<endl;
+    }
+    
+    return 0;
+}

reactions.cpp

@@ -55,7 +55,17 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
     }
     return exp(-cs*0.0001);
     }
-
+    
+double production_rate(double cs, double rcs_projectile, double rcs_product, const Material &target, const Config &c){
+    double t = target.number_density_cm2();
+    double res = 0.0;
+    if( std::abs(rcs_product - rcs_projectile)<5.0){
+        res = 0.0001*cs*t*(exp(-0.0001*rcs_product*t));
+    }
+    else{
+        res = cs*(exp(-0.0001*rcs_projectile*t) - exp(-0.0001*rcs_product*t))/(rcs_product-rcs_projectile);
+    }
+    return res;
 }
 
 #ifndef NUREX
@@ -84,3 +94,4 @@ double SigmaR_Kox(int Ap, int Zp, double E, int At, int Zt){
     }
 #endif
 
+} //end of catima namespace

reactions.h

@@ -22,6 +22,7 @@
 #include "catima/config.h"
 #include "catima/integrator.h"
 #include <cmath>
+#endif
 
 namespace catima{
     
@@ -49,17 +50,14 @@ namespace catima{
         double i = ii.integrate(f,0,t);
         return 1.0 - std::exp(-i*0.0001);
     }
-
-    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);   
-}
-#else
-
+    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);
+    double production_rate(double cs, double rcs_projectile, double rcs_product, const Material &target, const Config &c=default_config);
+    
+#ifndef NUREX
 double SigmaR_Kox(int Ap, int Zp, double E, int At, int Zt);
-
 inline double p_from_T(double T, double M=1.0){
     return M*sqrt(T*T + 2*T*atomic_mass_unit);
 }
-
 inline double Ecm_from_T_relativistic(double T, double Ap, double At){
     double mp = Ap*atomic_mass_unit;
     double mt = At*atomic_mass_unit;
@@ -72,3 +70,4 @@ inline double Ecm_from_T_relativistic(double T, double Ap, double At){
 
 #endif //NUREX
 #endif
+} // end of catime namespace

tests/test_reaction.cpp

@@ -48,7 +48,31 @@ const lest::test specification[] =
         EXPECT( (r > 0 && r<1.0) );
         EXPECT( (r2 > 0 && r2<1.0) );
         EXPECT( r2>r );
-
+    },
+    CASE("production"){
+        catima::Projectile proj{12,6,6,870};
+        auto c = catima::get_material(6);
+        c.thickness(0.1);
+        double r,r2;
+        
+        double cs = 70;
+        double rcsi = 870;
+        double rcso = 850;
+        
+        r = 0.0001*cs*c.number_density_cm2();
+        r2 = catima::production_rate(cs,rcsi,rcso, c);
+        EXPECT(r==approx(r2).R(0.01));
+        
+        r2 = catima::production_rate(cs,2,1, c);
+        EXPECT(r==approx(r2).R(0.001));
+        
+        r = catima::production_rate(cs,870,870, c);
+        r2 = catima::production_rate(cs,870,860, c);
+        EXPECT(r==approx(r2).R(0.001));
+        c.thickness(2.0);
+        r = catima::production_rate(cs,870,870, c);
+        r2 = catima::production_rate(cs,870,860, c);
+        EXPECT(r==approx(r2).R(0.001));
     }
 };