From 24e53e40a34d3c4b012e96c735e1d69e4d61217c Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Thu, 14 Dec 2017 15:29:23 +0100
Subject: [PATCH] material class change

---
 catima.cpp                  | 30 +++++-------
 catima.pyx                  | 10 +++-
 catimac.pxd                 |  5 +-
 constants.h                 |  8 +--
 storage.cpp                 |  4 +-
 structures.cpp              | 29 +++++------
 structures.h                | 17 +++++--
 tests/test_calculations.cpp | 14 ++++--
 tests/test_dedx_range.cpp   |  4 +-
 tests/test_structures.cpp   | 98 ++++++++++++++++++++++---------------
 tests/testutils.h           | 71 +++++++++++++++++++++++++++
 11 files changed, 198 insertions(+), 92 deletions(-)
 create mode 100644 tests/testutils.h

diff --git a/catima.cpp b/catima.cpp
index a971619..95ce76d 100644
--- a/catima.cpp
+++ b/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;
 }
 
 
diff --git a/catima.pyx b/catima.pyx
index be4b1b0..62368c9 100644
--- a/catima.pyx
+++ b/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
 
diff --git a/catimac.pxd b/catimac.pxd
index bbe1424..20b9569 100644
--- a/catimac.pxd
+++ b/catimac.pxd
@@ -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 +
diff --git a/constants.h b/constants.h
index b2bedb7..4353c67 100644
--- a/constants.h
+++ b/constants.h
@@ -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;
 
diff --git a/storage.cpp b/storage.cpp
index 07de1d7..3ab54c6 100644
--- a/storage.cpp
+++ b/storage.cpp
@@ -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
diff --git a/structures.cpp b/structures.cpp
index ce8b7cb..d013243 100644
--- a/structures.cpp
+++ b/structures.cpp
@@ -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){
     
diff --git a/structures.h b/structures.h
index 2f03cdb..bf11a6d 100644
--- a/structures.h
+++ b/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
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index 868834a..51ca74a 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -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));
     }
     
 };
diff --git a/tests/test_dedx_range.cpp b/tests/test_dedx_range.cpp
index 816ac91..d6d6df4 100644
--- a/tests/test_dedx_range.cpp
+++ b/tests/test_dedx_range.cpp
@@ -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;
diff --git a/tests/test_structures.cpp b/tests/test_structures.cpp
index 0d9ea69..42f8fb3 100644
--- a/tests/test_structures.cpp
+++ b/tests/test_structures.cpp
@@ -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));
     }
 };
 
diff --git a/tests/testutils.h b/tests/testutils.h
new file mode 100644
index 0000000..975fbcb
--- /dev/null
+++ b/tests/testutils.h
@@ -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;
+    }
+      
+}
+
+}