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