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

dx partial

Browse Source
This commit is contained in:
hrocho 2020-08-04 01:44:00 +02:00
parent 413741bf0d
commit 8ed9c469e5
6 changed files with 112 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
catima.cpp
View File

@ -261,7 +261,8 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
double t0 = std::min(range, t.thickness());
return (tt-t0)*(tt-t0)*angular_variance_spline.derivative(x);
};
res.sigma_x = sqrt(integrator.integrate(fx2,res.Eout, res.Ein));
res.sigma_x = integrator.integrate(fx2,res.Eout, res.Ein)/t.density()/t.density();
res.sigma_x = sqrt(res.sigma_x);
#ifdef REACTIONS
res.sp = nonreaction_rate(p,t,c);
#endif
@ -273,8 +274,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
double e = p.T;
res.total_result.Ein = e;
res.results.reserve(layers.num());
double z = 0;
for(auto&m:layers.get_materials()){
z += m.thickness_cm();
Result r = calculate(p,m,e,c);
e = r.Eout;
res.total_result.sigma_a += r.sigma_a*r.sigma_a;
@ -282,7 +284,9 @@ 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;
res.total_result.sigma_x += r.sigma_x*r.sigma_x;
double temp = r.sigma_a*(layers.thickness_cm() - z);
res.total_result.sigma_x += (r.sigma_x*r.sigma_x) + (temp*temp);
//res.total_result.sigma_x += r.sigma_x*r.sigma_x;
#ifdef REACTIONS
res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1;
#endif
@ -297,6 +301,7 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
else{
res.total_result.sigma_a = 0.0;
res.total_result.sigma_E = 0.0;
res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
}
return res;
}

2
pymodule/pycatima.cpp
View File

@ -139,7 +139,7 @@ PYBIND11_MODULE(pycatima,m){
.def("molar_mass",py::overload_cast<>(&Material::M, py::const_), "get mass")
.def("thickness",py::overload_cast<>(&Material::thickness, py::const_), "get thickness")
.def("thickness",py::overload_cast<double>(&Material::thickness), "set thickness")
.def("thickness_cm",&Material::thickness_cm,"set thickness in cm unit")
.def("thickness_cm",py::overload_cast<double>(&Material::thickness_cm),"set thickness in cm unit")
.def("I",py::overload_cast<>(&Material::I, py::const_), "get I")
.def("I",py::overload_cast<double>(&Material::I), "set I")
.def("__str__",&material_to_string);

16
structures.cpp
View File

@ -66,6 +66,22 @@ void Layers::add(Material m){
materials.push_back(m);
}
double Layers::thickness() const {
double sum = 0;
for(auto &m : materials){
sum += m.thickness();
}
return sum;
}
double Layers::thickness_cm() const {
double sum = 0;
for(auto &m : materials){
sum += m.thickness_cm();
}
return sum;
}
Layers operator+(const Layers &a, const Layers&b){
Layers res;
for(auto &e:a.materials){

10
structures.h
View File

@ -281,6 +281,16 @@ namespace catima{
*/
int num()const{return materials.size();};
/**
* @ return total thickness
*/
double thickness() const;
/**
* @ return total thickness in cm
*/
double thickness_cm() const;
Material& operator[](int i){return materials[i];}
friend Layers operator+(const Layers &a, const Layers&b);

69
tests/test_calculations.cpp
View File

@ -138,6 +138,7 @@ using namespace std;
catima::Material graphite({
{12.011,6,1},
});
graphite.density(2);
auto res = catima::calculate(p(1000),graphite);
CHECK(catima::dedx(p(1000), graphite) == approx(res.dEdxi).R(0.001) );
@ -183,6 +184,74 @@ using namespace std;
dif = res.tof - 0.038;
CHECK( fabs(dif) < 0.01);
}
TEST_CASE("angular scattering"){
catima::Projectile p{1,1,1,158.6};
catima::Material cu = catima::get_material(29);
cu.thickness_cm(0.02963);
auto res = catima::calculate(p,cu);
CHECK( 1000*res.sigma_a == approx(7.2,0.5));
cu.thickness_cm(0.2963);
res = catima::calculate(p,cu);
CHECK( 1000*res.sigma_a == approx(23.5,3));
catima::Layers ll;
cu.thickness_cm(0.2963/2.);
ll.add(cu);
ll.add(cu);
auto res2 = catima::calculate(p,ll);
CHECK( 1000*res2.total_result.sigma_a == approx(23.5,3));
CHECK( 1000*res2.total_result.sigma_a == approx(1000*res.sigma_a,0.05));
}
TEST_CASE("displacement test"){
catima::Projectile p{1,1,1,215};
catima::Material water({
{1.00794,1,2},
{15.9994,8,1}
});
water.density(1.0);
catima::Material water2({
{1.00794,1,2},
{15.9994,8,1}
});
water2.density(2.0);
water.thickness_cm(9.6);
auto res = catima::calculate(p,water);
CHECK( res.sigma_x == approx(0.1,0.03));
auto resb = catima::calculate(p,water2);
CHECK( res.sigma_x > resb.sigma_x);
water.thickness_cm(17.0);
auto res17 = catima::calculate(p,water);
CHECK( res17.sigma_x == approx(0.25,0.05));
water.thickness_cm(29.4);
auto res29 = catima::calculate(p,water);
CHECK( res29.sigma_x == approx(0.66,0.08));
catima::Layers ll;
water.thickness_cm(9.6);
ll.add(water);
auto res2 = catima::calculate(p,ll);
CHECK(res2.total_result.sigma_x == approx(0.1,0.03));
water.thickness_cm(7.4);
ll.add(water);
res2 = catima::calculate(p,ll);
CHECK(res2.total_result.sigma_x == approx(0.25,0.05));
CHECK(res2.total_result.sigma_x == approx(res17.sigma_x,0.001));
catima::Layers l;
water.thickness_cm(9.6/3);
l.add(water);
l.add(water);
l.add(water);
res2 = catima::calculate(p,l);
CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01));
}
TEST_CASE("result from stopped material"){
catima::Projectile p{12,6,6,1000};
catima::Material water({

8
tests/test_structures.cpp
View File

@ -42,8 +42,15 @@ using namespace std;
CHECK(water.thickness()==approx(1.0,0.0001));
water.thickness_cm(1.0);
CHECK(water.thickness()==approx(1.0,0.0001));
CHECK(water.thickness_cm()==approx(1.0,0.0001));
water.thickness_cm(2.0);
CHECK(water.thickness()==approx(2.0,0.0001));
CHECK(water.thickness_cm()==approx(2.0,0.0001));
water.density(2.0);
CHECK(water.thickness()==approx(2.0,0.0001));
CHECK(water.thickness_cm()==approx(1.0,0.0001));
water.thickness(1.0);
CHECK(water.thickness_cm()==approx(0.5,0.0001));
}
TEST_CASE("Material automatic atomic weight"){
catima::Material water({{0,1,2},{0,8,1}});
@ -85,6 +92,7 @@ using namespace std;
detector1.add(water2);
detector1.add(graphite);
CHECK(detector1.num() == 3);
CHECK(detector1.thickness() == approx(4.0,1e-6));
// check correct density and thickness
CHECK(detector1[0].density()==1.8);
CHECK(detector1[0].thickness()==1.0);