Merge pull request #4 from hrosiak/python

Python

catima.pyx

@@ -1,16 +1,41 @@
+"""
+    catima python module
+    ~~~~~~~~~~~
+    This module provides interface to the catima c++ library
+    :copyright: (c) 2017 by Andrej Prochazka
+    :licence: GNU Affero General Public License, see LICENCE for more details
+"""
+
 cimport catimac
 from enum import IntEnum
 import numpy
 
 cdef class Material:
     cdef catimac.Material cbase
-    def __cinit__(self, elements):
+    def __cinit__(self, elements=None, thickness=None, density=None):
         self.cbase = catimac.Material()
-        if(elements and isinstance(elements[0],int)):
+        if(elements and (isinstance(elements[0],float) or isinstance(elements[0],int))):
             self.cbase.add_element(elements[0],elements[1],elements[2])
         if(elements and isinstance(elements[0],list)):
             for e in elements:
                 self.cbase.add_element(e[0],e[1],e[2])
+        if(not thickness is None):
+            self.thickness(thickness)
+        if(not density is None):
+            self.density(density)
+
+    cdef from_c(self, catimac.Material &other):
+        self.cbase = other
+    
+    cdef catimac.Material getc(self):
+        cdef catimac.Material res
+        res = self.cbase 
+        return res
+    
+    def copy(self):
+        res = Material()
+        res.cbase = self.cbase
+        return res
 
     def add_element(self, a, z , s):
         self.cbase.add_element(a, z, s)
@@ -36,6 +61,35 @@ cdef class Material:
         else:
             return self.cbase.thickness(val)
 
+class material(IntEnum):
+            PLASTIC = 201
+            AIR = 202
+            CH2 = 203
+            LH2 = 204
+            LD2 = 205
+            WATER = 206
+            DIAMOND = 207
+            GLASS = 208
+            ALMG3 = 209
+            ARCO2_30 = 210
+            CF4 = 211
+            ISOBUTANE = 212
+            KAPTON = 213
+            MYLAR = 214
+            NAF = 215
+            P10 = 216
+            POLYOLEFIN = 217
+            CMO2 = 218
+            SUPRASIL = 219
+            HAVAR = 220
+
+def get_material(int matid):
+    res = Material()
+    cdef catimac.Material cres = catimac.get_material(matid);
+    res.from_c(cres)
+    return res
+
+
 cdef class Target:
     cdef catimac.Target cbase
 
@@ -49,34 +103,56 @@ cdef class Target:
         return self.cbase.Z
 
 cdef class Layers:
-    cdef catimac.Layers cbase
+    cdef public:
-    
+        materials
-    def __cinit__(self):
+    def __init__(self):
-        self.cbase = catimac.Layers()
+        self.materials=[]
-        self.materials = []
         
     def add(self,Material m):
-        self.cbase.add(m.cbase)
+        self.materials.append(m.copy())
-        self.materials.append(m)
         
     def num(self):
-        return self.cbase.num()
+        return len(self.materials)
+
+    def get(self, key):
+        return self.materials[key]
         
     def __getitem__(self, key):
-        if(isinstance(key,int)):
+        if(isinstance(key,int) and key<self.num()):
-            return self.materials[key]
+            return self.get(key)
+        return None
+    
+    def __add__(self, other):
+        res = Layers()
+        for e in self.materials:
+            res.add(e)
+            
+        if(isinstance(other,Layers)):
+            for e in other.materials:
+                res.add(e)
+        if(isinstance(other,Material)):
+            res.add(other.copy())
+        return res
+    
+    cdef catimac.Layers getc(self):
+        cdef catimac.Layers res
+        #for l in self.materials:
+        #    res.add(l.getc())
+        return res
 
 cdef class Projectile:
     cdef catimac.Projectile cbase
-    def __cinit__(self, a, z, t=None,q=None):
+    def __cinit__(self, A, Z, Q=None,T=None):
-        self.cbase.A = a
+        self.cbase.A = A
-        self.cbase.Z = z
+        self.cbase.Z = Z
-        self.cbase.Q = z
+        self.cbase.Q = Z
-        if(q):
+        if(Q):
-            self.cbase.Q = q
+            self.cbase.Q = Q
-        if(t):
+        if(T):
-            self.cbase.T = t
+            self.cbase.T = T
-    def T(self,val):
+    def T(self,val=None):
+        if(val is None):
+            return self.cbase.T
         self.cbase.T = val;
     def __call__(self,val=None):
         if(val is None):
@@ -88,6 +164,8 @@ cdef class Projectile:
         return self.cbase.A
     def Z(self):
         return self.cbase.Z
+    def Q(self):
+        return self.cbase.Q
 
 cdef class Result:
     cdef public double Ein
@@ -112,6 +190,24 @@ cdef class Result:
         self.sigma_r=0.0
         self.tof=0.0
 
+    def get_dict(self):
+        return {"Ein":self.Ein,
+                "Eout":self.Eout,
+                "Eloss":self.Eloss,
+                "range":self.range,
+                "dEdxi":self.dEdxi,
+                "dEdxo":self.dEdxo,
+                "sigma_E":self.sigma_E,
+                "sigma_a":self.sigma_a,
+                "sigma_r":self.sigma_r,
+                "tof":self.tof,
+                }
+
+    def __getitem__(self,key):
+        d = self.get_dict()
+        if(key in d):
+            return d[key]
+
     cdef setc(self,catimac.Result &val):
         self.Ein=val.Ein
         self.Eout=val.Eout
@@ -124,6 +220,29 @@ cdef class Result:
         self.sigma_r=val.sigma_r
         self.tof=val.tof
 
+cdef class MultiResult:
+    cdef public Result total_result
+    cdef public results
+    cdef public total
+
+    def __init__(self):
+        self.total_result = Result()
+        self.results = []
+        self.total = {}
+
+    cdef setc(self, catimac.MultiResult &val):
+        self.total_result.setc(val.total_result)
+        for e in val.results:
+            self.results.append(e)
+        self.total = self.total_result.get_dict()
+
+    def __getitem__(self,key):
+        if(isinstance(key,int) and key<len(self.results)):
+            return self.results[key]
+        if(isinstance(key,str) and key in self.total):
+            return self.total[key]
+        return None
+
 class z_eff_type(IntEnum):
     none = 0,
     atima = 1
@@ -176,6 +295,30 @@ def calculate(Projectile projectile, Material material, energy = None, Config co
     res.setc(cres)
     return res
 
+cdef catimac.Layers get_clayers(Layers layers):
+    cdef catimac.Layers res
+    cdef catimac.Material m
+    for l in layers.materials:
+        m = get_cmaterial(l)
+        res.add(m)
+    return res
+
+cdef catimac.Material get_cmaterial(Material material):
+    cdef catimac.Material res
+    res = material.cbase
+    return res
+
+def calculate_layers(Projectile projectile, Layers layers, energy = None,  Config config = default_config):
+    cdef catimac.Layers clayers
+    clayers = catimac.Layers()
+    clayers = get_clayers(layers)
+    if(not energy is None):
+        projectile.T(energy)
+    cdef catimac.MultiResult cres = catimac.calculate(projectile.cbase, clayers, config.cbase)
+    res = MultiResult()
+    res.setc(cres)
+    return res
+
 def range(Projectile projectile, Material material, energy = None, Config config = default_config):  
     if(isinstance(energy,numpy.ndarray)):
         res = numpy.empty(energy.size)

catimac.pxd

@@ -1,5 +1,13 @@
+"""
+    catima cython
+    ~~~~~~~~~~~~~~~~~
+    :copyright: (c) 2017 by Andrej Prochazka
+    :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
 
 cdef extern from "catima/structures.h" namespace "catima":
     cdef struct Target:
@@ -24,12 +32,14 @@ cdef extern from "catima/structures.h" namespace "catima":
         double sigma_r
         double tof
 
-    cdef cppclass MultiResult
+    cdef cppclass MultiResult:
+        vector[Result] results
+        Result total_result
 
     cdef cppclass Material:
         Material() except +
         void add_element(double , int , double )
-        pair[Target,double] getElement(int) 
+        pair[Target,double] get_element(int) 
         int ncomponents()
         double M()
         double density()
@@ -45,6 +55,9 @@ cdef extern from "catima/structures.h" namespace "catima":
         Material& operator[](int i)
         Layers& operator=(const Layers& other)
 
+cdef extern from "catima/material_database.h" namespace "catima":
+    cdef Material get_material(int)
+
 cdef extern from "catima/config.h" namespace "catima":
     cdef struct Config:
         char z_effective;

tests/test.py

@@ -0,0 +1,193 @@
+import sys
+sys.path.insert(0,"../build")
+import unittest
+import catima
+
+class TestStructures(unittest.TestCase):
+
+    def test_Projectile(self):
+        p = catima.Projectile(238,92)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),92)
+
+        p = catima.Projectile(238,92,90)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),90)
+        p.T(1000)
+        self.assertEqual(p.T(),1000)
+        self.assertEqual(p(),1000)
+        p(500)
+        self.assertEqual(p.T(),500)
+        self.assertEqual(p(),500)
+
+        p = catima.Projectile(238,92,90, T=100)
+        self.assertEqual(p.T(),100)
+    
+    def test_Material(self):
+        mat = catima.Material()
+        mat.add_element(12,6,1)
+        self.assertEqual(mat.ncomponents(),1)
+        mat.add_element(1,1,2)
+        self.assertEqual(mat.ncomponents(),2)
+
+        mat2 = catima.Material([12.01,6,1])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertEqual(mat2.molar_mass(),12.01)
+
+        mat3 = catima.Material([12,6,1])
+        self.assertEqual(mat3.ncomponents(),1)
+        self.assertEqual(mat3.molar_mass(),12)
+
+        water = catima.Material([[1,1,2],[16,8,1]])
+        self.assertEqual(water.molar_mass(),18)
+
+        mat2 = catima.Material([0,6,1])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertAlmostEqual(mat2.molar_mass(),12,1)
+
+        mat5 = catima.Material([0,6,1],density=1.9, thickness=0.5)
+        self.assertEqual(mat5.ncomponents(),1)
+        self.assertEqual(mat5.thickness(),0.5)
+        self.assertEqual(mat5.density(),1.9)
+
+        # copy
+        mat3.density(1.8)
+        matc = mat3.copy()
+        self.assertEqual(matc.ncomponents(),1)
+        self.assertEqual(matc.molar_mass(),12)
+        self.assertEqual(matc.density(),1.8)
+        mat3.density(2.0)
+        self.assertEqual(matc.density(),1.8)
+        self.assertEqual(mat3.density(),2.0)
+
+    
+    def test_default_material(self):
+        m1 = catima.get_material(6);
+        self.assertAlmostEqual(m1.molar_mass(),12,1)
+        self.assertEqual(m1.ncomponents(),1)
+        self.assertAlmostEqual(m1.density(),2.0,1)
+
+        m2 = catima.get_material(catima.material.WATER)
+        self.assertEqual(m2.ncomponents(),2)
+        self.assertAlmostEqual(m2.molar_mass(),18,1)
+        self.assertAlmostEqual(m2.density(),1.0,1)
+
+    def test_layers(self):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p10 = catima.get_material(catima.material.P10)
+        p10.thickness(0.01)
+        n2 = catima.get_material(7)
+        n2.thickness(0.02)
+
+        mat= catima.Layers()
+        self.assertEqual(mat.num(),0)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),1)
+        self.assertAlmostEqual(mat[0].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        self.assertAlmostEqual(mat[0].density(),2.0,1)
+        
+        mat.add(p10)
+        self.assertEqual(mat.num(),2)
+
+        graphite.thickness(1.0)
+        graphite.density(1.8)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),3)
+        self.assertAlmostEqual(mat[2].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        self.assertAlmostEqual(mat[0].density(),2.0,1)
+        self.assertAlmostEqual(mat[2].thickness(),1.0,1)
+        self.assertAlmostEqual(mat[2].density(),1.8,1)
+        mat[2].thickness(1.2)
+        mat[2].density(1.9)
+        self.assertAlmostEqual(mat.materials[2].thickness(),1.2,1)
+        self.assertAlmostEqual(mat.materials[2].density(),1.9,1)
+        #self.assertAlmostEqual(mat.materials[0].thickness(),0.5,1)
+        #self.assertAlmostEqual(mat.materials[0].density(),2.0,1)
+        self.assertEqual(mat[3],None)
+        self.assertEqual(mat["a"],None)
+
+        mat2 = catima.Layers()
+        mat2.add(n2)
+        self.assertEqual(mat2.num(),1)
+
+        mats = mat2 + mat
+        self.assertEqual(mats.num(),4)
+        self.assertAlmostEqual(mats[0].molar_mass(),14,1)
+        self.assertEqual(mats[0].thickness(),0.02)
+        self.assertAlmostEqual(mats[1].molar_mass(),12,1)
+        self.assertAlmostEqual(mats[3].molar_mass(),12,1)
+
+        n2.thickness(0.5)
+        mats = mats + n2
+        self.assertEqual(mats.num(),5)
+        self.assertAlmostEqual(mats[0].molar_mass(),14,1)
+        self.assertEqual(mats[0].thickness(),0.02)
+        self.assertAlmostEqual(mats[4].molar_mass(),14,1)
+        self.assertEqual(mats[4].thickness(),0.5)
+
+    def test_material_calculation(self):
+        water = catima.get_material(catima.material.WATER)
+        p = catima.Projectile(1,1)
+        
+        p(1000)
+        res = catima.calculate(p,water)
+        res2 = catima.dedx(p,water)
+        self.assertAlmostEqual(res.dEdxi,2.23,1)
+        self.assertAlmostEqual(res["dEdxi"],2.23,1)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+        res = catima.calculate(p(500),water)
+        res2 = catima.dedx(p,water)
+        self.assertAlmostEqual(res.dEdxi,2.76,1)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+
+        res = catima.calculate(p(9),water)
+        res2 = catima.dedx(p,water)
+        self.assertAlmostEqual(res.dEdxi,51.17,1)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+        res = catima.calculate(p(9),water)
+        res = catima.calculate(p(9),water)
+        self.assertAlmostEqual(res.dEdxi,51.17,1)
+    
+    def test_eout(self):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p = catima.Projectile(12,6)
+        res = catima.calculate(p(1000),graphite)
+        res2 = catima.energy_out(p(1000),graphite)
+        self.assertAlmostEqual(res.Eout,997.077,1)
+        self.assertAlmostEqual(res["Eout"],997.077,1)
+        self.assertAlmostEqual(res.Eout,res2,3)
+
+    def test_layer_calculateion(self):
+        p = catima.Projectile(12,6)
+        water = catima.get_material(catima.material.WATER)
+        water.thickness(10.0)
+        graphite = catima.get_material(6)
+        graphite.thickness(1.0)
+
+        mat = catima.Layers()
+        mat.add(water)
+        mat.add(graphite)
+        res = catima.calculate_layers(p(1000),mat)
+        self.assertEqual(len(res.results),2)
+        self.assertAlmostEqual(res.total_result.Eout,926.3,1)
+        self.assertAlmostEqual(res.total_result.sigma_a,0.00269,1)
+        self.assertAlmostEqual(res["Eout"],926.3,1)
+        self.assertAlmostEqual(res["sigma_a"],0.00269,4)
+        self.assertAlmostEqual(res["tof"],0.402,2)
+        self.assertAlmostEqual(res["Eloss"],884,0)
+
+        self.assertAlmostEqual(res[0]["Eout"],932,0)
+        self.assertAlmostEqual(res[1]["Eout"],926,0)
+        self.assertAlmostEqual(res[0]["sigma_a"],0.00258,4)
+        self.assertAlmostEqual(res[1]["sigma_a"],0.000774,4)
+        self.assertAlmostEqual(res[0]["range"],107.1,0)
+        self.assertAlmostEqual(res[1]["range"],110.7,0)
+
+if __name__ == "__main__":
+    unittest.main()

tests/test_structures.cpp

@@ -49,7 +49,7 @@ const lest::test specification[] =
                 EXPECT(water2.M()==18);
             }
             SECTION("equal operator check"){
-              EXPECT(water==water);
+              EXPECT(water==water2);
               EXPECT(!(water==graphite));
             }
         }