mirror of
https://github.com/gwm17/catima.git
synced 2024-11-26 12:08:52 -05:00
material class change
This commit is contained in:
parent
8a08fe5c0f
commit
24e53e40a3
30
catima.cpp
30
catima.cpp
|
@ -26,47 +26,41 @@ bool operator==(const Config &a, const Config&b){
|
|||
|
||||
double dedx(Projectile &p, double T, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
Target t;
|
||||
double w=0;
|
||||
if(T<=0)return 0.0;
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto res = mat.get_element(i);
|
||||
t = res.first; //struct of target
|
||||
w = res.second; //number of atoms of the element
|
||||
auto t = mat.get_element(i);
|
||||
w = mat.weight_fraction(i);
|
||||
p.T = T;
|
||||
sum += t.A*w*dedx(p,t);
|
||||
sum += w*dedx(p,t);
|
||||
}
|
||||
return sum/mat.M();
|
||||
return sum;
|
||||
}
|
||||
|
||||
double domega2dx(Projectile &p, double T, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
Target t;
|
||||
double w=0;
|
||||
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto res = mat.get_element(i);
|
||||
t = res.first; //struct of target
|
||||
w = res.second; //number of atoms of the element
|
||||
auto t= mat.get_element(i);
|
||||
w = mat.weight_fraction(i);
|
||||
p.T = T;
|
||||
sum += t.A*w*dedx_variance(p,t);
|
||||
sum += w*dedx_variance(p,t);
|
||||
}
|
||||
return sum/mat.M();
|
||||
return sum;
|
||||
}
|
||||
|
||||
double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
|
||||
double sum = 0;
|
||||
Target t;
|
||||
double w=0;
|
||||
|
||||
for(int i=0;i<mat.ncomponents();i++){
|
||||
auto res = mat.get_element(i);
|
||||
t = res.first; //struct of target
|
||||
w = res.second; //number of atoms of the element
|
||||
auto t = mat.get_element(i);
|
||||
w = mat.weight_fraction(i);
|
||||
p.T = T;
|
||||
sum += t.A*w*angular_scattering_variance(p,t);
|
||||
sum += w*angular_scattering_variance(p,t);
|
||||
}
|
||||
return sum/mat.M();
|
||||
return sum;
|
||||
}
|
||||
|
||||
|
||||
|
|
10
catima.pyx
10
catima.pyx
|
@ -19,6 +19,7 @@ cdef class Material:
|
|||
if(elements and isinstance(elements[0],list)):
|
||||
for e in elements:
|
||||
self.cbase.add_element(e[0],e[1],e[2])
|
||||
self.cbase.calculate()
|
||||
if(not thickness is None):
|
||||
self.thickness(thickness)
|
||||
if(not density is None):
|
||||
|
@ -82,6 +83,8 @@ class material(IntEnum):
|
|||
CMO2 = 218
|
||||
SUPRASIL = 219
|
||||
HAVAR = 220
|
||||
STEEL = 221
|
||||
METHANE = 222
|
||||
|
||||
def get_material(int matid):
|
||||
res = Material()
|
||||
|
@ -93,14 +96,17 @@ def get_material(int matid):
|
|||
cdef class Target:
|
||||
cdef catimac.Target cbase
|
||||
|
||||
def __cinit__(self,a,z):
|
||||
def __cinit__(self,a,z,stn):
|
||||
self.cbase.A = a
|
||||
self.cbase.Z = z
|
||||
self.cbase.stn = stn
|
||||
|
||||
def A(self):
|
||||
return self.cbase.A
|
||||
def Z(self):
|
||||
return self.cbase.Z
|
||||
def stn(self):
|
||||
return self.cbase.stn
|
||||
|
||||
cdef class Layers:
|
||||
cdef public:
|
||||
|
@ -440,7 +446,7 @@ def storage_info():
|
|||
matter = []
|
||||
for j in range(data.m.ncomponents()):
|
||||
e = data.m.get_element(j)
|
||||
matter.append([e.first.A,e.first.Z,e.second])
|
||||
matter.append([e.A,e.Z,e.stn])
|
||||
res.append({"projectile":[data.p.A,data.p.Z],"matter":matter})
|
||||
return res
|
||||
|
||||
|
|
|
@ -5,7 +5,6 @@
|
|||
:licence: GNU Affero General Public License, see LICENCE for more details
|
||||
"""
|
||||
|
||||
from libcpp.pair cimport pair
|
||||
from libcpp.vector cimport vector
|
||||
from libcpp cimport bool
|
||||
|
||||
|
@ -13,6 +12,7 @@ cdef extern from "catima/structures.h" namespace "catima":
|
|||
cdef struct Target:
|
||||
double A
|
||||
int Z
|
||||
double stn
|
||||
|
||||
cdef struct Projectile:
|
||||
double A
|
||||
|
@ -39,13 +39,14 @@ cdef extern from "catima/structures.h" namespace "catima":
|
|||
cdef cppclass Material:
|
||||
Material() except +
|
||||
void add_element(double , int , double )
|
||||
pair[Target,double] get_element(int)
|
||||
Target get_element(int)
|
||||
int ncomponents()
|
||||
double M()
|
||||
double density()
|
||||
void density(double val)
|
||||
double thickness()
|
||||
void thickness(double val)
|
||||
void calculate()
|
||||
|
||||
cdef cppclass Layers:
|
||||
Layers() except +
|
||||
|
|
|
@ -10,10 +10,12 @@ constexpr int max_datapoints = 500; // how many datapoints between logEmin and l
|
|||
constexpr int max_storage_data = 50; // number of datapoints which can be stored in cache
|
||||
|
||||
/// required integration precision (relative units)
|
||||
constexpr double int_eps_range = 0.01;
|
||||
constexpr double int_eps_range_str = 0.01;
|
||||
/*
|
||||
constexpr double int_eps_range = 0.001;
|
||||
constexpr double int_eps_range_str = 0.001;
|
||||
constexpr double int_eps_ang_str = 0.01;
|
||||
constexpr double int_eps_tof = 0.01;
|
||||
*/
|
||||
|
||||
constexpr double thin_target_limit = 1 - 1e-3;
|
||||
|
||||
|
|
|
@ -1,5 +1,6 @@
|
|||
#include "structures.h"
|
||||
#include "catima/nucdata.h"
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
namespace catima{
|
||||
|
@ -17,7 +18,7 @@ bool operator==(const Material &a, const Material&b){
|
|||
if(a.density() != b.density())return false;
|
||||
if(a.ncomponents() != b.ncomponents())return false;
|
||||
for(int i=0;i<a.ncomponents();i++){
|
||||
if(a.stn[i] != b.stn[i])return false;
|
||||
if(a.atoms[i].stn != b.atoms[i].stn)return false;
|
||||
if(a.atoms[i].A != b.atoms[i].A)return false;
|
||||
if(a.atoms[i].Z != b.atoms[i].Z)return false;
|
||||
}
|
||||
|
@ -31,9 +32,13 @@ Material::Material(const std::array<double,2> &list){
|
|||
*/
|
||||
Material::Material(std::initializer_list<std::array<double,3>>list,double _density):rho(_density){
|
||||
std::initializer_list<std::array<double,3>>::iterator it;
|
||||
atoms.reserve(list.size());
|
||||
for ( it=list.begin(); it!=list.end(); ++it){
|
||||
add_element( (*it)[0],(*it)[1],(*it)[2]);
|
||||
add_element((*it)[0],(*it)[1],(*it)[2]);
|
||||
}
|
||||
|
||||
calculate(); // calculate if needed, ie average molar mass
|
||||
|
||||
}
|
||||
|
||||
Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th){
|
||||
|
@ -42,24 +47,20 @@ Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th)
|
|||
|
||||
void Material::add_element(double _a, int _z, double _stn){
|
||||
double a = (_a>0)?_a:element_atomic_weight(_z);
|
||||
stn.push_back(_stn);
|
||||
atoms.push_back({a,_z});
|
||||
atoms.push_back({a,_z,_stn});
|
||||
molar_mass += _stn*a;
|
||||
}
|
||||
|
||||
|
||||
std::pair<Target,double> Material::get_element(int i) const{
|
||||
return std::pair<Target,double>(atoms[i],stn[i]);
|
||||
}
|
||||
|
||||
/*
|
||||
void Material::calculate(){
|
||||
molar_mass = 0;
|
||||
for(int i=0;i<ncomponents();i++){
|
||||
molar_mass += stn[i]*a[i];
|
||||
if(std::all_of(atoms.begin(),atoms.end(),[](const Target &t){return t.stn<1.0;})){
|
||||
double sum = 0;
|
||||
for(auto& e: atoms){
|
||||
sum+= e.stn/e.A;
|
||||
}
|
||||
molar_mass = 1.0/sum;
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
Layers& Layers::operator=(const Layers& other){
|
||||
|
||||
|
|
17
structures.h
17
structures.h
|
@ -51,6 +51,7 @@ namespace catima{
|
|||
struct Target{
|
||||
double A=0;
|
||||
int Z=0;
|
||||
double stn=1.0;
|
||||
};
|
||||
|
||||
/**
|
||||
|
@ -65,7 +66,6 @@ namespace catima{
|
|||
double th=0;
|
||||
double molar_mass=0;
|
||||
std::vector<Target>atoms;
|
||||
std::vector<double>stn;
|
||||
|
||||
public:
|
||||
Material(){};
|
||||
|
@ -91,6 +91,9 @@ namespace catima{
|
|||
Material(std::initializer_list<std::array<double,3>>list,double _density=0.0);
|
||||
//Material(const std::array<double,2> &list);
|
||||
|
||||
/**
|
||||
* calculates internal variables if needed
|
||||
*/
|
||||
void calculate();
|
||||
|
||||
/**
|
||||
|
@ -104,14 +107,20 @@ namespace catima{
|
|||
/**
|
||||
* returns i-th element of the Material as a std::pair of Target and corresponding stoichiometric number
|
||||
* @param i - index of element to return
|
||||
* @return std::pair of Target and corresponding stoichiometric number
|
||||
* @return Target class
|
||||
*/
|
||||
std::pair<Target,double> get_element(int i) const;
|
||||
Target get_element(int i) const {return atoms[i];};
|
||||
|
||||
/**
|
||||
* return weight fraction of i-th element
|
||||
* @return weight fraction
|
||||
*/
|
||||
double weight_fraction(int i) const {return (atoms[i].stn<1.0)?atoms[i].stn:atoms[i].stn*atoms[i].A/M();};
|
||||
|
||||
/**
|
||||
* @return number of components in Material
|
||||
*/
|
||||
int ncomponents() const {return stn.size();}
|
||||
int ncomponents() const {return atoms.size();}
|
||||
|
||||
/**
|
||||
* @return returns Molar Mass of the Material
|
||||
|
|
|
@ -1,9 +1,9 @@
|
|||
#include "lest.hpp"
|
||||
#include "testutils.h"
|
||||
#include <math.h>
|
||||
using namespace std;
|
||||
|
||||
using catima::approx;
|
||||
#include "catima/catima.h"
|
||||
using lest::approx;
|
||||
bool rcompare(double a, double b,double eps){
|
||||
if(fabs((a-b)/fabs(b))<eps){
|
||||
return true;
|
||||
|
@ -333,9 +333,13 @@ const lest::test specification[] =
|
|||
auto air = catima::get_material(catima::material::Air);
|
||||
air.thickness(0.500);
|
||||
auto res = catima::calculate(p(350),air);
|
||||
EXPECT(res.Eout == approx(345.6).epsilon(1e-2));
|
||||
EXPECT(res.sigma_a == approx(0.0013).epsilon(1e-2));
|
||||
EXPECT(res.sigma_E == approx(0.12).epsilon(1e-2));
|
||||
EXPECT(res.Eout == approx(345.6).epsilon(1.0));
|
||||
EXPECT(res.sigma_a == approx(0.0013).epsilon(1e-4));
|
||||
EXPECT(res.sigma_E == approx(0.12).epsilon(1e-3));
|
||||
|
||||
auto water = catima::get_material(catima::material::Water);
|
||||
auto res2 = catima::calculate(p(600),water,600);
|
||||
EXPECT(res2.dEdxi == approx(92.5).epsilon(2));
|
||||
}
|
||||
|
||||
};
|
||||
|
|
|
@ -44,8 +44,8 @@ double rdif(double v1, double v2){
|
|||
void comp_dedx(catima::Projectile p, catima::Material t, double epsilon = 0.001, bool fout=false, const char* fname=""){
|
||||
double dif;
|
||||
//struct atima_results results;
|
||||
int tz = t.get_element(0).first.Z;
|
||||
double ta = t.get_element(0).first.A;
|
||||
int tz = t.get_element(0).Z;
|
||||
double ta = t.get_element(0).A;
|
||||
bool res;
|
||||
cout<<"-----------------------------"<<endl;
|
||||
cout<<"projectile: A = "<<p.A<<", Z = "<<p.Z<<endl;
|
||||
|
|
|
@ -1,6 +1,8 @@
|
|||
#include "lest.hpp"
|
||||
#include "testutils.h"
|
||||
#include <math.h>
|
||||
using namespace std;
|
||||
using catima::approx;
|
||||
|
||||
#include "catima/catima.h"
|
||||
#include "catima/material_database.h"
|
||||
|
@ -57,28 +59,25 @@ const lest::test specification[] =
|
|||
CASE("Material automatic atomic weight"){
|
||||
catima::Material water({{0,1,2},{0,8,1}});
|
||||
catima::Material graphite(0,6);
|
||||
EXPECT(water.get_element(0).first.A == 1.00794);
|
||||
EXPECT(water.get_element(1).first.A == 15.9994);
|
||||
EXPECT(graphite.get_element(0).first.A == 12.0107);
|
||||
EXPECT(water.get_element(0).A == 1.00794);
|
||||
EXPECT(water.get_element(1).A == 15.9994);
|
||||
EXPECT(graphite.get_element(0).A == 12.0107);
|
||||
EXPECT(water.M()>16);
|
||||
EXPECT(graphite.M()>12);
|
||||
},
|
||||
CASE("default materials"){
|
||||
catima::Material m = catima::get_material(6);
|
||||
EXPECT(m.get_element(0).first.A == 12.0107);
|
||||
EXPECT(m.get_element(0).first.Z == 6);
|
||||
EXPECT(m.get_element(0).A == 12.0107);
|
||||
EXPECT(m.get_element(0).Z == 6);
|
||||
EXPECT(m.density() == 2.0);
|
||||
EXPECT(m.M() == 12.0107);
|
||||
|
||||
m = catima::get_material(catima::material::Water);
|
||||
EXPECT(m.get_element(0).first.A == 1.00794);
|
||||
EXPECT(m.get_element(0).first.Z == 1);
|
||||
EXPECT(m.get_element(1).first.A == 15.9994);
|
||||
EXPECT(m.get_element(1).first.Z == 8);
|
||||
EXPECT(m.get_element(0).A == 1.00794);
|
||||
EXPECT(m.get_element(0).Z == 1);
|
||||
EXPECT(m.get_element(1).A == 15.9994);
|
||||
EXPECT(m.get_element(1).Z == 8);
|
||||
EXPECT(m.density() == 1.0);
|
||||
EXPECT(m.M() > 16.0);
|
||||
|
||||
|
||||
},
|
||||
CASE("Layers"){
|
||||
catima::Material water2;
|
||||
|
@ -102,10 +101,10 @@ const lest::test specification[] =
|
|||
EXPECT(detector1[0].thickness()==1.0);
|
||||
EXPECT(detector1[1].density()==1.0);
|
||||
EXPECT(detector1[1].thickness()==2.0);
|
||||
EXPECT(detector1[0].get_element(0).first.Z == 6);
|
||||
EXPECT(detector1[0].get_element(0).first.A == 12);
|
||||
EXPECT(detector1[1].get_element(0).first.A == 1);
|
||||
EXPECT(detector1[1].get_element(0).first.Z == 1);
|
||||
EXPECT(detector1[0].get_element(0).Z == 6);
|
||||
EXPECT(detector1[0].get_element(0).A == 12);
|
||||
EXPECT(detector1[1].get_element(0).A == 1);
|
||||
EXPECT(detector1[1].get_element(0).Z == 1);
|
||||
EXPECT(detector1[2].density() == 1.8);
|
||||
detector1[1].density(1.2);
|
||||
EXPECT(detector1[1].density()==1.2);
|
||||
|
@ -163,43 +162,62 @@ const lest::test specification[] =
|
|||
});
|
||||
EXPECT(mat2.ncomponents()==1);
|
||||
EXPECT(mat3.ncomponents()==1);
|
||||
EXPECT(mat3.get_element(0).first.A==12.01);
|
||||
EXPECT(mat3.get_element(0).A==12.01);
|
||||
EXPECT(mat4.ncomponents()==1);
|
||||
EXPECT(mat4.get_element(0).first.A==12.01);
|
||||
EXPECT(mat4.get_element(0).A==12.01);
|
||||
EXPECT(mat5.ncomponents()==2);
|
||||
EXPECT(mat5.get_element(0).first.A==12.01);
|
||||
EXPECT(mat5.get_element(0).first.Z==6);
|
||||
EXPECT(mat5.get_element(1).first.A==16.0);
|
||||
EXPECT(mat5.get_element(1).second==2);
|
||||
EXPECT(mat5.get_element(0).A==12.01);
|
||||
EXPECT(mat5.get_element(0).Z==6);
|
||||
EXPECT(mat5.get_element(1).A==16.0);
|
||||
EXPECT(mat5.get_element(1).stn==2);
|
||||
|
||||
catima::Material mat6;
|
||||
mat6 = mat5;
|
||||
EXPECT(mat5==mat6);
|
||||
EXPECT(mat5.ncomponents()==mat6.ncomponents());
|
||||
EXPECT(mat5.get_element(0).first.A==mat6.get_element(0).first.A);
|
||||
EXPECT(mat5.get_element(1).first.A==mat6.get_element(1).first.A);
|
||||
EXPECT(mat5.get_element(0).A==mat6.get_element(0).A);
|
||||
EXPECT(mat5.get_element(1).A==mat6.get_element(1).A);
|
||||
mat5.add_element(12,6,1);
|
||||
EXPECT(mat5.ncomponents()==mat6.ncomponents()+1);
|
||||
|
||||
},
|
||||
CASE("int and double stn check"){
|
||||
catima::Projectile p{1,1,1,1000};
|
||||
catima::Material mat1({
|
||||
{12.01, 6, 1},
|
||||
{16.00, 8, 1}
|
||||
CASE("fraction vs stn init"){
|
||||
catima::Projectile p{12,6};
|
||||
catima::Material water1({
|
||||
{0, 1, 0.111894},
|
||||
{0, 8, 0.888106}
|
||||
});
|
||||
catima::Material mat2({
|
||||
{12.01, 6, 0.5},
|
||||
{16.00, 8, 0.5}
|
||||
catima::Material water2({
|
||||
{0, 1, 2},
|
||||
{0, 8, 1}
|
||||
});
|
||||
mat1.thickness(1.0);
|
||||
mat2.thickness(1.0);
|
||||
auto res1 = catima::calculate(p(1000),mat1);
|
||||
auto res2 = catima::calculate(p(1000),mat2);
|
||||
EXPECT(res1.dEdxi == res2.dEdxi);
|
||||
EXPECT(res1.range == res2.range);
|
||||
EXPECT(res1.sigma_a == res2.sigma_a);
|
||||
EXPECT(res1.sigma_r == res2.sigma_r);
|
||||
water1.thickness(1.0);
|
||||
water2.thickness(1.0);
|
||||
auto res1 = catima::calculate(p(600),water1);
|
||||
auto res2 = catima::calculate(p(600),water2);
|
||||
EXPECT(res1.dEdxi == approx(res2.dEdxi,0.001));
|
||||
EXPECT(res1.range == approx(res2.range).R(1e-2));
|
||||
EXPECT(res1.sigma_a == approx(res2.sigma_a).R(1e-2));
|
||||
EXPECT(res1.sigma_r == approx(res2.sigma_r).R(1e-2));
|
||||
},
|
||||
CASE("fraction calculation"){
|
||||
catima::Material water1({
|
||||
{0, 1, 0.111894},
|
||||
{0, 8, 0.888106}
|
||||
});
|
||||
catima::Material water2({
|
||||
{0, 1, 2},
|
||||
{0, 8, 1}
|
||||
});
|
||||
|
||||
EXPECT(water1.weight_fraction(0)==0.111894);
|
||||
EXPECT(water2.weight_fraction(0)==approx(water1.weight_fraction(0)).R(1e-3));
|
||||
EXPECT(water1.weight_fraction(1)==0.888106);
|
||||
EXPECT(water2.weight_fraction(1)==approx(water1.weight_fraction(1)).R(1e-3));
|
||||
EXPECT(water2.M()==approx(18).epsilon(0.1));
|
||||
|
||||
EXPECT(water1.M()==approx(6.0,0.1));
|
||||
EXPECT(water2.M()==approx(18,0.1));
|
||||
}
|
||||
};
|
||||
|
||||
|
|
71
tests/testutils.h
Normal file
71
tests/testutils.h
Normal file
|
@ -0,0 +1,71 @@
|
|||
namespace catima{
|
||||
|
||||
class approx
|
||||
{
|
||||
public:
|
||||
explicit approx ( double magnitude, double eps=1.0)
|
||||
: epsilon_ { eps }
|
||||
, magnitude_{ magnitude } {}
|
||||
|
||||
approx( approx const & other ) = default;
|
||||
|
||||
approx operator()( double magnitude, double eps = 1.0 )
|
||||
{
|
||||
approx approx ( magnitude, eps);
|
||||
return approx;
|
||||
}
|
||||
|
||||
double magnitude() const { return magnitude_; }
|
||||
|
||||
approx & epsilon( double epsilon ) { epsilon_ = epsilon; return *this; }
|
||||
approx & R( double relative ) { epsilon_ = relative*magnitude_; return *this; }
|
||||
|
||||
|
||||
friend bool operator == ( double lhs, approx const & rhs )
|
||||
{
|
||||
return std::abs( lhs - rhs.magnitude_) < rhs.epsilon_;
|
||||
}
|
||||
|
||||
friend bool operator == ( approx const & lhs, double rhs ) { return operator==( rhs, lhs ); }
|
||||
friend bool operator != ( double lhs, approx const & rhs ) { return !operator==( lhs, rhs ); }
|
||||
friend bool operator != ( approx const & lhs, double rhs ) { return !operator==( rhs, lhs ); }
|
||||
|
||||
friend bool operator <= ( double lhs, approx const & rhs ) { return lhs < rhs.magnitude_ || lhs == rhs; }
|
||||
friend bool operator <= ( approx const & lhs, double rhs ) { return lhs.magnitude_ < rhs || lhs == rhs; }
|
||||
friend bool operator >= ( double lhs, approx const & rhs ) { return lhs > rhs.magnitude_ || lhs == rhs; }
|
||||
friend bool operator >= ( approx const & lhs, double rhs ) { return lhs.magnitude_ > rhs || lhs == rhs; }
|
||||
|
||||
//private:
|
||||
double epsilon_;
|
||||
double magnitude_;
|
||||
};
|
||||
|
||||
std::ostream & operator<<(std::ostream &os, approx const &a){
|
||||
using lest::to_string;
|
||||
return os<<to_string(a.magnitude_)<<" +- "<<a.epsilon_;
|
||||
}
|
||||
|
||||
|
||||
bool rdiff(double a, double b,double eps){
|
||||
if(fabs((a-b)/fabs(b))<eps){
|
||||
return true;
|
||||
}
|
||||
else{
|
||||
std::cout<<"\033[1;31m"<<a<<" == "<<b<<"\033[0m"<<std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
bool diff(double a, double b,double eps){
|
||||
if(fabs((a-b))<eps){
|
||||
return true;
|
||||
}
|
||||
else{
|
||||
std::cout<<"\033[1;31m"<<a<<" == "<<b<<"\033[0m"<<std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
Loading…
Reference in New Issue
Block a user