Pull request for fixes / documentation edits related to the profiles

docs/examples/examples.rst

@@ -8,6 +8,10 @@ Here's an overview of some of the different modes of operation and models for TA
 .. toctree::
     :maxdepth: 1
 
+    profileexp
+
+    profilepl
+
     profileuniform
 
     profilemodel

docs/examples/profileexp.rst

@@ -0,0 +1,58 @@
+***********************************************
+Model with exponential density profile and uniform abundances
+***********************************************
+
+TARDIS can be used to compute a synthetic spectrum for a model with a
+user-specified density and chosen set of abundances. One simple (but
+flexible) class of model is when the density is parameterised as an
+exponential function of velocity, as described below.
+
+
+Exponential density profile
+=========================
+
+In this mode, the density profile (function of velocity and time since
+explosion) is assumed to follow a functional form:
+
+.. math::
+
+    \rho (v, t_{exp}) = \rho_0 (t_{0} / t_{exp})^{3} \exp( -v / v_0)
+
+defined by reference density, velocity and time parameters. These
+parameters are set in the input yaml file, specifically in the "structure"
+subsection of the "model" section, under the "density" heading (see
+example below).
+
+
+Uniform abundances
+============================
+
+For a model with an exponential density profile, a set of
+uniform abundances can be supplied directly in the input (yaml)
+file. Elemental 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
+
+will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
+
+.. note::
+
+    The mass fractions must sum to one. If mass fractions are supplied that do not sum to one, TARDIS will
+    renormalise all the supplied abundances and print a "WARNING" message.
+
+
+TARDIS input file
+=================
+
+Here is an example of an input file that sets up an exponential density
+profile with a uniform set of
+abundances:  
+
+.. literalinclude:: tardis_configv1_density_exponential_test.yml
+    :language: yaml
+

docs/examples/profilepl.rst

@@ -0,0 +1,58 @@
+***********************************************
+Model with power law density profile and uniform abundances
+***********************************************
+
+TARDIS can be used to compute a synthetic spectrum for a model with a
+user-specified density and chosen set of abundances. One simple (but
+flexible) class of model is when the density is parameterised as an
+power law function of velocity, as described below.
+
+
+Exponential density profile
+=========================
+
+In this mode, the density profile (function of velocity and time since
+explosion) is assumed to follow a functional form:
+
+.. math::
+
+    \rho (v, t_{exp}) = \rho_0 (t_{0} / t_{exp})^{3} ( -v / v_{0})^{\rm exponent}
+
+This form is defined by reference density, velocity and time parameters, and the
+"exponent", each of which is set in the input file ("structure"
+subsection of the "model" section, under the "density" heading; see
+example below). 
+
+
+Uniform abundances
+============================
+
+For a model with an power law density profile, a set of
+uniform abundances can be supplied directly in the input (yaml)
+file. Elemental 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
+
+will set the mass fraction of silicon (Z=14) to 0.6 and sulphur (Z=16) to 0.4.
+
+.. note::
+
+    The mass fractions must sum to one. If mass fractions are supplied that do not sum to one, TARDIS will
+    renormalise all the supplied abundances and print a "WARNING" message.
+
+
+TARDIS input file
+=================
+
+Here is an example of an input file that sets up a power law density
+profile with a uniform set of
+abundances:  
+
+.. literalinclude:: tardis_configv1_density_power_law_test.yml
+    :language: yaml
+

docs/examples/tardis_configv1_density_exponential_test.yml

@@ -0,0 +1 @@
+../../tardis/io/tests/data/tardis_configv1_density_exponential_test.yml

docs/examples/tardis_configv1_density_power_law_test.yml

@@ -0,0 +1 @@
+../../tardis/io/tests/data/tardis_configv1_density_power_law_test.yml

tardis/io/config_reader.py

@@ -47,15 +47,15 @@ def parse_quantity_linspace(quantity_linspace_dictionary, add_one=True):
      'stop': 10000 km/s,
      'num': 1000}
 
-    Parameters:
-    -----------
+    Parameters
+    ----------
 
     quantity_linspace_dictionary: ~dict
 
     add_one: boolean, default: True
 
-    Returns:
-    --------
+    Returns
+    -------
 
     ~np.array
 
@@ -85,36 +85,56 @@ def parse_spectral_bin(spectral_bin_boundary_1, spectral_bin_boundary_2):
     return spectrum_start_wavelength, spectrum_end_wavelength
 
 
-
-def calculate_exponential_densities(velocities, velocity_0, rho_0, exponent):
+def calc_exponential_density(velocities, v_0, rho0):
     """
-
-    This function computes a descret exponential density profile.
-    :math:`\\rho = \\rho_0 \\times \\left( \\frac{v_0}{v} \\right)^n`
+    This function computes the exponential density profile.
+    :math:`\\rho = \\rho_0 \\times \\exp \\left( -\\frac{v}{v_0} \\right)`
 
     Parameters
     ----------
 
-    velocities : Array like list
-                velocities in km/s
+    velocities : ~astropy.Quantity
+        Array like velocity profile
+    velocity_0 : ~astropy.Quantity
+        reference velocity
+    rho0 : ~astropy.Quantity
+        reference density
 
-    velocity_0 : ~float
-        Velocity at the inner boundary
+    Returns
+    -------
 
+    densities : ~astropy.Quantity
 
-    rho_0 : ~float
-        density at velocity_0
+    """
+    densities =  rho0 * np.exp(-(velocities / v_0))
+    return densities
 
+
+def calc_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
     -------
 
-    Array like density structure
+    densities : ~astropy.Quantity
 
     """
-    densities = rho_0 * np.power((velocity_0 / velocities), exponent)
+    densities = rho_0 * np.power((velocities / velocity_0), exponent)
     return densities
 
 
@@ -387,26 +407,70 @@ def parse_branch85(density_dict, v_inner, v_outer, time_explosion):
 
     density_parser['branch85_w7'] = parse_branch85
 
-    def parse_exponential(density_dict, v_inner, v_outer, time_explosion):
+    def parse_power_law(density_dict, v_inner, v_outer, time_explosion):
         time_0 = density_dict.pop('time_0', 19.9999584)
         if isinstance(time_0, basestring):
-            time_0 = parse_quantity(time_0).to('s').value
+            time_0 = parse_quantity(time_0).to('s')
         else:
-            logger.debug('time_0 not supplied for density branch85 - using sensible default %g', time_0)
+            logger.debug('time_0 not supplied for density powerlaw - using sensible default %g', time_0)
         try:
-            rho_0 = float(density_dict.pop('rho_0'))
+            rho_0 = density_dict.pop('rho_0')
+            if isinstance(rho_0, basestring):
+                rho_0 = parse_quantity(rho_0)
+            else:
+                raise KeyError
         except KeyError:
-            rho_0 = 1e-2
+            rho_0 = parse_quantity('1e-2 g/cm^3')
             logger.warning('rho_o was not given in the config! Using %g', rho_0)
         try:
             exponent = density_dict.pop('exponent')
         except KeyError:
             exponent = 2
             logger.warning('exponent was not given in the config file! Using %f', exponent)
+        try:
+            v_0 = density_dict.pop('v_0')
+            if isinstance(v_0, basestring):
+                v_0 = parse_quantity(v_0).to('cm/s')
+
+        except KeyError:
+            v_0 = parse_quantity('1 cm/s')
+            logger.warning('v_0 was not given in the config file! Using %f km/s', v_0)
 
+
         velocities = 0.5 * (v_inner + v_outer)
-        densities = calculate_exponential_densities(velocities, v_inner[0], rho_0, exponent)
+        densities = calc_power_law_density(velocities, v_0, rho_0, exponent)
+        densities = calculate_density_after_time(densities, time_0, time_explosion)
+        return densities
 
+    density_parser['power_law'] = parse_power_law
+
+    def parse_exponential(density_dict, v_inner, v_outer, time_explosion):
+        time_0 = density_dict.pop('time_0', 19.9999584)
+        if isinstance(time_0, basestring):
+            time_0 = parse_quantity(time_0).to('s')
+        else:
+            logger.debug('time_0 not supplied for density exponential - using sensible default %g', time_0)
+        try:
+            rho_0 = density_dict.pop('rho_0')
+            if isinstance(rho_0, basestring):
+                rho_0 = parse_quantity(rho_0)
+            else:
+                raise KeyError
+        except KeyError:
+            rho_0 = parse_quantity('1e-2 g/cm^3')
+            logger.warning('rho_o was not given in the config! Using %g', rho_0)
+        try:
+            v_0 = density_dict.pop('v_0')
+            if isinstance(v_0, basestring):
+                v_0 = parse_quantity(v_0).to('cm/s')
+
+        except KeyError:
+            v_0 = parse_quantity('1 cm/s')
+            logger.warning('v_0 was not given in the config file! Using %f km/s', v_0)
+
+        velocities = 0.5 * (v_inner + v_outer)
+        densities = calc_exponential_density(velocities, v_0, rho_0)
+        densities = calculate_density_after_time(densities, time_0, time_explosion)
         return densities
 
     density_parser['exponential'] = parse_exponential

tardis/io/tests/data/tardis_configv1_density_exponential_test.yml

@@ -0,0 +1,52 @@
+tardis_config_version: v1.0
+
+supernova:
+    luminosity_requested: 9.44 log_lsun
+    time_explosion: 13 day
+
+atom_data: kurucz_atom_pure.h5
+
+model:
+    structure:
+        type: specific
+
+
+        velocity:
+            start : 1.1e4 km/s
+            stop : 20000 km/s
+            num: 20
+
+
+        density:
+            type : exponential
+            time_0: 2. day
+            rho_0: 6.e-10 g/cm^3
+            v_0: 3000.  km/s
+
+
+    abundances:
+        type: uniform
+        O: 0.19
+        Mg: 0.03
+        Si: 0.52
+        S: 0.19
+        Ar: 0.04
+        Ca: 0.03
+
+plasma:
+    ionization: nebular
+    excitation: dilute-lte
+    radiative_rates_type: dilute-blackbody
+    line_interaction_type: scatter
+
+montecarlo:
+    seed: 23111963
+    no_of_packets : 2.0e+5
+    iterations: 30
+    last_no_of_packets: 5.0e+5
+    no_of_virtual_packets: 5
+
+spectrum:
+    start: 500 angstrom
+    stop: 20000 angstrom
+    num: 10000