diff --git a/catima.cpp b/catima.cpp index 053e091..2b62098 100644 --- a/catima.cpp +++ b/catima.cpp @@ -87,6 +87,23 @@ double dedx_from_range(Projectile &p, double T, const Material &t, const Config return p.A/range_spline.derivative(T); } +std::vector dedx_from_range(Projectile &p, const std::vector &T, const Material &t, const Config &c){ + auto data = _storage.Get(p,t,c); + Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); + + std::vector dedx; + dedx.reserve(T.size()); + for(auto e:T){ + if(e energy_out(Projectile &p, const std::vector &T, const Material &t, const Config &c){ + auto data = _storage.Get(p,t,c); + Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); + + std::vector eout; + eout.reserve(T.size()); + for(auto e:T){ + if(e +#include // #define NDEBUG #include "catima/build_config.h" @@ -70,6 +71,15 @@ namespace catima{ */ double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c=default_config); + /** + * returns the dEdx calculated from range spline as derivative + * @param p - Projectile + * @param T - energy vector + * @param mat - Material + * @return range + */ + std::vector dedx_from_range(Projectile &p, const std::vector &T, const Material &t, const Config &c=default_config); + /** * returns the range straggling of the Projectile in Material from spline * @param p - Projectile @@ -144,6 +154,15 @@ namespace catima{ */ double energy_out(Projectile &p, double T, const Material &t, const Config &c=default_config); + /** + * calculates outcoming energy + * @p - Projectile + * @t - Material + * @param T - incoming energy vector + * @return outcoming energy after the material in Mev/u + */ + std::vector energy_out(Projectile &p, const std::vector &T, const Material &t, const Config &c=default_config); + /** * calculates all observables for projectile passing material * @param p - Projectile diff --git a/catima.pyx b/catima.pyx index f07df73..88afac6 100644 --- a/catima.pyx +++ b/catima.pyx @@ -7,6 +7,7 @@ """ cimport catimac +from libcpp.vector cimport vector from enum import IntEnum import numpy @@ -391,13 +392,13 @@ def projectile_range(Projectile projectile, Material material, energy = None, Co def dedx_from_range(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - res[i] = catimac.dedx_from_range(projectile.cbase, e, material.cbase, config.cbase) - return res + res = catimac.dedx_from_range(projectile.cbase, energy.tolist(), material.cbase, config.cbase) + return numpy.asarray(res); + if(isinstance(energy,list)): + return catimac.dedx_from_range(projectile.cbase, energy, material.cbase, config.cbase) if(energy is None): energy = projectile.T() - return catimac.dedx_from_range(projectile.cbase, energy, material.cbase, config.cbase); + return catimac.dedx_from_range(projectile.cbase, energy, material.cbase, config.cbase); def domega2de(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): @@ -424,7 +425,7 @@ def dedx(Projectile projectile, Material material, energy = None, Config config res = numpy.empty(energy.size) for i,e in enumerate(energy): res[i] = catimac.dedx(projectile.cbase, e, material.cbase, config.cbase) - return res + return res if(energy is None): energy = projectile.T() return catimac.dedx(projectile.cbase, energy, material.cbase, config.cbase) @@ -441,13 +442,13 @@ def domega2dx(Projectile projectile, Material material, energy = None, Config co def energy_out(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - res[i] = catimac.energy_out(projectile.cbase, e, material.cbase, config.cbase) - return res + res = catimac.energy_out(projectile.cbase, energy.tolist(), material.cbase, config.cbase) + return numpy.asarray(res) + if(isinstance(energy,list)): + return catimac.energy_out(projectile.cbase, energy, material.cbase, config.cbase) if(energy is None): energy = projectile.T() - return catimac.energy_out(projectile.cbase, energy, material.cbase, config.cbase) + return catimac.energy_out(projectile.cbase, energy, material.cbase, config.cbase) def sezi_dedx_e(Projectile projectile, Target t): return catimac.sezi_dedx_e(projectile.cbase, t.cbase) diff --git a/catimac.pxd b/catimac.pxd index 149bfcf..8edc2ad 100644 --- a/catimac.pxd +++ b/catimac.pxd @@ -73,7 +73,9 @@ cdef extern from "catima/catima.h" namespace "catima": cdef double domega2dx(Projectile &p, double T, const Material &mat, const Config &c) cdef double range(Projectile &p, double T, const Material &t, const Config &c); cdef double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c); + cdef vector[double] dedx_from_range(Projectile &p, vector[double] &T, const Material &t, const Config &c); cdef double energy_out(Projectile &p, double T, const Material &t, const Config &c); + cdef vector[double] energy_out(Projectile &p, vector[double] &T, const Material &t, const Config &c); cdef double domega2de(Projectile &p, double T, const Material &t, const Config &c); cdef double da2de(Projectile &p, double T, const Material &t, const Config &c); diff --git a/tests/test.py b/tests/test.py index 9b5d317..7978ad6 100644 --- a/tests/test.py +++ b/tests/test.py @@ -144,8 +144,6 @@ class TestStructures(unittest.TestCase): def test_material_calculation(self): water = catima.get_material(catima.material.WATER) - water2 = catima.get_material(catima.material.WATER) - water2.I(78.0) p = catima.Projectile(1,1) p(1000) @@ -154,15 +152,10 @@ class TestStructures(unittest.TestCase): self.assertAlmostEqual(res.dEdxi,2.23,1) self.assertAlmostEqual(res["dEdxi"],2.23,1) self.assertAlmostEqual(res.dEdxi,res2,3) - res3 = catima.calculate(p,water2) - self.assertTrue(res.dEdxi>res3.dEdxi) - res = catima.calculate(p(500),water) res2 = catima.dedx(p,water) self.assertAlmostEqual(res.dEdxi,2.76,1) self.assertAlmostEqual(res.dEdxi,res2,3) - res3 = catima.calculate(p,water2) - self.assertTrue(res.dEdxi>res3.dEdxi) res = catima.calculate(p(9),water) res2 = catima.dedx(p,water) @@ -186,8 +179,32 @@ class TestStructures(unittest.TestCase): self.assertAlmostEqual(res.Eout,997.077,1) self.assertAlmostEqual(res["Eout"],997.077,1) self.assertAlmostEqual(res.Eout,res2,3) + + def test_eout_list(self): + graphite = catima.get_material(6) + graphite.thickness(0.5) + p = catima.Projectile(12,6) + energies = [100,500,1000] + res = catima.calculate(p(1000),graphite) + self.assertAlmostEqual(res.Eout,997.077,1) + res2 = catima.energy_out(p,graphite,energy=energies) + self.assertEqual(len(res2),len(energies)) + self.assertAlmostEqual(res2[2], 997.077,1) + self.assertAlmostEqual(res2[0], catima.calculate(p(energies[0]),graphite).Eout ,1) + self.assertAlmostEqual(res2[1], catima.calculate(p(energies[1]),graphite).Eout ,1) + + def test_dedx_from_range_list(self): + graphite = catima.get_material(6) + graphite.thickness(0.5) + p = catima.Projectile(12,6) + energies = [100,500,1000] + res2 = catima.dedx_from_range(p,graphite,energy=energies) + self.assertEqual(len(res2),len(energies)) + for i,e in enumerate(energies): + r = catima.dedx_from_range(p, graphite, energy=e) + self.assertAlmostEqual(res2[i], r, 0.1) - def test_layer_calculateion(self): + def test_layer_calculation(self): p = catima.Projectile(12,6) water = catima.get_material(catima.material.WATER) water.thickness(10.0) @@ -261,9 +278,8 @@ class TestStructures(unittest.TestCase): graphite = catima.get_material(6) graphite.thickness(1.0) data = catima.get_data(p, water) - self.assertEqual(catima.max_storage_data,50) # assuming 50, this has to be changed manually + self.assertEqual(catima.max_storage_data,100) # assuming 50, this has to be changed manually r = catima.storage_info() - print(r) #self.assertAlmostEqual(catima.da2de(p,water,et[100]),data[2][100],6) #self.assertAlmostEqual(catima.da2de(p,water,et[400]),data[2][400],6) diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 824a06d..cbaa505 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -358,6 +358,34 @@ const lest::test specification[] = c.z_effective = z_eff_type::atima14; EXPECT(z_eff_atima14(92,1900,13) == z_effective(p_u(1900.),t,c)); #endif + }, + CASE("vector_inputs"){ + catima::Projectile p{12,6,6,1000}; + catima::Material water({ + {1.00794,1,2}, + {15.9994,8,1} + }); + catima::Material graphite; + graphite.add_element(12,6,1); + graphite.density(2.0); + graphite.thickness(0.5); + + auto res = catima::calculate(p,graphite); + EXPECT( res.Eout == approx(997.07,01)); + + std::vector energies{100,500,1000}; + auto res2 = catima::energy_out(p,energies, graphite); + EXPECT(res2.size()==energies.size()); + EXPECT(res2[2] == approx(997.07,01)); + EXPECT(res2[0] == approx(catima::energy_out(p,energies[0],graphite),0.1)); + EXPECT(res2[1] == approx(catima::energy_out(p,energies[1],graphite),0.1)); + EXPECT(res2[2] == approx(catima::energy_out(p,energies[2],graphite),0.1)); + + auto res3 = catima::dedx_from_range(p,energies,graphite); + EXPECT(res3.size()==energies.size()); + EXPECT(res3[0] == approx(catima::dedx_from_range(p,energies[0],graphite),0.1)); + EXPECT(res3[1] == approx(catima::dedx_from_range(p,energies[1],graphite),0.1)); + EXPECT(res3[2] == approx(catima::dedx_from_range(p,energies[2],graphite),0.1)); } }; diff --git a/tests/test_storage.cpp b/tests/test_storage.cpp index 3d5b2ef..93b36fd 100644 --- a/tests/test_storage.cpp +++ b/tests/test_storage.cpp @@ -86,8 +86,7 @@ const lest::test specification[] = EXPECT(catima::_storage.get_index()==4); }, - CASE("test maximum storage"){ - auto maxdata = catima::max_storage_data; + CASE("test maximum storage"){ // this test assumes max storage = 50 catima::Projectile p{12,6,6,1000}; catima::Material water({ {1,1,2}, @@ -99,18 +98,18 @@ const lest::test specification[] = }); catima::_storage.Reset(); EXPECT(catima::_storage.get_index()==0); - for(int i=1;i