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

Merge pull request #15 from hrosiak/material

Browse Source 
material class change
This commit is contained in:
Andrej Prochazka 2017-12-14 15:36:28 +01:00 committed by GitHub
parent 8a08fe5c0f 24e53e40a3
commit aff32c24bf
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 198 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
catima.cpp
View File

@ -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
View File

@ -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
catimac.pxd
View File

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

8
constants.h
View File

@ -2,7 +2,7 @@
#define CONSTANTS_H
namespace catima {
constexpr double Ezero = 9E-4; // lowest E to calculate, below taken as 0
constexpr double logEmin = -3; // log of minimum energy
constexpr double logEmax = 5.0; // log of max energy
@ -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;

4
storage.cpp
View File

@ -37,8 +37,8 @@ void Data::Add(const Projectile &p, const Material &t, Config c){
#if(0)
*index = dp;
index->range = calculate_range(p,t,c);
index->range_straggling = calculate_range_straggling(p,t,c);
index->angular_variance = calculate_angular_variance(p,t,c);
index->range_straggling = calculate_range_straggling(p,t,c);
index->angular_variance = calculate_angular_variance(p,t,c);
#else
*index = calculate_DataPoint(p,t,c);
#endif

29
structures.cpp
View File

@ -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
View File

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

14
tests/test_calculations.cpp
View File

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

4
tests/test_dedx_range.cpp
View File

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

98
tests/test_structures.cpp
View File

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