Merge pull request #756 from vg3095/radioactive_decay

[TEP007] [WIP] Isotope Abundances class

tardis/atomic.py

@@ -493,7 +493,7 @@ def prepare_atom_data(self, selected_atomic_numbers, line_interaction_type='scat
             self.lines = self.lines[self.lines['ion_number'] <= max_ion_number]
 
         # self.lines.sort(['wavelength', 'line_id'], inplace=True)
-        self.lines.sort(['wavelength'], inplace=True)
+        self.lines.sort_values(by='wavelength', inplace=True)
         self.lines.set_index('line_id', inplace=True)
 
 
@@ -672,4 +672,3 @@ def get_collision_matrix(self, species, t_electrons):
                                            t_electrons.reshape((1, 1, t_electrons.shape[0]))) * \
                       self.g_ratio_matrices[species].reshape((no_of_levels, no_of_levels, 1))
         return c_ul_matrix + c_lu_matrix.transpose(1, 0, 2)
-

tardis/io/decay.py

@@ -0,0 +1,89 @@
+import pandas as pd
+from pyne import nucname, material
+from astropy import units as u
+
+class IsotopeAbundances(pd.DataFrame):
+
+    @property
+    def _constructor(self):
+        return IsotopeAbundances
+
+    def _update_material(self):
+        self.comp_dicts = [{}] * len(self.columns)
+        for (atomic_number, mass_number), abundances in self.iterrows():
+            nuclear_symbol = '%s%d'.format(nucname.name(atomic_number),
+                                           mass_number)
+            for i in xrange(len(self.columns)):
+                self.comp_dicts[i][nuclear_symbol] = abundances[i]
+
+    @classmethod
+    def from_materials(cls, materials):
+        multi_index_tuples = set([])
+        for material in materials:
+            multi_index_tuples.update([cls.id_to_tuple(key)
+                                       for key in material.keys()])
+
+        index = pd.MultiIndex.from_tuples(
+            multi_index_tuples, names=['atomic_number', 'mass_number'])
+
+
+        abundances = pd.DataFrame(index=index, columns=xrange(len(materials)))
+
+        for i, material in enumerate(materials):
+            for key, value in material.items():
+                abundances.loc[cls.id_to_tuple(key), i] = value
+
+        return abundances
+
+
+
+
+    @staticmethod
+    def id_to_tuple(atomic_id):
+        return nucname.znum(atomic_id), nucname.anum(atomic_id)
+
+
+    def to_materials(self):
+        """
+        Convert DataFrame to a list of materials interpreting the MultiIndex as
+        atomic_number and mass_number
+
+        Returns
+        -------
+            : ~list
+            list of pyne Materialss
+        :return:
+        """
+
+        comp_dicts = [{}] * len(self.columns)
+        for (atomic_number, mass_number), abundances in self.iterrows():
+            nuclear_symbol = '{0:s}{1:d}'.format(nucname.name(atomic_number),
+                                           mass_number)
+            for i in xrange(len(self.columns)):
+                comp_dicts[i][nuclear_symbol] = abundances[i]
+        return [material.Material(comp_dict) for comp_dict in comp_dicts]
+
+
+
+    def decay(self, t):
+        """
+        Decay the Model
+
+        Parameters
+        ----------
+
+        t: ~float or ~astropy.units.Quantity
+            if float it will be understood as days
+
+        Returns:
+            : decayed abundances
+        """
+
+        materials = self.to_materials()
+        t_second = u.Quantity(t, u.day).to(u.s).value
+
+        decayed_materials = [item.decay(t_second) for item in materials]
+
+        df = IsotopeAbundances.from_materials(decayed_materials)
+        df.sort_index(inplace=True)
+        return df 

tardis/io/tests/test_decay.py

@@ -0,0 +1,20 @@
+import pytest
+import pandas as pd
+
+from tardis.io.decay import IsotopeAbundances
+from numpy.testing import assert_almost_equal
+
+@pytest.fixture
+def simple_abundance_model():
+    index = pd.MultiIndex.from_tuples([(28, 56)],
+                                      names=['atomic_number', 'mass_number'])
+    return IsotopeAbundances([[1.0, 1.0]], index=index)
+
+def test_simple_decay(simple_abundance_model):
+    decayed_abundance = simple_abundance_model.decay(100)
+    assert_almost_equal(decayed_abundance.ix[26, 56][0], 0.55752)
+    assert_almost_equal(decayed_abundance.ix[26, 56][1], 0.55752)
+    assert_almost_equal(decayed_abundance.ix[27, 56][0], 0.4423791)
+    assert_almost_equal(decayed_abundance.ix[27, 56][1], 0.4423791)
+    assert_almost_equal(decayed_abundance.ix[28, 56][0], 1.1086e-05)
+    assert_almost_equal(decayed_abundance.ix[28, 56][1], 1.1086e-05)

tardis/montecarlo/tests/test_formal_integral.py

@@ -49,7 +49,7 @@ def test_trapezoid_integration(clib, N):
     func = clib.trapezoid_integration
     func.restype = c_double
     h = 1.
-
+    N = int(N)
     func.argtypes = [
             ndpointer(c_double),
             c_double,

tardis_env27.yml

@@ -1,20 +1,25 @@
 name: tardis
+
+channels:
+    - defaults
+    - conda-forge
+
 dependencies:
 - python=2.7
-- numpy=1.10
-- scipy=0.17.0
-- pandas=0.16.2
-- pytables=3.2.2
-- h5py=2.5
-- matplotlib=1.4.3
-- astropy=1.1.2
-- numexpr=2.4.4
+- numpy=1.12
+- scipy=0.18
+- pandas=0.20
+- pytables
+- h5py=2.6
+- matplotlib=2.0
+- astropy=1.3
+- numexpr=2.6
 - Cython=0.21
 - networkx=1.10
 - pytest=3.0
 - pyyaml=3.12
 - jsonschema=2.5.1
-
+- pyne=0.5.3
 
 # RTD requirements
 - sphinx=1.5.1