diff --git a/docs/examples/abundancecust.rst b/docs/examples/abundancecust.rst index 19475608657..19e340647de 100644 --- a/docs/examples/abundancecust.rst +++ b/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,8 +46,8 @@ 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 @@ -55,3 +55,67 @@ 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=' ') \ No newline at end of file diff --git a/docs/examples/abundanceuni.rst b/docs/examples/abundanceuni.rst index 7df5f52d834..07017fb9648 100644 --- a/docs/examples/abundanceuni.rst +++ b/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:: diff --git a/docs/examples/tardis_configv1_abundance_uniform_example.yml b/docs/examples/tardis_configv1_abundance_uniform_example.yml index f354293b387..77472625871 100644 --- a/docs/examples/tardis_configv1_abundance_uniform_example.yml +++ b/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 diff --git a/docs/examples/tardis_configv1_isotope_abundance_cust_example.yml b/docs/examples/tardis_configv1_isotope_abundance_cust_example.yml new file mode 100644 index 00000000000..740d4148ae4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/docs/examples/tardis_model_abund.csv b/docs/examples/tardis_model_abund.csv new file mode 100644 index 00000000000..46200ee5b77 --- /dev/null +++ b/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 diff --git a/docs/scripts.rst b/docs/scripts.rst new file mode 100644 index 00000000000..1d62e45ca63 --- /dev/null +++ b/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/ diff --git a/setup.py b/setup.py index 24ace215f7a..ba47ae50217 100755 --- a/setup.py +++ b/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 diff --git a/tardis/io/decay.py b/tardis/io/decay.py new file mode 100644 index 00000000000..d37a92808c0 --- /dev/null +++ b/tardis/io/decay.py @@ -0,0 +1,122 @@ +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 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 = [{}] * 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 + + 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() + #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 diff --git a/tardis/io/model_reader.py b/tardis/io/model_reader.py index 7ada3550da6..4c95ae842b6 100644 --- a/tardis/io/model_reader.py +++ b/tardis/io/model_reader.py @@ -4,6 +4,7 @@ from numpy import recfromtxt, genfromtxt import pandas as pd from astropy import units as u +from pyne import nucname import logging # Adding logging support @@ -11,7 +12,6 @@ from tardis.util import parse_quantity - class ConfigurationError(Exception): pass @@ -42,12 +42,21 @@ def read_density_file(filename, filetype): """ file_parsers = {'artis': read_artis_density, - 'simple_ascii': read_simple_ascii_density} + 'simple_ascii': read_simple_ascii_density, + 'tardis_model': read_cmfgen_density} + + electron_densities = None + temperature = None + if filetype == 'tardis_model': + (time_of_model, velocity, + unscaled_mean_densities, electron_densities, temperature) = read_cmfgen_density(filename) + else: + (time_of_model, velocity, + unscaled_mean_densities) = file_parsers[filetype](filename) - (time_of_model, velocity, - unscaled_mean_densities) = file_parsers[filetype](filename) v_inner = velocity[:-1] v_outer = velocity[1:] + invalid_volume_mask = (v_outer - v_inner) <= 0 if invalid_volume_mask.sum() > 0: message = "\n".join(["cell {0:d}: v_inner {1:s}, v_outer " @@ -59,7 +68,7 @@ def read_density_file(filename, filetype): raise ConfigurationError("Invalid volume of following cell(s):\n" "{:s}".format(message)) - return time_of_model, velocity, unscaled_mean_densities + return time_of_model, velocity, unscaled_mean_densities, electron_densities, temperature def read_abundances_file(abundance_filename, abundance_filetype, inner_boundary_index=None, outer_boundary_index=None): @@ -85,9 +94,17 @@ def read_abundances_file(abundance_filename, abundance_filetype, """ file_parsers = {'simple_ascii': read_simple_ascii_abundances, - 'artis': read_simple_ascii_abundances} + 'artis': read_simple_ascii_abundances, + 'tardis_model': read_simple_isotope_abundances} + + isotope_abundance = pd.DataFrame() + if abundance_filetype == 'tardis_model': + index, abundances, isotope_abundance = read_simple_isotope_abundances( + abundance_filename) + else: + index, abundances = file_parsers[abundance_filetype]( + abundance_filename) - index, abundances = file_parsers[abundance_filetype](abundance_filename) if outer_boundary_index is not None: outer_boundary_index_m1 = outer_boundary_index - 1 else: @@ -95,9 +112,53 @@ def read_abundances_file(abundance_filename, abundance_filetype, index = index[inner_boundary_index:outer_boundary_index] abundances = abundances.ix[:, slice(inner_boundary_index, outer_boundary_index_m1)] abundances.columns = np.arange(len(abundances.columns)) - return index, abundances + return index, abundances, isotope_abundance +def read_uniform_abundances(abundances_section, no_of_shells): + """ + Parameters + ---------- + + abundances_section: ~config.model.abundances + no_of_shells: int + + Returns + ------- + abundance: ~pandas.DataFrame + isotope_abundance: ~pandas.DataFrame + """ + abundance = pd.DataFrame(columns=np.arange(no_of_shells), + index=pd.Index(np.arange(1, 120), + name='atomic_number'), + dtype=np.float64) + + isotope_index = pd.MultiIndex( + [[]] * 2, [[]] * 2, names=['atomic_number', 'mass_number']) + isotope_abundance = pd.DataFrame(columns=np.arange(no_of_shells), + index=isotope_index, + dtype=np.float64) + + for element_symbol_string in abundances_section: + if element_symbol_string == 'type': + continue + try: + if element_symbol_string in nucname.name_zz: + z = nucname.name_zz[element_symbol_string] + abundance.ix[z] = float( + abundances_section[element_symbol_string]) + else: + mass_no = nucname.anum(element_symbol_string) + z = nucname.znum(element_symbol_string) + isotope_abundance.loc[(z, mass_no), :] = float( + abundances_section[element_symbol_string]) + + except RuntimeError as err: + raise RuntimeError( + "Abundances are not defined properly in config file : {}".format(err.args)) + + return abundance, isotope_abundance + def read_simple_ascii_density(fname): """ Reading a density file of the following structure (example; lines starting with a hash will be ignored): @@ -183,6 +244,54 @@ def read_artis_density(fname): return time_of_model, velocity, mean_density +def read_cmfgen_density(fname): + """ + Reading a density file of the following structure (example; lines starting with a hash will be ignored): + The first density describes the mean density in the center of the model and is not used. + The file consists of a header row and next row contains unit of the respective attributes + velocity densities electron_densities temperature + km/s g/cm^3 /cm^3 K + 871.66905 4.2537191e-09 2.5953807e+14 7.6395577 + 877.44269 4.2537191e-09 2.5953807e+14 7.6395577 + + Rest columns contain abundances of elements and isotopes + + Parameters + ---------- + + fname: str + filename or path with filename + + + Returns + ------- + + time_of_model: ~astropy.units.Quantity + time at which the model is valid + + velocity: ~np.ndarray + mean_density: ~np.ndarray + electron_densities: ~np.ndarray + temperature: ~np.ndarray + + """ + df = pd.read_csv(fname, comment='#', delimiter='\s+', skiprows=[0, 2]) + + with open(fname) as fh: + for row_index, line in enumerate(fh): + if row_index == 0: + time_of_model_string = line.strip().replace('t0:', '') + time_of_model = parse_quantity(time_of_model_string) + elif row_index == 2: + quantities = line.split() + + velocity = u.Quantity(df['velocity'].values, quantities[0]).to('cm/s') + temperature = u.Quantity(df['temperature'].values, quantities[1])[1:] + mean_density = u.Quantity(df['densities'].values, quantities[2])[1:] + electron_densities = u.Quantity( + df['electron_densities'].values, quantities[3])[1:] + + return time_of_model, velocity, mean_density, electron_densities, temperature def read_simple_ascii_abundances(fname): """ @@ -212,3 +321,58 @@ def read_simple_ascii_abundances(fname): abundances = pd.DataFrame(data[1:,1:].transpose(), index=np.arange(1, data.shape[1])) return index, abundances + + +def read_simple_isotope_abundances(fname, delimiter='\s+'): + """ + Reading an abundance file of the following structure (example; lines starting with hash will be ignored): + The first line of abundances describe the abundances in the center of the model and are not used. + First 4 columns contain values related to velocity, density, electron_density and temperature. + From 5th column onwards, abundances of elements and isotopes begin. + The file consists of a header row and next row contains unit of the respective attributes + Since abundance fractions are unitless , its unit row is filled with ones + Example + velocity...temperature C O Ni56 + km/s.........K 1 1 1 + ...................... 0.4 0.3 0.2 + + Parameters + ---------- + + fname: str + filename or path with filename + + Returns + ------- + + index: ~np.ndarray + abundances: ~pandas.DataFrame + isotope_abundance: ~pandas.MultiIndex + """ + df = pd.read_csv(fname, comment='#', + delimiter=delimiter, skiprows=[0, 2]) + df = df.transpose() + + abundance = pd.DataFrame(columns=np.arange(df.shape[1] - 1), + index=pd.Index([], + name='atomic_number'), + dtype=np.float64) + + isotope_index = pd.MultiIndex( + [[]] * 2, [[]] * 2, names=['atomic_number', 'mass_number']) + isotope_abundance = pd.DataFrame(columns=np.arange(df.shape[1] - 1), + index=isotope_index, + dtype=np.float64) + + #First 4 columns related to density parser (e.g. velocity) + for element_symbol_string in df.index[4:]: + if element_symbol_string in nucname.name_zz: + z = nucname.name_zz[element_symbol_string] + abundance.loc[z, :] = df.loc[element_symbol_string].tolist()[1:] + else: + z = nucname.znum(element_symbol_string) + mass_no = nucname.anum(element_symbol_string) + isotope_abundance.loc[( + z, mass_no), :] = df.loc[element_symbol_string].tolist()[1:] + + return abundance.index, abundance, isotope_abundance diff --git a/tardis/io/schemas/model.yml b/tardis/io/schemas/model.yml index e317d21b68b..c6e76ed0b65 100644 --- a/tardis/io/schemas/model.yml +++ b/tardis/io/schemas/model.yml @@ -104,6 +104,7 @@ definitions: enum: - simple_ascii - artis + - tardis_model description: file type v_inner_boundary: type: quantity diff --git a/tardis/io/setup_package.py b/tardis/io/setup_package.py index f70fa88268e..63de29adebb 100644 --- a/tardis/io/setup_package.py +++ b/tardis/io/setup_package.py @@ -1,4 +1,4 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst def get_package_data(): - return {'tardis.io.tests':['data/*.dat', 'data/*.yml']} + return {'tardis.io.tests': ['data/*.dat', 'data/*.yml', 'data/*.csv']} diff --git a/tardis/io/tests/data/tardis_configv1_isotope_uniabund.yml b/tardis/io/tests/data/tardis_configv1_isotope_uniabund.yml new file mode 100755 index 00000000000..adf02e1d765 --- /dev/null +++ b/tardis/io/tests/data/tardis_configv1_isotope_uniabund.yml @@ -0,0 +1,49 @@ +tardis_config_version: v1.0 + +supernova: + luminosity_requested: 2.8e9 solLum + time_explosion: 13 day + +atom_data: kurucz_atom_pure_simple.h5 +model: + + structure: + type: specific + + velocity: + start: 1.1e4 km/s + stop: 2.0e4 km/s + num: 20 + + density: + type: branch85_w7 + + abundances: + type: uniform + O: 0.19 + Mg: 0.03 + Si: 0.42 + S: 0.19 + Ar: 0.04 + Ca: 0.03 + Ni56: 0.05 + Ni58: 0.05 + + +plasma: + ionization: lte + excitation: lte + radiative_rates_type: dilute-blackbody + line_interaction_type: macroatom + +montecarlo: + seed: 23111963 + no_of_packets : 2.0e+5 + iterations: 5 + last_no_of_packets: 5.0e+5 + no_of_virtual_packets: 5 + +spectrum: + start: 500 angstrom + stop: 20000 angstrom + num: 10000 diff --git a/tardis/io/tests/data/tardis_configv1_tardis_model_format.yml b/tardis/io/tests/data/tardis_configv1_tardis_model_format.yml new file mode 100755 index 00000000000..b3b6a54f6c3 --- /dev/null +++ b/tardis/io/tests/data/tardis_configv1_tardis_model_format.yml @@ -0,0 +1,36 @@ +tardis_config_version: v1.0 + +supernova: + luminosity_requested: 2.8e9 solLum + time_explosion: 13 day + +atom_data: kurucz_cd23_chianti_H_He.h5 +model: + + structure: + type: file + filename: tardis_model_format.csv + filetype: tardis_model + + abundances: + type: file + filename: tardis_model_format.csv + filetype: tardis_model + +plasma: + ionization: lte + excitation: lte + radiative_rates_type: dilute-blackbody + line_interaction_type: macroatom + +montecarlo: + seed: 23111963 + no_of_packets : 2.0e+5 + iterations: 5 + last_no_of_packets: 5.0e+5 + no_of_virtual_packets: 5 + +spectrum: + start: 500 angstrom + stop: 20000 angstrom + num: 10000 diff --git a/tardis/io/tests/data/tardis_model_format.csv b/tardis/io/tests/data/tardis_model_format.csv new file mode 100644 index 00000000000..4e7dffee474 --- /dev/null +++ b/tardis/io/tests/data/tardis_model_format.csv @@ -0,0 +1,13 @@ +t0: 0.976 day +velocity temperature densities electron_densities C O Mg Si Ni56 Ni58 +km/s K g/cm^3 /cm^3 1 1 1 1 1 1 +871.66905 76395.577 4.2537191E-09 2.60E+14 0 0 0 0.6 0.4 0 +877.44269 76395.577 4.2537191E-09 2.60E+14 0 0 0 0.1 0.5 0.4 +894.99407 76395.631 4.2537191E-09 2.60E+14 0 0 0 0.3 0 0.7 +931.1571 76396.057 4.2537191E-09 2.60E+14 0 0.2 0.8 0 0 0 +990.30752 76399.042 4.2537271E-09 2.60E+14 0 0.3 0.7 0 0 0 +1050.8676 76411.983 4.2539537E-09 2.60E+14 0 0.2 0.8 0 0 0 +1115.1545 76459.592 4.2563604E-09 2.60E+14 0 0.2 0.8 0 0 0 +1183.3741 76633.367 4.2683079E-09 2.61E+14 0 0.2 0.8 0 0 0 +1255.767 77312.12 4.290997E-09 2.64E+14 0.5 0.5 0 0 0 0 +1332.5886 79602.375 4.3396835E-09 2.72E+14 0.5 0.5 0 0 0 0 diff --git a/tardis/io/tests/test_decay.py b/tardis/io/tests/test_decay.py new file mode 100644 index 00000000000..0cd1ef5bef4 --- /dev/null +++ b/tardis/io/tests/test_decay.py @@ -0,0 +1,36 @@ +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) + +@pytest.fixture +def raw_abundance_simple(): + abundances = pd.DataFrame([[0.2, 0.2], [0.1, 0.1]], index=[28, 30]) + abundances.index.rename('atomic_number', inplace=True) + return abundances + +def test_abundance_merge(simple_abundance_model, raw_abundance_simple): + decayed_df = simple_abundance_model.decay(100) + isotope_df = decayed_df.as_atoms() + combined_df = decayed_df.merge(raw_abundance_simple, normalize=False) + + assert_almost_equal(combined_df.loc[28][0], raw_abundance_simple.loc[28][0] + isotope_df.loc[28][0]) + assert_almost_equal(combined_df.loc[28][1], raw_abundance_simple.loc[28][1] + isotope_df.loc[28][1]) + assert_almost_equal(combined_df.loc[30][1], raw_abundance_simple.loc[30][1]) + assert_almost_equal(combined_df.loc[26][0], isotope_df.loc[26][0]) \ No newline at end of file diff --git a/tardis/io/tests/test_model_reader.py b/tardis/io/tests/test_model_reader.py index fd2719b0435..8bb7f08517e 100644 --- a/tardis/io/tests/test_model_reader.py +++ b/tardis/io/tests/test_model_reader.py @@ -5,8 +5,9 @@ import pytest import tardis -from tardis.io.model_reader import (read_artis_density, -read_simple_ascii_abundances) +from tardis.io.config_reader import Configuration +from tardis.io.model_reader import ( + read_artis_density, read_simple_ascii_abundances, read_simple_isotope_abundances, read_uniform_abundances, read_cmfgen_density) data_path = os.path.join(tardis.__path__[0], 'io', 'tests', 'data') @@ -18,6 +19,19 @@ def artis_density_fname(): def artis_abundances_fname(): return os.path.join(data_path, 'artis_abundances.dat') + +@pytest.fixture +def tardis_model_fname(): + return os.path.join(data_path, 'tardis_model_format.csv') + + +@pytest.fixture +def isotope_uniform_abundance(): + config_path = os.path.join( + data_path, 'tardis_configv1_isotope_uniabund.yml') + config = Configuration.from_yaml(config_path) + return config.model.abundances + def test_simple_read_artis_density(artis_density_fname): time_of_model, velocity, mean_density = read_artis_density(artis_density_fname) @@ -32,3 +46,35 @@ def test_read_simple_ascii_abundances(artis_abundances_fname): index, abundances = read_simple_ascii_abundances(artis_abundances_fname) assert len(abundances.columns) == 69 assert np.isclose(abundances[23].ix[2], 2.672351e-08 , atol=1.e-12) + + +def test_read_simple_isotope_abundances(tardis_model_fname): + index, abundances, isotope_abundance = read_simple_isotope_abundances( + tardis_model_fname) + assert np.isclose(abundances.loc[6, 8], 0.5, atol=1.e-12) + assert np.isclose(abundances.loc[12, 5], 0.8, atol=1.e-12) + assert np.isclose(abundances.loc[14, 1], 0.3, atol=1.e-12) + assert np.isclose(isotope_abundance.loc[(28, 56), 0], 0.5, atol=1.e-12) + assert np.isclose(isotope_abundance.loc[(28, 58), 1], 0.7, atol=1.e-12) + + +def test_read_uniform_abundances(isotope_uniform_abundance): + abundances, isotope_abundance = read_uniform_abundances( + isotope_uniform_abundance, 20) + assert np.isclose(abundances.loc[8, 2], 0.19, atol=1.e-12) + assert np.isclose(abundances.loc[20, 5], 0.03, atol=1.e-12) + assert np.isclose(isotope_abundance.loc[(28, 56), 15], 0.05, atol=1.e-12) + assert np.isclose(isotope_abundance.loc[(28, 58), 2], 0.05, atol=1.e-12) + + +def test_simple_read_cmfgen_density(tardis_model_fname): + time_of_model, velocity, mean_density, electron_densities, temperature = read_cmfgen_density( + tardis_model_fname) + + assert np.isclose(0.976 * u.day, time_of_model, atol=1e-7 * u.day) + assert np.isclose(mean_density[4], 4.2539537e-09 * u.g / u.cm**3, atol=1.e-6 + * u.g / u.cm**3) + assert np.isclose(electron_densities[5], 2.6e+14 * u.cm**-3, atol=1.e-6 + * u.cm**-3) + assert len(mean_density) == 9 + assert len(velocity) == len(mean_density) + 1 diff --git a/tardis/model/base.py b/tardis/model/base.py index eda7363c135..43a36058e8f 100644 --- a/tardis/model/base.py +++ b/tardis/model/base.py @@ -4,9 +4,10 @@ import pandas as pd from astropy import constants, units as u -from tardis.util import quantity_linspace, element_symbol2atomic_number -from tardis.io.model_reader import read_density_file, read_abundances_file +from tardis.util import quantity_linspace +from tardis.io.model_reader import read_density_file, read_abundances_file, read_uniform_abundances from tardis.io.util import HDFWriterMixin +from tardis.io.decay import IsotopeAbundances from density import HomologousDensity logger = logging.getLogger(__name__) @@ -59,11 +60,10 @@ class Radial1DModel(HDFWriterMixin): """ hdf_properties = ['t_inner', 'w', 't_radiative', 'v_inner', 'v_outer', 'homologous_density'] hdf_name = 'model' - - def __init__(self, velocity, homologous_density, abundance, time_explosion, - t_inner, luminosity_requested=None, t_radiative=None, - dilution_factor=None, v_boundary_inner=None, - v_boundary_outer=None): + + def __init__(self, velocity, homologous_density, abundance, isotope_abundance, + time_explosion, t_inner, luminosity_requested=None, t_radiative=None, + dilution_factor=None, v_boundary_inner=None, v_boundary_outer=None, electron_densities=None): self._v_boundary_inner = None self._v_boundary_outer = None self._velocity = None @@ -73,6 +73,11 @@ def __init__(self, velocity, homologous_density, abundance, time_explosion, self.homologous_density = homologous_density self._abundance = abundance self.time_explosion = time_explosion + self._electron_densities = electron_densities + + self.raw_abundance = self._abundance + self.raw_isotope_abundance = isotope_abundance + if t_inner is None: if luminosity_requested is not None: self.t_inner = ((luminosity_requested / @@ -185,6 +190,9 @@ def abundance(self): # abundance has one element less than velocity # ix stop index is inclusive stop -= 2 + if not self.raw_isotope_abundance.empty: + self._abundance = self.raw_isotope_abundance.decay( + self.time_explosion).merge(self.raw_abundance) abundance = self._abundance.ix[:, start:stop] abundance.columns = range(len(abundance.columns)) return abundance @@ -282,6 +290,8 @@ def from_config(cls, config): time_explosion = config.supernova.time_explosion.cgs structure = config.model.structure + electron_densities = None + temperature = None if structure.type == 'specific': velocity = quantity_linspace(structure.velocity.start, structure.velocity.stop, @@ -294,7 +304,7 @@ def from_config(cls, config): structure_fname = os.path.join(config.config_dirname, structure.filename) - time_0, velocity, density_0 = read_density_file( + time_0, velocity, density_0, electron_densities, temperature = read_density_file( structure_fname, structure.filetype) density_0 = density_0.insert(0, 0) homologous_density = HomologousDensity(density_0, time_0) @@ -304,7 +314,9 @@ def from_config(cls, config): # v boundaries. no_of_shells = len(velocity) - 1 - if config.plasma.initial_t_rad > 0 * u.K: + if temperature: + t_radiative = temperature + elif config.plasma.initial_t_rad > 0 * u.K: t_radiative = np.ones(no_of_shells) * config.plasma.initial_t_rad else: t_radiative = None @@ -317,17 +329,11 @@ def from_config(cls, config): t_inner = config.plasma.initial_t_inner abundances_section = config.model.abundances - if abundances_section.type == 'uniform': - abundance = pd.DataFrame(columns=np.arange(no_of_shells), - index=pd.Index(np.arange(1, 120), - name='atomic_number'), - dtype=np.float64) + isotope_abundance = pd.DataFrame() - for element_symbol_string in abundances_section: - if element_symbol_string == 'type': - continue - z = element_symbol2atomic_number(element_symbol_string) - abundance.ix[z] = float(abundances_section[element_symbol_string]) + if abundances_section.type == 'uniform': + abundance, isotope_abundance = read_uniform_abundances( + abundances_section, no_of_shells) elif abundances_section.type == 'file': if os.path.isabs(abundances_section.filename): @@ -336,27 +342,31 @@ def from_config(cls, config): abundances_fname = os.path.join(config.config_dirname, abundances_section.filename) - index, abundance = read_abundances_file(abundances_fname, - abundances_section.filetype) + index, abundance, isotope_abundance = read_abundances_file(abundances_fname, + abundances_section.filetype) abundance = abundance.replace(np.nan, 0.0) abundance = abundance[abundance.sum(axis=1) > 0] - norm_factor = abundance.sum(axis=0) + norm_factor = abundance.sum(axis=0) + isotope_abundance.sum(axis=0) if np.any(np.abs(norm_factor - 1) > 1e-12): logger.warning("Abundances have not been normalized to 1." " - normalizing") abundance /= norm_factor + isotope_abundance /= norm_factor + isotope_abundance = IsotopeAbundances(isotope_abundance) return cls(velocity=velocity, homologous_density=homologous_density, abundance=abundance, + isotope_abundance=isotope_abundance, time_explosion=time_explosion, t_radiative=t_radiative, t_inner=t_inner, luminosity_requested=luminosity_requested, dilution_factor=None, v_boundary_inner=structure.get('v_inner_boundary', None), - v_boundary_outer=structure.get('v_outer_boundary', None)) + v_boundary_outer=structure.get('v_outer_boundary', None), + electron_densities=electron_densities) diff --git a/tardis/model/tests/test_base.py b/tardis/model/tests/test_base.py index 9053c7423de..06122c78ff0 100644 --- a/tardis/model/tests/test_base.py +++ b/tardis/model/tests/test_base.py @@ -6,7 +6,7 @@ from tardis.io.config_reader import Configuration from tardis.model import Radial1DModel - +from tardis.io.decay import IsotopeAbundances def data_path(filename): return os.path.abspath(os.path.join('tardis/io/tests/data/', filename)) @@ -205,6 +205,35 @@ def test_ascii_reader_exponential_law(): assert_almost_equal(model.density[i].value, mdens) assert model.density[i].unit == u.Unit(expected_unit) + +@pytest.fixture +def simple_isotope_abundance(): + index = pd.MultiIndex.from_tuples([(6, 14), (12, 28)], + names=['atomic_number', 'mass_number']) + abundance = [[0.2] * 20] * 2 + return IsotopeAbundances(abundance, index=index) + + +def test_model_decay(simple_isotope_abundance): + filename = 'tardis_configv1_verysimple.yml' + config = Configuration.from_yaml(data_path(filename)) + model = Radial1DModel.from_config(config) + + model.raw_isotope_abundance = simple_isotope_abundance + decayed = simple_isotope_abundance.decay( + model.time_explosion).as_atoms() + norm_factor = 1.4 + + assert_almost_equal( + model.abundance.loc[8][0], model.raw_abundance.loc[8][0] / norm_factor, decimal=4) + assert_almost_equal(model.abundance.loc[14][0], ( + model.raw_abundance.loc[14][0] + decayed.loc[14][0]) / norm_factor, decimal=4) + assert_almost_equal(model._abundance.loc[12][5], ( + model.raw_abundance.loc[12][5] + decayed.loc[12][5]) / norm_factor, decimal=4) + assert_almost_equal( + model.abundance.loc[6][12], (decayed.loc[6][12]) / norm_factor, decimal=4) + + ### # Save and Load ### diff --git a/tardis/plasma/properties/ion_population.py b/tardis/plasma/properties/ion_population.py index 3ee8fe5913c..2a3f8d12ba1 100644 --- a/tardis/plasma/properties/ion_population.py +++ b/tardis/plasma/properties/ion_population.py @@ -184,10 +184,11 @@ class IonNumberDensity(ProcessingPlasmaProperty): outputs = ('ion_number_density', 'electron_densities') latex_name = ('N_{i,j}','n_{e}',) - def __init__(self, plasma_parent, ion_zero_threshold=1e-20): + def __init__(self, plasma_parent, ion_zero_threshold=1e-20, electron_densities=None): super(IonNumberDensity, self).__init__(plasma_parent) self.ion_zero_threshold = ion_zero_threshold self.block_ids = None + self._electron_densities = electron_densities @staticmethod def calculate_with_n_electron(phi, partition_function, @@ -223,31 +224,39 @@ def _calculate_block_ids(phi): return calculate_block_ids_from_dataframe(phi) def calculate(self, phi, partition_function, number_density): - n_e_convergence_threshold = 0.05 - n_electron = number_density.sum(axis=0) - n_electron_iterations = 0 - - while True: + if self._electron_densities is None: + n_e_convergence_threshold = 0.05 + n_electron = number_density.sum(axis=0) + n_electron_iterations = 0 + + while True: + ion_number_density, self.block_ids = \ + self.calculate_with_n_electron( + phi, partition_function, number_density, n_electron, + self.block_ids, self.ion_zero_threshold) + ion_numbers = ion_number_density.index.get_level_values(1).values + ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) + new_n_electron = (ion_number_density.values * ion_numbers).sum( + axis=0) + if np.any(np.isnan(new_n_electron)): + raise PlasmaIonizationError('n_electron just turned "nan" -' + ' aborting') + n_electron_iterations += 1 + if n_electron_iterations > 100: + logger.warn('n_electron iterations above 100 ({0}) -' + ' something is probably wrong'.format( + n_electron_iterations)) + if np.all(np.abs(new_n_electron - n_electron) + / n_electron < n_e_convergence_threshold): + break + n_electron = 0.5 * (new_n_electron + n_electron) + else: + n_electron = self._electron_densities ion_number_density, self.block_ids = \ self.calculate_with_n_electron( - phi, partition_function, number_density, n_electron, + phi, partition_function, number_density, n_electron, self.block_ids, self.ion_zero_threshold) - ion_numbers = ion_number_density.index.get_level_values(1).values - ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) - new_n_electron = (ion_number_density.values * ion_numbers).sum( - axis=0) - if np.any(np.isnan(new_n_electron)): - raise PlasmaIonizationError('n_electron just turned "nan" -' - ' aborting') - n_electron_iterations += 1 - if n_electron_iterations > 100: - logger.warn('n_electron iterations above 100 ({0}) -' - ' something is probably wrong'.format( - n_electron_iterations)) - if np.all(np.abs(new_n_electron - n_electron) - / n_electron < n_e_convergence_threshold): - break - n_electron = 0.5 * (new_n_electron + n_electron) + return ion_number_density, n_electron class IonNumberDensityHeNLTE(ProcessingPlasmaProperty): @@ -270,10 +279,11 @@ class IonNumberDensityHeNLTE(ProcessingPlasmaProperty): 'helium_population_updated') latex_name = ('N_{i,j}','n_{e}',) - def __init__(self, plasma_parent, ion_zero_threshold=1e-20): + def __init__(self, plasma_parent, ion_zero_threshold=1e-20, electron_densities=None): super(IonNumberDensityHeNLTE, self).__init__(plasma_parent) self.ion_zero_threshold = ion_zero_threshold self.block_ids = None + self._electron_densities = electron_densities def update_he_population(self, helium_population, n_electron, number_density): @@ -291,14 +301,46 @@ def update_he_population(self, helium_population, n_electron, def calculate(self, phi, partition_function, number_density, helium_population): - n_e_convergence_threshold = 0.05 - n_electron = number_density.sum(axis=0) - n_electron_iterations = 0 - while True: + if self._electron_densities is None: + n_e_convergence_threshold = 0.05 + n_electron = number_density.sum(axis=0) + n_electron_iterations = 0 + while True: + ion_number_density, self.block_ids = \ + IonNumberDensity.calculate_with_n_electron( + phi, partition_function, number_density, n_electron, + self.block_ids, self.ion_zero_threshold) + helium_population_updated = self.update_he_population( + helium_population, n_electron, number_density) + ion_number_density.ix[2].ix[0].update(helium_population_updated.ix[ + 0].sum(axis=0)) + ion_number_density.ix[2].ix[1].update(helium_population_updated.ix[ + 1].sum(axis=0)) + ion_number_density.ix[2].ix[2].update(helium_population_updated.ix[ + 2].ix[0]) + ion_numbers = ion_number_density.index.get_level_values(1).values + ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) + new_n_electron = (ion_number_density.values * ion_numbers).sum( + axis=0) + if np.any(np.isnan(new_n_electron)): + raise PlasmaIonizationError('n_electron just turned "nan" -' + ' aborting') + n_electron_iterations += 1 + if n_electron_iterations > 100: + logger.warn('n_electron iterations above 100 ({0}) -' + ' something is probably wrong'.format( + n_electron_iterations)) + if np.all(np.abs(new_n_electron - n_electron) + / n_electron < n_e_convergence_threshold): + break + n_electron = 0.5 * (new_n_electron + n_electron) + else: + n_electron = self._electron_densities ion_number_density, self.block_ids = \ IonNumberDensity.calculate_with_n_electron( - phi, partition_function, number_density, n_electron, + phi, partition_function, number_density, n_electron, self.block_ids, self.ion_zero_threshold) + helium_population_updated = self.update_he_population( helium_population, n_electron, number_density) ion_number_density.ix[2].ix[0].update(helium_population_updated.ix[ @@ -307,20 +349,4 @@ def calculate(self, phi, partition_function, number_density, 1].sum(axis=0)) ion_number_density.ix[2].ix[2].update(helium_population_updated.ix[ 2].ix[0]) - ion_numbers = ion_number_density.index.get_level_values(1).values - ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) - new_n_electron = (ion_number_density.values * ion_numbers).sum( - axis=0) - if np.any(np.isnan(new_n_electron)): - raise PlasmaIonizationError('n_electron just turned "nan" -' - ' aborting') - n_electron_iterations += 1 - if n_electron_iterations > 100: - logger.warn('n_electron iterations above 100 ({0}) -' - ' something is probably wrong'.format( - n_electron_iterations)) - if np.all(np.abs(new_n_electron - n_electron) - / n_electron < n_e_convergence_threshold): - break - n_electron = 0.5 * (new_n_electron + n_electron) return ion_number_density, n_electron, helium_population_updated diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index b9c3af21798..cc7cb7c0978 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -14,6 +14,7 @@ nlte_properties, helium_nlte_properties, helium_numerical_nlte_properties, helium_lte_properties, detailed_j_blues_properties, detailed_j_blues_inputs) from tardis.plasma.exceptions import PlasmaConfigError + from tardis.plasma.properties import ( LevelBoltzmannFactorNLTE, JBluesBlackBody, @@ -21,7 +22,8 @@ JBluesDetailed, RadiationFieldCorrection, StimulatedEmissionFactor, - HeliumNumericalNLTE) + HeliumNumericalNLTE, + IonNumberDensity) logger = logging.getLogger(__name__) @@ -148,6 +150,16 @@ def assemble_plasma(config, model, atom_data=None): heating_rate_data_file=config.plasma.heating_rate_data_file) else: plasma_modules += helium_lte_properties + + if model._electron_densities: + electron_densities = pd.Series(model._electron_densities.cgs.value) + if config.plasma.helium_treatment == 'numerical-nlte': + property_kwargs[IonNumberDensityHeNLTE] = dict( + electron_densities=electron_densities) + else: + property_kwargs[IonNumberDensity] = dict( + electron_densities=electron_densities) + kwargs['helium_treatment'] = config.plasma.helium_treatment plasma = BasePlasma(plasma_properties=plasma_modules, diff --git a/tardis/plasma/tests/test_tardis_model_density_config.py b/tardis/plasma/tests/test_tardis_model_density_config.py new file mode 100644 index 00000000000..22d0a1130b5 --- /dev/null +++ b/tardis/plasma/tests/test_tardis_model_density_config.py @@ -0,0 +1,34 @@ +import pytest +import os +from tardis.io.config_reader import Configuration +from tardis.model import Radial1DModel +from tardis.plasma.standard_plasmas import assemble_plasma +from numpy.testing import assert_almost_equal + +data_path = os.path.join('tardis', 'io', 'tests', 'data') + + +@pytest.fixture +def tardis_model_density_config(): + filename = 'tardis_configv1_tardis_model_format.yml' + return Configuration.from_yaml(os.path.join(data_path, filename)) + + +@pytest.fixture() +def raw_model(tardis_model_density_config): + return Radial1DModel.from_config(tardis_model_density_config) + + +@pytest.fixture() +def raw_plasma(tardis_model_density_config, raw_model, kurucz_atomic_data): + return assemble_plasma(tardis_model_density_config, raw_model, kurucz_atomic_data) + + +def test_electron_densities(raw_plasma): + assert_almost_equal(raw_plasma.electron_densities[8], 2.72e+14) + assert_almost_equal(raw_plasma.electron_densities[3], 2.6e+14) + + +def test_t_rad(raw_plasma): + assert_almost_equal(raw_plasma.t_rad[5], 76459.592) + assert_almost_equal(raw_plasma.t_rad[3], 76399.042) diff --git a/tardis/scripts/__init__.py b/tardis/scripts/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/scripts/cmfgen2tardis.py b/tardis/scripts/cmfgen2tardis.py new file mode 100644 index 00000000000..cef4ee18809 --- /dev/null +++ b/tardis/scripts/cmfgen2tardis.py @@ -0,0 +1,109 @@ +import os +import sys +import argparse +import numpy as np +import pandas as pd + +from tardis.atomic import AtomData + + +atomic_dataset = AtomData.from_hdf5() + + +def get_atomic_number(element): + index = -1 + for atomic_no, row in atomic_dataset.atom_data.iterrows(): + if element in row['name']: + index = atomic_no + break + return index + + +def extract_file_block(f): + qty = [] + + for line in f: + items = line.split() + if items: + qty.extend(np.array(items).astype(np.float64)) + else: + break + + qty = np.array(qty) + # Convention in CMFGEN files is different from TARDIS files. + # CMFGEN stores velocity values in decreasing order, while + # TARDIS stores it in ascending order, so alongwith + # velocity, all columns will be reversed. + return qty[::-1] + + +def convert_format(file_path): + quantities_row = [] + prop_list = ['Velocity', 'Density', 'Electron density', 'Temperature'] + with open(file_path, 'r') as f: + for line in f: + items = line.replace('(', '').replace(')', '').split() + n = len(items) + + if 'data points' in line: + abundances_df = pd.DataFrame(columns=np.arange(int(items[n - 1])), + index=pd.Index([], + name='element'), + dtype=np.float64) + if any(prop in line for prop in prop_list): + quantities_row.append(items[n - 1].replace('gm', 'g')) + if 'Time' in line: + time_of_model = float(items[n - 1]) + if 'Velocity' in line: + velocity = extract_file_block(f) + if 'Density' in line: + density = extract_file_block(f) + if 'Electron density' in line: + electron_density = extract_file_block(f) + if 'Temperature' in line: + temperature = extract_file_block(f) + + if 'mass fraction\n' in line: + element_string = items[0] + atomic_no = get_atomic_number(element_string.capitalize()) + element_symbol = atomic_dataset.atom_data.loc[atomic_no]['symbol'] + + #Its a Isotope + if n == 4: + element_symbol += items[1] + + abundances = extract_file_block(f) + abundances_df.loc[element_symbol] = abundances + + density_df = pd.DataFrame.from_records( + [velocity, temperature * 10**4, density, electron_density]).transpose() + density_df.columns = ['velocity', 'temperature', + 'densities', 'electron_densities'] + quantities_row += abundances_df.shape[0] * [1] + return abundances_df.transpose(), density_df, time_of_model, quantities_row + + +def parse_file(args): + abundances_df, density_df, time_of_model, quantities_row = convert_format( + args.input_path) + + filename = os.path.splitext(os.path.basename(args.input_path))[0] + save_fname = '.'.join((filename, 'csv')) + resultant_df = pd.concat([density_df, abundances_df], axis=1) + resultant_df.columns = pd.MultiIndex.from_tuples( + zip(resultant_df.columns, quantities_row)) + save_file_path = os.path.join(args.output_path, save_fname) + with open(save_file_path, 'w') as f: + f.write(" ".join(('t0:', str(time_of_model), 'day'))) + f.write("\n") + + resultant_df.to_csv(save_file_path, index=False, sep=' ', mode='a') + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument('input_path', help='Path to a CMFGEN file') + parser.add_argument( + 'output_path', help='Path to store converted TARDIS format files') + args = parser.parse_args() + parse_file(args) diff --git a/tardis/tests/test_util.py b/tardis/tests/test_util.py index c5d51697b6c..a027a4ffdd0 100644 --- a/tardis/tests/test_util.py +++ b/tardis/tests/test_util.py @@ -1,5 +1,6 @@ import pytest import numpy as np +import os from astropy import units as u from io import StringIO @@ -8,8 +9,14 @@ calculate_luminosity, intensity_black_body, savitzky_golay, species_tuple_to_string, species_string_to_tuple, parse_quantity, element_symbol2atomic_number, atomic_number2element_symbol, - reformat_element_symbol, quantity_linspace) + reformat_element_symbol, quantity_linspace, convert_abundances_format) +data_path = os.path.join('tardis', 'io', 'tests', 'data') + + +@pytest.fixture +def artis_abundances_fname(): + return os.path.join(data_path, 'artis_abundances.dat') def test_malformed_species_error(): malformed_species_error = MalformedSpeciesError('He') @@ -203,3 +210,11 @@ def test_quantity_linspace(start, stop, num, expected): with pytest.raises(ValueError): quantity_linspace(u.Quantity(0.5, 'eV'), '0.6 eV', 3) + + +def test_convert_abundances_format(artis_abundances_fname): + abundances = convert_abundances_format(artis_abundances_fname) + assert np.isclose(abundances.loc[3, 'O'], 1.240199e-08, atol=1.e-12) + assert np.isclose(abundances.loc[1, 'Co'], 2.306023e-05, atol=1.e-12) + assert np.isclose(abundances.loc[69, 'Ni'], 1.029928e-17, atol=1.e-12) + assert np.isclose(abundances.loc[2, 'C'], 4.425876e-09, atol=1.e-12) diff --git a/tardis/util.py b/tardis/util.py index 402a95fd4ba..6dd260bb0ef 100644 --- a/tardis/util.py +++ b/tardis/util.py @@ -1,8 +1,10 @@ # Utilities for TARDIS from astropy import units as u, constants +from pyne import nucname import numexpr as ne import numpy as np +import pandas as pd import os import yaml import re @@ -446,4 +448,14 @@ def quantity_linspace(start, stop, num, **kwargs): raise ValueError('Both start and stop need to be quantities with a ' 'unit attribute') - return np.linspace(start.value, stop.to(start.unit).value, num, **kwargs) * start.unit \ No newline at end of file + return np.linspace(start.value, stop.to(start.unit).value, num, **kwargs) * start.unit + + +def convert_abundances_format(fname, delimiter='\s+'): + df = pd.read_csv(fname, delimiter=delimiter, comment='#', header=None) + #Drop shell index column + df.drop(df.columns[0], axis=1, inplace=True) + #Assign header row + df.columns = [nucname.name(i) + for i in range(1, df.shape[1] + 1)] + return df diff --git a/tardis_env27.yml b/tardis_env27.yml index 0eeb28654a8..3f73b494b77 100644 --- a/tardis_env27.yml +++ b/tardis_env27.yml @@ -8,7 +8,7 @@ channels: dependencies: - python=2.7 - numpy=1.12 -- scipy=0.18 +- scipy=1.0 - pandas=0.20 - pytables - h5py=2.6 @@ -21,7 +21,6 @@ dependencies: - pyyaml=3.12 - jsonschema=2.5.1 - pyne=0.5.3 -- pyside=1.2 - graphviz=2.38 - pygraphviz