From b38fe06d1fffcd7b3ce7f1209803c41cd8dbbbe2 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@hrocho.localdomain>
Date: Wed, 5 Aug 2020 01:10:03 +0200
Subject: [PATCH] backup2

---
 calculations.cpp            |  2 +-
 catima.cpp                  | 30 ++++++++++++++++++++++--------
 tests/test_calculations.cpp |  5 +++--
 3 files changed, 26 insertions(+), 11 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index bb71d0e..8915122 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -486,7 +486,7 @@ double angular_scattering_variance(Projectile &p, Target &t){
     return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
 }
 
-/// radioation lengths are taken frm Particle Data Group 2014
+/// radiation lengths are taken from Particle Data Group 2014
 double radiation_length(int z, double m){
     double lr = 0;
     if(z==1){return 63.04;}
diff --git a/catima.cpp b/catima.cpp
index 1044188..9dc285c 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -259,23 +259,35 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
 	    double range = range_spline(T);
         double tt = range - range_spline(x);      
         double t0 = std::min(range, t.thickness());
-	    return (tt-t0)*(tt-t0)*angular_variance_spline.derivative(x);
+	    return (tt-t0)*(tt-t0)*da2dx(p,x, t, c);
             };
-           
+    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();
+            };           
     
-    res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density()/t.density();    
-    //res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4)/t.density()/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,4);    
     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)*angular_variance_spline.derivative(x);
+	    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();
             };
     
-    res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4)/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-6, 1e-3,4);
+    p.T = T;
     #ifdef REACTIONS
     res.sp = nonreaction_rate(p,t,c);
     #endif
@@ -298,7 +310,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
         res.total_result.tof += r.tof;
         res.total_result.Eout = r.Eout;        
         double a2 = res.total_result.sigma_a;
-        res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov) + (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x;
+        res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov) 
+                                 + (a2*m.thickness_cm()*m.thickness_cm()) 
+                                 + r.sigma_x*r.sigma_x;
         //res.total_result.sigma_x += (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x;
         res.total_result.cov += a2*m.thickness_cm() + r.cov;
         res.total_result.sigma_a += r.sigma_a*r.sigma_a;
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index 60cb286..d921583 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -217,6 +217,7 @@ using namespace std;
                 {1.00794,1,2},
                 {15.9994,8,1}
                 });
+      water2.thickness_cm(9.6);
       water2.density(2.0);
 
       water.thickness_cm(9.6);
@@ -250,7 +251,7 @@ using namespace std;
       l.add(water);
       l.add(water);
       l.add(water);
-      res2 = catima::calculate(p,l);
+      res2 = catima::calculate(p(215),l);
       CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01));
 
       catima::Layers lll;
@@ -259,7 +260,7 @@ using namespace std;
       lll.add(water);
       lll.add(water);
       lll.add(water);
-      res2 = catima::calculate(p,lll);
+      res2 = catima::calculate(p(215),lll);
       CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
 
     }