From 4c54ed0dd0c00081a402c2a1ad956017a7aa7455 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Tue, 31 Jul 2018 17:40:25 +0200
Subject: [PATCH] reactions v1

---
 CMakeLists.txt          |  2 ++
 build_config.in         |  2 +-
 catima.cpp              |  6 +++---
 config.h                |  2 +-
 constants.h             |  2 +-
 reactions.cpp           | 39 +++++++++++++++++++++++++++++++++++----
 reactions.h             | 20 ++++++++++++++++++--
 structures.h            |  2 +-
 tests/CMakeLists.txt    |  5 ++---
 tests/test_reaction.cpp |  2 ++
 10 files changed, 66 insertions(+), 16 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0d67399..043e80f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,6 +11,7 @@ option(GENERATE_DATA "make data tables generator" OFF)
 option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
 option(DOCS "build documentation (requires doxygen)" OFF)
 option(GLOBAL "build with global, sources are required" OFF)
+option(REACTIONS "enable/disable nuclear reaction rate" ON)
 option(APPS "build catima applications" ON)
 
 ######## build type ############
@@ -148,6 +149,7 @@ if(EXAMPLES)
     #add_subdirectory("examples")
 endif(EXAMPLES)
 if(TESTS)
+enable_testing()
 add_subdirectory("tests")
 endif(TESTS)
 
diff --git a/build_config.in b/build_config.in
index 2f2a6d0..b770cb4 100644
--- a/build_config.in
+++ b/build_config.in
@@ -4,7 +4,7 @@
 #cmakedefine THIN_TARGET_APPROXIMATION
 #cmakedefine GSL_INTEGRATION
 #cmakedefine GLOBAL
-
+#cmakedefine REACTIONS
 #cmakedefine NUREX
 
 #endif
diff --git a/catima.cpp b/catima.cpp
index 9ce17bb..73b08a8 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -9,7 +9,7 @@
 #include "catima/storage.h"
 #include "catima/nucdata.h"
 #include "catima/calculations.h"
-#ifdef NUREX
+#ifdef REACTIONS
 #include "catima/reactions.h"
 #endif
 
@@ -254,7 +254,7 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
     }
     res.sigma_r = sqrt(range_straggling_spline(T));
     res.Eloss = (res.Ein - res.Eout)*p.A;
-    #ifdef NUREX
+    #ifdef REACTIONS
     res.sp = nonreaction_rate(p,t,c);
     #endif
     return res;
@@ -274,7 +274,7 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
         res.total_result.sigma_E += r.sigma_E*r.sigma_E; 
         res.total_result.tof += r.tof;
         res.total_result.Eout = r.Eout;
-        #ifdef NUREX
+        #ifdef REACTIONS
         res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1;
         #endif
         res.results.push_back(r);
diff --git a/config.h b/config.h
index b78a5e4..0f584e3 100644
--- a/config.h
+++ b/config.h
@@ -60,7 +60,7 @@ namespace catima{
         //char z_effective=z_eff_type::atima14;
         unsigned char dedx = 0;
         unsigned char dedx_straggling = omega::atima;
-        #ifdef NUREX
+        #ifdef REACTIONS
         unsigned char skip=skip_none;
         #else
         unsigned char skip=skip_calculation::skip_reactions;
diff --git a/constants.h b/constants.h
index 1ac2479..fcf712c 100644
--- a/constants.h
+++ b/constants.h
@@ -20,7 +20,7 @@ constexpr double int_eps_tof = 0.01;
 */
 
 constexpr double thin_target_limit = 1 - 1e-3;
-#ifdef NUREX
+#ifdef REACTIONS
 constexpr double emin_reaction = 30.0;
 constexpr bool reactions = true;
 #else
diff --git a/reactions.cpp b/reactions.cpp
index 906b002..2499f8f 100644
--- a/reactions.cpp
+++ b/reactions.cpp
@@ -1,12 +1,18 @@
 #include "catima/reactions.h"
-
-#ifdef NUREX
-#include "nurex/Parametrization.h"
 #include "catima/catima.h"
 #include "catima/abundance_database.h"
 #include "catima/storage.h"
 #include <cmath>
 #include <iostream>
+
+#ifdef NUREX
+#include "nurex/Parametrization.h"
+using nurex::SigmaR_Kox;
+#else
+using catima::SigmaR_Kox;
+#endif
+
+
 namespace catima{
     
 double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c){
@@ -27,7 +33,7 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
         double e = energy_out(projectile.T, th, range_spline);
         for(unsigned int i = 0;i<target.ncomponents();i++){
             stn_sum += target.molar_fraction(i);
-            sum += target.molar_fraction(i)*nurex::SigmaR_Kox(ap,zp,e,at,zt); 
+            sum += target.molar_fraction(i)*SigmaR_Kox(ap,zp,e,at,zt); 
         }
         return sum/stn_sum;
     };
@@ -52,4 +58,29 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
 
 }
 
+#ifndef NUREX
+double SigmaR_Kox(int Ap, int Zp, double E, int At, int Zt){
+    constexpr double rc = 1.3;
+    constexpr double r0 = 1.1;
+    constexpr double a = 1.85;
+    constexpr double c1 = 2-(10.0/(1.5*1.5*1.5*1.5*1.5));
+    double Ap13 = pow(Ap,1.0/3.0);
+    double At13 = pow(At,1.0/3.0);
+    double D = 5.0*(At-2*Zt)*Zp/(Ap*At);
+    double Bc = Zp*Zt/(rc*(Ap13+At13));
+    double logE = std::log10(E);
+    double c = 0;
+    if(logE < 1.5){
+        c = c1*std::pow(logE/1.5,3);
+    }
+    else{
+        c = (-10.0/std::pow(logE,5)) + 2;
+    }
+    double Rs = r0 * ((a*Ap13*At13)/(Ap13+At13)-c)+D;
+    double Rv = r0 * (Ap13 + At13);
+    double Ri = Rv + Rs;
+    double Ecm = Ecm_from_T_relativistic(E,Ap,At);
+    return 10.0*PI*Ri*Ri*(1-(Bc/Ecm));
+    }
 #endif
+
diff --git a/reactions.h b/reactions.h
index 6c14b07..66b2937 100644
--- a/reactions.h
+++ b/reactions.h
@@ -50,8 +50,24 @@ namespace catima{
         return 1.0 - std::exp(-i*0.0001);
     }
 
-    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);
-    
+    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);   
+}
+#else
+
+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;
+    double plab= p_from_T(T,Ap);
+    double elab = sqrt(plab*plab + mp*mp);
+    double ecm = sqrt(mp*mp + mt*mt + 2*elab*mt);
+    double pcm = plab * mt / ecm;
+    return sqrt(pcm*pcm+mp*mp)-mp;
 }
 
 #endif //NUREX
diff --git a/structures.h b/structures.h
index b004610..ebafb17 100644
--- a/structures.h
+++ b/structures.h
@@ -227,7 +227,7 @@ namespace catima{
         double sigma_a=0.0;
         double sigma_r=0.0;
         double tof=0.0;   
-        #ifdef NUREX
+        #ifdef REACTIONS
         double sp = 1.0;
         #endif
     };
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1502833..745bd70 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,4 +1,3 @@
-enable_testing()
 add_definitions(-Dlest_FEATURE_COLOURISE=1)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/tests)
 
@@ -10,11 +9,11 @@ foreach(entry ${CATIMA_TESTS})
     add_test(${entry} ${PROJECT_BINARY_DIR}/tests/${entry})
 endforeach(entry in ${CATIMA_TESTS})
 
-if(NUREX)
+if(REACTIONS)
     add_executable(test_reaction test_reaction.cpp)
     target_link_libraries(test_reaction catima)
     add_test(test_reaction ${PROJECT_BINARY_DIR}/tests/test_reaction)
-endif(NUREX)
+endif(REACTIONS)
 
 set(C_TESTS test_c)
 MESSAGE( STATUS "CMAKE_C_COMPILER: " ${CMAKE_C_COMPILER} )
diff --git a/tests/test_reaction.cpp b/tests/test_reaction.cpp
index dc1e4e4..7ea291b 100644
--- a/tests/test_reaction.cpp
+++ b/tests/test_reaction.cpp
@@ -12,6 +12,8 @@ const lest::test specification[] =
     CASE("reaction"){
         catima::Projectile proj{12,6,6,870};
         auto c = catima::get_material(6);
+        auto h = catima::get_material(1);
+        catima::Material water({{0,8,2},{0,1,1}},1.0);
         c.thickness(2.0);
         double r;