gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-22 18:28:51 -05:00
Code Issues Projects Releases Wiki Activity

dxrc2

Browse Source
This commit is contained in:
hrocho 2020-08-06 01:28:19 +02:00
parent e2c3db8405
commit afd735d11d
3 changed files with 25 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
catima.cpp
View File

@ -270,39 +270,24 @@ 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();
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

2
integrator.h
View File

@ -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;

17
tests/test_calculations.cpp
View File

@ -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));
for(int n=2;n<10;n++){
catima::Layers lll;
water.thickness_cm(29.4/4);
lll.add(water);
lll.add(water);
lll.add(water);
lll.add(water);
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"){