From 2be81a4544d758e11ce1223d2aecc9f6913a2e7d Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Mon, 13 Aug 2018 18:08:54 +0200
Subject: [PATCH] production rate type

---
 CMakeLists.txt          |  6 +++---
 README.md               |  2 +-
 reactions.cpp           | 13 ++++++++++++-
 reactions.h             | 13 ++++++-------
 tests/test_reaction.cpp | 26 +++++++++++++++++++++++++-
 5 files changed, 47 insertions(+), 13 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5851ea7..5d3ee38 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/README.md b/README.md
index f3c0e83..1e08fd4 100644
--- a/README.md
+++ b/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)
diff --git a/reactions.cpp b/reactions.cpp
index 74efb1d..127670f 100644
--- a/reactions.cpp
+++ b/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
diff --git a/reactions.h b/reactions.h
index 66b2937..f40de79 100644
--- a/reactions.h
+++ b/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
diff --git a/tests/test_reaction.cpp b/tests/test_reaction.cpp
index 58c13dd..6627123 100644
--- a/tests/test_reaction.cpp
+++ b/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));
     }
 };