From afd735d11d4791118184021a84b346671ab2d532 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@hrocho.localdomain>
Date: Thu, 6 Aug 2020 01:28:19 +0200
Subject: [PATCH] dxrc2

---
 catima.cpp                  | 33 +++++++++------------------------
 integrator.h                |  2 +-
 tests/test_calculations.cpp | 23 +++++++++++++++--------
 3 files changed, 25 insertions(+), 33 deletions(-)

diff --git a/catima.cpp b/catima.cpp
index be00a11..de2b670 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -269,40 +269,25 @@ 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;
-
-    auto fx2 = [&](double x)->double{ 
-	    double range = range_spline(T);
-        double tt = range - range_spline(x);      
-        double t0 = std::min(range, t.thickness());
-	    return (tt-t0)*(tt-t0)*da2dx(p(x),t, c);
-            };
+    
+    double rrange = std::min(res.range/t.density(), t.thickness_cm());
     auto fx2p = [&](double x)->double{ 
 	    double range = range_spline(T);
         double e =energy_out(T,x*t.density(),range_spline);
-        double t0 = std::min(range/t.density(), t.thickness_cm());
-	    return (t0-x)*(t0-x)*da2dx(p(e), t, c)*t.density();
+	    return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
             };           
     
-    //res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density();    
-    res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,2);    
+    //res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3, 1e-3,1);    
+    res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange);    
     res.sigma_x = sqrt(res.sigma_x);
 
-    auto fx1 = [&](double x)->double{ 
-	    double range = range_spline(T);
-        double tt = range - range_spline(x);      
-        double t0 = std::min(range, t.thickness());
-	    return (t0-tt)*da2dx(p(x), t, c);
-            };
     auto fx1p = [&](double x)->double{ 
-	    double range = range_spline(T);
-        double e =energy_out(T,x*t.density(),range_spline);
-        double t0 = std::min(range/t.density(), t.thickness_cm());
-	    return (t0-x)*da2dx(p(e), t, c)*t.density();
+        double e =energy_out(T,x*t.density(),range_spline);        
+	    return (rrange-x)*da2dx(p(e), t, c)*t.density();
             };
     
-    //res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4);
-    res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,2);
-    p.T = T;
+    //res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,1);
+    res.cov = integrator.integrate(fx1p,0, rrange);
     #ifdef REACTIONS
     res.sp = nonreaction_rate(p,t,c);
     #endif
diff --git a/integrator.h b/integrator.h
index a039055..e9225d8 100644
--- a/integrator.h
+++ b/integrator.h
@@ -62,7 +62,7 @@ using integrator_type = IntegratorGSL;
 #else
 using integrator_type = GaussLegendreIntegration<8>;
 #endif
-using integrator_adaptive_type = GaussKronrodIntegration<21>;
+using integrator_adaptive_type = GaussKronrodIntegration<15>;
 
 extern integrator_type integrator;
 extern integrator_adaptive_type integrator_adaptive;
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index e34c767..15881df 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -254,14 +254,21 @@ using namespace std;
       res2 = catima::calculate(p(215),l);
       CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.03));
 
-      catima::Layers lll;
-      water.thickness_cm(29.4/4);
-      lll.add(water);
-      lll.add(water);
-      lll.add(water);
-      lll.add(water);
-      res2 = catima::calculate(p(215),lll);
-      CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
+      for(int n=2;n<10;n++){
+        catima::Layers lll;
+        water.thickness_cm(29.4/n);
+        for(int i=0;i<n;i++)lll.add(water);        
+        res2 = catima::calculate(p(215),lll);
+        CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
+        }
+      
+
+      /// position sigma is taken from range thickness if particle stoppped inside
+      water.thickness_cm(30);
+      auto r1 = catima::calculate(p,water);
+      water.thickness_cm(40);
+      auto r2 = catima::calculate(p,water);
+      CHECK(r1.sigma_x == approx(r2.sigma_x).R(0.01));
 
     }
     TEST_CASE("result from stopped material"){