Merge pull request #801 from tardis-sn/decay

Decay of isotopes branch

docs/examples/abundancecust.rst

@@ -2,8 +2,8 @@
 Using a custom stratified composition
 *************************************
 
-Overview
-========
+ASCII Format
+=============
 
 To use a stratified ejecta composition in TARDIS, the elemental abundances may
 be specified on a per-cell basis via an external ascii file (similar to setting
@@ -46,12 +46,76 @@ The example file shown here has three simple layers:
    including details of how to develop your own dataset to suit your
    needs, please contact us.
 
-TARDIS input file
-=================
+TARDIS ascii input file
+=======================
 
 If you create a correctly formatted abundance profile file (called "abund.dat"
 in this example), you can use it in TARDIS by putting the following lines in
 the model section of the yaml file:
 
 .. literalinclude:: tardis_configv1_abundance_cust_example.yml
     :language: yaml
+
+
+CSV Format
+==========
+
+In this format, both elemental and isotopic abundances may
+be specified on a per-cell basis via an external csv file. A csv file that could
+work on a mesh with ten cells should be formatted like this:
+
+.. csv-table:: Example
+    :file: tardis_model_abund.csv
+    :delim: space 
+    :header-rows: 1
+
+In this file:
+
+- Header row contains element and isotopes symbol 
+- the remaining entries in each row give the set of elemental and isotopic abundances.
+
+The abundances are specified as mass fractions (i.e. the sum of columns
+in each row should be 1.0). The mass fractions specified will be adopted directly in
+the TARDIS calculations.
+
+The example file shown here has three simple layers: 
+
+- an innermost region that is composed of Si and two Nickel Isotopes Ni56 and Ni58
+
+- a middle region that is composed of O and Mg
+
+- an outer region that is composed of C and O.
+
+.. note::
+
+    Suppose you specify Elemental and Isotopic abundances for the same element. For ex-
+    :code:`Ni` and :code:`Ni56`. 
+    Here, Ni will refer to the stable Nickel, i.e. (Z=26, A=58).
+
+.. warning::
+
+   The example given here is to show the format only. It is not a
+   realistic model. In any real calculation better resolution
+   (i.e. more grid points) should be used.
+
+TARDIS csv input file
+=====================
+
+If you create a correctly formatted abundance profile file (called "tardis_model_abund.csv"
+in this example), you can use it in TARDIS by putting the following lines in
+the model section of the yaml file:
+
+.. literalinclude:: tardis_configv1_isotope_abundance_cust_example.yml
+    :language: yaml   
+
+Convert ascii abundance file format to csv format
+=================================================
+
+If you want to convert an ASCII abundance file(say "abund.dat") to CSV format, you can use 
+:code:`convert_abundances_format` function for it. Here is an example to demonstrate this:
+
+.. code:: python
+
+    from tardis.util import convert_abundances_format
+    df = convert_abundances_format('abund.dat')
+    df.to_csv('converted_abund.csv', index=False, sep=' ')  

docs/examples/abundanceuni.rst

@@ -8,16 +8,28 @@ Overview
 The simplest possibility for specifying an ejecta composition is to use a
 uniform abundance pattern in all cells specified in the density profile setup
 step. These constant abundances  can be supplied directly in the input (yaml)
-file. Elemental abundances are set in the "abundances" subsection of the "model"
+file. Elemental and Isotopic abundances are set in the "abundances" subsection of the "model"
 section, following the "type: uniform" specifier (see example input
 file below). They are specified as mass fractions. E.g.
 
 .. code-block:: none
 
     Si: 0.6
-    S: 0.4
+    S: 0.2
+    Ni56: 0.2
 
-will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
+will set the mass fraction of silicon (Z=14) to 0.6, sulphur (Z=16) to 0.2 and Nickel(Z=26, A=56) to 0.2.
+
+.. note::
+
+    Suppose you specify Elemental and Isotopic abundance for the same element. For ex-
+
+    .. code-block:: none
+
+        Ni: 0.8
+        Ni56: 0.2
+
+    Here, Ni will refer to the stable Nickel, i.e. (Z=26, A=58).
 
 .. note::
 

docs/examples/tardis_configv1_abundance_uniform_example.yml

@@ -5,5 +5,7 @@ model:
         Mg: 0.03
         Si: 0.52
         S: 0.19
-        Ar: 0.04
+        Ar: 0.02
         Ca: 0.03
+        Ni56: 0.01
+        Ni58: 0.01

docs/examples/tardis_configv1_isotope_abundance_cust_example.yml

@@ -0,0 +1,5 @@
+model:
+    abundances:
+        type: file
+        filename: tardis_model_abund.csv
+        filetype: tardis_model

docs/examples/tardis_model_abund.csv

@@ -0,0 +1,11 @@
+C O Mg Si Ni56 Ni58
+0 0 0 0.6 0.4 0
+0 0 0 0.1 0.5 0.4
+0 0 0 0.3 0 0.7
+0 0.2 0.8 0 0 0
+0 0.3 0.7 0 0 0
+0 0.2 0.8 0 0 0
+0 0.2 0.8 0 0 0
+0 0.2 0.8 0 0 0
+0.5 0.5 0 0 0 0
+0.5 0.5 0 0 0 0

docs/scripts.rst

@@ -0,0 +1,19 @@
+***************
+Running scripts
+***************
+
+Convert CMFGEN files to TARDIS file format
+==========================================
+
+This script takes a CMFGEN file as an input , and an output path to save converted files in new TARDIS file format. 
+CSV file contains abundances of both elements and isotopes.
+DAT file contains values of velocity, density, electron_density and temperature.  
+
+Format of command - :code:`cmfgen2tardis /path/to/input_file path/to/output/`  
+
+Example, for how to use-  
+
+
+.. code-block:: bash
+
+   $ cmfgen2tardis tardis_example/DDC15_SN_HYDRO_DATA_0.976d tardis_example/

setup.py

@@ -107,7 +107,8 @@
     entry_points[hook_ep] = ['%s = %s' % (hook_name, hook_func)]
 
 entry_points['console_scripts'] = [
-    'tardis_test_runner = tardis.tests.integration_tests.runner:run_tests'
+    'tardis_test_runner = tardis.tests.integration_tests.runner:run_tests',
+    'cmfgen2tardis = tardis.scripts.cmfgen2tardis:main'
 ]
 
 #from Cython.Build import cythonize

tardis/io/decay.py

@@ -0,0 +1,123 @@
+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 = [dict() for i in xrange(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(data=0.0, 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 cls(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 = [dict() for i in xrange(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]
+        for i in range(len(materials)):
+            materials[i].update(decayed_materials[i])
+        df = IsotopeAbundances.from_materials(materials)
+        df.sort_index(inplace=True)
+        return df 
+
+    def as_atoms(self):
+        """
+        Merge Isotope dataframe according to atomic number 
+
+        Returns:
+            : merged isotope abundances
+        """
+
+        return self.groupby('atomic_number').sum()
+
+    def merge(self, other, normalize=True):
+        """
+        Merge Isotope dataframe with abundance passed as parameter 
+
+        Parameters
+        ----------
+        other: pd.DataFrame 
+        normalize : bool
+            If true, resultant dataframe will be normalized
+
+        Returns:
+            : merged abundances
+        """
+        isotope_abundance = self.as_atoms()
+        isotope_abundance = isotope_abundance.fillna(0.0)
+        #Merge abundance and isotope dataframe
+        modified_df = isotope_abundance.add(other, fill_value=0)
+
+        if normalize:
+            norm_factor = modified_df.sum(axis=0)
+            modified_df /= norm_factor
+
+        return modified_df