Restructure Radial1DModel

tardis/io/config_reader.py

@@ -9,8 +9,6 @@
 import pandas as pd
 
 import tardis
-from tardis.io.model_reader import (
-    read_density_file, calculate_density_after_time, read_abundances_file)
 from tardis.io import config_validator
 from tardis.io.util import YAMLLoader, yaml_load_file
 from tardis import atomic
@@ -51,59 +49,6 @@ def parse_spectral_bin(spectral_bin_boundary_1, spectral_bin_boundary_2):
     return spectrum_start_wavelength, spectrum_end_wavelength
 
 
-def calculate_exponential_density(velocities, v_0, rho0):
-    """
-    This function computes the exponential density profile.
-    :math:`\\rho = \\rho_0 \\times \\exp \\left( -\\frac{v}{v_0} \\right)`
-
-    Parameters
-    ----------
-
-    velocities : ~astropy.Quantity
-        Array like velocity profile
-    velocity_0 : ~astropy.Quantity
-        reference velocity
-    rho0 : ~astropy.Quantity
-        reference density
-
-    Returns
-    -------
-
-    densities : ~astropy.Quantity
-
-    """
-    densities = rho0 * np.exp(-(velocities / v_0))
-    return densities
-
-
-def calculate_power_law_density(velocities, velocity_0, rho_0, exponent):
-    """
-
-    This function computes a descret exponential density profile.
-    :math:`\\rho = \\rho_0 \\times \\left( \\frac{v}{v_0} \\right)^n`
-
-    Parameters
-    ----------
-
-    velocities : ~astropy.Quantity
-        Array like velocity profile
-    velocity_0 : ~astropy.Quantity
-        reference velocity
-    rho0 : ~astropy.Quantity
-        reference density
-    exponent : ~float
-        exponent used in the powerlaw
-
-    Returns
-    -------
-
-    densities : ~astropy.Quantity
-
-    """
-    densities = rho_0 * np.power((velocities / velocity_0), exponent)
-    return densities
-
-
 def parse_model_file_section(model_setup_file_dict, time_explosion):
     def parse_artis_model_setup_files(model_file_section_dict, time_explosion):
 
@@ -784,6 +729,8 @@ def from_config_dict(cls, config_dict, atom_data=None, test_parser=False,
         else:
             validated_config_dict = config_dict
 
+        validated_config_dict.config_dirname = config_dirname
+
         #First let's see if we can find an atom_db anywhere:
         if test_parser:
             atom_data = None
@@ -817,132 +764,10 @@ def from_config_dict(cls, config_dict, atom_data=None, test_parser=False,
             validated_config_dict['supernova']['luminosity_nu_end'] = (
                 np.inf * u.Hz)
 
-        validated_config_dict['supernova']['time_explosion'] = (
-            validated_config_dict['supernova']['time_explosion'].cgs)
-
-        validated_config_dict['supernova']['luminosity_requested'] = (
-            validated_config_dict['supernova']['luminosity_requested'].cgs)
-
-        #Parsing the model section
-
-        model_section = validated_config_dict['model']
-        v_inner = None
-        v_outer = None
-        mean_densities = None
-        abundances = None
-
-
-
-        structure_section = model_section['structure']
-
-        if structure_section['type'] == 'specific':
-            velocity = model_section['structure']['velocity']
-            start, stop, num = velocity['start'], velocity['stop'], \
-                               velocity['num']
-            num += 1
-            velocities = quantity_linspace(start, stop, num)
-
-            v_inner, v_outer = velocities[:-1], velocities[1:]
-            mean_densities = parse_density_section(
-                model_section['structure']['density'], v_inner, v_outer,
-                validated_config_dict['supernova']['time_explosion']).cgs
-
-        elif structure_section['type'] == 'file':
-            if os.path.isabs(structure_section['filename']):
-                structure_fname = structure_section['filename']
-            else:
-                structure_fname = os.path.join(config_dirname,
-                                               structure_section['filename'])
-
-            v_inner, v_outer, mean_densities, inner_boundary_index, \
-            outer_boundary_index = read_density_file(
-                structure_fname, structure_section['filetype'],
-                validated_config_dict['supernova']['time_explosion'],
-                structure_section['v_inner_boundary'],
-                structure_section['v_outer_boundary'])
-
-        r_inner = validated_config_dict['supernova']['time_explosion'] * v_inner
-        r_outer = validated_config_dict['supernova']['time_explosion'] * v_outer
-        r_middle = 0.5 * (r_inner + r_outer)
-
-        structure_section['v_inner'] = v_inner.cgs
-        structure_section['v_outer'] = v_outer.cgs
-        structure_section['mean_densities'] = mean_densities.cgs
-        no_of_shells = len(v_inner)
-        structure_section['no_of_shells'] = no_of_shells
-        structure_section['r_inner'] = r_inner.cgs
-        structure_section['r_outer'] = r_outer.cgs
-        structure_section['r_middle'] = r_middle.cgs
-        structure_section['volumes'] = ((4. / 3) * np.pi * \
-                                       (r_outer ** 3 -
-                                        r_inner ** 3)).cgs
-
-
-        #### TODO the following is legacy code and should be removed
-        validated_config_dict['structure'] = \
-            validated_config_dict['model']['structure']
-        # ^^^^^^^^^^^^^^^^
-
-
-        abundances_section = model_section['abundances']
-
-        if abundances_section['type'] == 'uniform':
-            abundances = pd.DataFrame(columns=np.arange(no_of_shells),
-                                      index=pd.Index(np.arange(1, 120), name='atomic_number'), dtype=np.float64)
-
-            for element_symbol_string in abundances_section:
-                if element_symbol_string == 'type': continue
-                z = element_symbol2atomic_number(element_symbol_string)
-                abundances.ix[z] = float(abundances_section[element_symbol_string])
-
-        elif abundances_section['type'] == 'file':
-            if os.path.isabs(abundances_section['filename']):
-                abundances_fname = abundances_section['filename']
-            else:
-                abundances_fname = os.path.join(config_dirname,
-                                                abundances_section['filename'])
-
-            index, abundances = read_abundances_file(abundances_fname,
-                                                     abundances_section['filetype'],
-                                                     inner_boundary_index, outer_boundary_index)
-            if len(index) != no_of_shells:
-                raise ConfigurationError('The abundance file specified has not the same number of cells'
-                                         'as the specified density profile')
-
-        abundances = abundances.replace(np.nan, 0.0)
-
-        abundances = abundances[abundances.sum(axis=1) > 0]
-
-        norm_factor = abundances.sum(axis=0)
-
-        if np.any(np.abs(norm_factor - 1) > 1e-12):
-            logger.warning("Abundances have not been normalized to 1. - normalizing")
-            abundances /= norm_factor
-
-        validated_config_dict['abundances'] = abundances
-
-
 
         ########### DOING PLASMA SECTION ###############
         plasma_section = validated_config_dict['plasma']
 
-        if plasma_section['initial_t_inner'] < 0.0 * u.K:
-            luminosity_requested = validated_config_dict['supernova']['luminosity_requested']
-            plasma_section['t_inner'] = ((luminosity_requested /
-                                          (4 * np.pi * r_inner[0] ** 2 *
-                                           constants.sigma_sb)) ** .25).to('K')
-            logger.info('"initial_t_inner" is not specified in the plasma '
-                        'section - initializing to %s with given luminosity',
-                        plasma_section['t_inner'])
-        else:
-            plasma_section['t_inner'] = plasma_section['initial_t_inner']
-
-        if plasma_section['initial_t_rad'] > 0 * u.K:
-            plasma_section['t_rads'] = np.ones(no_of_shells) * \
-                                       plasma_section['initial_t_rad']
-        else:
-            plasma_section['t_rads'] = None
-
         if plasma_section['disable_electron_scattering'] is False:
             logger.debug("Electron scattering switched on")
             validated_config_dict['montecarlo']['sigma_thomson'] = 6.652486e-25 / (u.cm ** 2)

tardis/io/model_reader.py

@@ -11,88 +11,41 @@
 
 from tardis.util import parse_quantity
 
+
 class ConfigurationError(Exception):
     pass
 
 
-def read_density_file(density_filename, density_filetype, time_explosion,
-                      v_inner_boundary=0.0, v_outer_boundary=np.inf):
+def read_density_file(filename, filetype):
     """
     read different density file formats
 
     Parameters
     ----------
 
-    density_filename: ~str
+    filename: ~str
         filename or path of the density file
 
-    density_filetype: ~str
+    filetype: ~str
         type of the density file
 
-    time_explosion: ~astropy.units.Quantity
-        time since explosion used to scale the density
+    Returns
+    -------
+    time_of_model: ~astropy.units.Quantity
+        time at which the model is valid
+
+    velocity: ~np.ndarray
+        the array containing the velocities
+
+    unscaled_mean_densities: ~np.ndarray
+        the array containing the densities
 
     """
     file_parsers = {'artis': read_artis_density,
                     'simple_ascii': read_simple_ascii_density}
 
-    time_of_model, index, v_inner, v_outer, unscaled_mean_densities = file_parsers[density_filetype](density_filename)
-    mean_densities = calculate_density_after_time(unscaled_mean_densities, time_of_model, time_explosion)
-
-    if v_inner_boundary > v_outer_boundary:
-        raise ConfigurationError('v_inner_boundary > v_outer_boundary '
-                                 '({0:s} > {1:s}). unphysical!'.format(
-            v_inner_boundary, v_outer_boundary))
-
-    if (not np.isclose(v_inner_boundary, 0.0 * u.km / u.s,
-                       atol=1e-8 * u.km / u.s)
-        and v_inner_boundary > v_inner[0]):
-
-        if v_inner_boundary > v_outer[-1]:
-            raise ConfigurationError('Inner boundary selected outside of model')
-
-        inner_boundary_index = v_inner.searchsorted(v_inner_boundary) - 1
-        # check for zero volume of designated first cell
-        if np.isclose(v_inner_boundary, v_inner[inner_boundary_index + 1],
-                      atol=1e-8 * u.km / u.s) and (v_inner_boundary <=
-                                                   v_inner[inner_boundary_index + 1]):
-            inner_boundary_index += 1
-
-    else:
-        inner_boundary_index = None
-        v_inner_boundary = v_inner[0]
-        logger.warning("v_inner_boundary requested too small for readin file."
-                       " Boundary shifted to match file.")
-
-    if not np.isinf(v_outer_boundary) and v_outer_boundary < v_outer[-1]:
-        outer_boundary_index = v_outer.searchsorted(v_outer_boundary) + 1
-    else:
-        outer_boundary_index = None
-        v_outer_boundary = v_outer[-1]
-        logger.warning("v_outer_boundary requested too large for readin file. Boundary shifted to match file.")
-
-
-    v_inner = v_inner[inner_boundary_index:outer_boundary_index]
-    v_inner[0] = v_inner_boundary
-
-    v_outer = v_outer[inner_boundary_index:outer_boundary_index]
-    v_outer[-1] = v_outer_boundary
-
-    mean_densities = mean_densities[inner_boundary_index:outer_boundary_index]
-
-    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 "
-                             "{2:s}".format(i, v_inner_i, v_outer_i) for i,
-                             v_inner_i, v_outer_i in
-                             zip(np.arange(len(v_outer))[invalid_volume_mask],
-                                 v_inner[invalid_volume_mask],
-                                 v_outer[invalid_volume_mask])])
-        raise ConfigurationError("Invalid volume of following cell(s):\n"
-                                 "{:s}".format(message))
-
-    return (v_inner, v_outer, mean_densities,
-            inner_boundary_index, outer_boundary_index)
+    time_of_model,velocity, unscaled_mean_densities = file_parsers[filetype](filename)
+    return time_of_model, velocity, unscaled_mean_densities
 
 def read_abundances_file(abundance_filename, abundance_filetype,
                          inner_boundary_index=None, outer_boundary_index=None):
@@ -161,12 +114,12 @@ def read_simple_ascii_density(fname):
         time_of_model_string = fh.readline().strip()
         time_of_model = parse_quantity(time_of_model_string)
 
-    data = recfromtxt(fname, skip_header=1, names=('index', 'velocity', 'density'), dtype=(int, float, float))
+    data = recfromtxt(fname, skip_header=1,
+                      names=('velocity', 'density'), dtype=(float, float))
     velocity = (data['velocity'] * u.km / u.s).to('cm/s')
-    v_inner, v_outer = velocity[:-1], velocity[1:]
     mean_density = (data['density'] * u.Unit('g/cm^3'))[1:]
 
-    return time_of_model, data['index'], v_inner, v_outer, mean_density
+    return time_of_model, velocity, mean_density
 
 def read_artis_density(fname):
     """
@@ -211,9 +164,8 @@ def read_artis_density(fname):
 
     velocity = u.Quantity(artis_model['velocities'], 'km/s').to('cm/s')
     mean_density = u.Quantity(10 ** artis_model['mean_densities_0'], 'g/cm^3')[1:]
-    v_inner, v_outer = velocity[:-1], velocity[1:]
 
-    return time_of_model, artis_model['index'], v_inner, v_outer, mean_density
+    return time_of_model, velocity, mean_density
 
 
 
@@ -245,33 +197,3 @@ def read_simple_ascii_abundances(fname):
     abundances = pd.DataFrame(data[1:,1:].transpose(), index=np.arange(1, data.shape[1]))
 
     return index, abundances
-
-
-
-
-
-def calculate_density_after_time(densities, time_0, time_explosion):
-    """
-    scale the density from an initial time of the model to the time of the explosion by ^-3
-
-    Parameters:
-    -----------
-
-    densities: ~astropy.units.Quantity
-        densities
-
-    time_0: ~astropy.units.Quantity
-        time of the model
-
-    time_explosion: ~astropy.units.Quantity
-        time to be scaled to
-
-    Returns:
-    --------
-
-    scaled_density
-    """
-
-    return densities * (time_explosion / time_0) ** -3
-
-