Radiative recombination (#933)

* Remove deprecated ion_cx data

* Remove unused attributes from atom_data

* Add photoionization cross sections to atom_data

* Check if continuum interaction species are in selected atoms

* Add calculation of Saha factor

* Add photoionization properties needed for recombination rates

* Add calculation of spontaneous radiative recombination rates

* Fix zero division in Saha factor calculation

* Improve index naming scheme in docstrings

* Add description of zeta_data in docstring

tardis/io/atom_data/base.py

@@ -44,24 +44,24 @@ class AtomData(object):
 
     levels: pandas.DataFrame
         A DataFrame containing the *levels data* with:
-            index: no index;
+            index: numerical index;
             columns: atomic_number, ion_number, level_number, energy[eV], g[1], metastable.
 
     lines: pandas.DataFrame
         A DataFrame containing the *lines data* with:
-            index: no index;
+            index: numerical index;
             columns: line_id, atomic_number, ion_number, level_number_lower, level_number_upper,
                 wavelength[angstrom], nu[Hz], f_lu[1], f_ul[1], B_ul[?], B_ul[?], A_ul[1/s].
 
     macro_atom_data:
         A DataFrame containing the *macro atom data* with:
-            index: no_index;
+            index: numerical index;
             columns: atomic_number, ion_number, source_level_number, destination_level_number,
                 transition_line_id, transition_type, transition_probability;
 
     macro_atom_references:
         A DataFrame containing  the *macro atom references* with:
-            index: no_index;
+            index: numerical index;
             columns: atomic_number, ion_number, source_level_number, count_down, count_up, count_total.
 
         Refer to the docs: http://tardis.readthedocs.io/en/latest/physics/plasma/macroatom.html
@@ -74,10 +74,20 @@ class AtomData(object):
     collision_data_temperatures: np.array
         An array with the collision temperatures.
 
-    zeta_data: ?
+    zeta_data:
+        A DataFrame containing the *zeta data* for the
+        nebular ionization calculation
+        (i.e., the fraction of recombinations that go directly to the
+        ground state) with:
+            index: atomic_number, ion_charge;
+            columns: temperatures[K]
+
     synpp_refs: ?
-    ion_cx_tx_data: ?
-    ion_cx_sp_data: ?
+
+    photoionization_data: pandas.DataFrame
+        A DataFrame containing the *photoionization data* with:
+            index: numerical index;
+            columns: atomic_number, ion_number, level_number, nu[Hz], x_sect[cm^2]
 
     Attributes
     -------------
@@ -91,10 +101,9 @@ class AtomData(object):
     collision_data_temperatures: numpy.array
     zeta_data: pandas.DataFrame
     synpp_refs: pandas.DataFrame
-    ion_cx_tx_data: pandas.DataFrame
-    ion_cx_sp_data: pandas.DataFrame
     symbol2atomic_number: OrderedDict
     atomic_number2symbol OrderedDict
+    photoionization_data: pandas.DataFrame
 
     Methods
     --------
@@ -118,8 +127,8 @@ class AtomData(object):
             "collision_data",
             "collision_data_temperatures",
             "synpp_refs",
-            "ion_cx_th_data",
-            "ion_cx_sp_data"]
+            "photoionization_data"
+    ]
 
     # List of tuples of the related dataframes.
     # Either all or none of the related dataframes must be given
@@ -179,12 +188,11 @@ def from_hdf(cls, fname=None):
 
         return atom_data
 
-    def __init__(
-            self, atom_data, ionization_data, levels=None, lines=None,
-            macro_atom_data=None, macro_atom_references=None,
-            zeta_data=None, collision_data=None,
-            collision_data_temperatures=None, synpp_refs=None,
-            ion_cx_th_data=None, ion_cx_sp_data=None):
+    def __init__(self, atom_data, ionization_data, levels=None, lines=None,
+                 macro_atom_data=None, macro_atom_references=None,
+                 zeta_data=None, collision_data=None,
+                 collision_data_temperatures=None, synpp_refs=None,
+                 photoionization_data=None):
 
         self.prepared = False
 
@@ -229,8 +237,8 @@ def __init__(
         self.collision_data_temperatures = collision_data_temperatures
 
         self.synpp_refs = synpp_refs
-        self.ion_cx_th_data = ion_cx_th_data
-        self.ion_cx_sp_data = ion_cx_sp_data
+
+        self.photoionization_data = photoionization_data
 
         self._check_related()
 
@@ -317,9 +325,6 @@ def prepare_atom_data(
                 self.levels_index.loc[tmp_lines_upper2level_idx].
                 astype(np.int64).values)
 
-        self.atom_ion_index = None
-        self.levels_index2atom_ion_index = None
-
         if (
                 self.macro_atom_data_all is not None and
                 not line_interaction_type == 'scatter'):

tardis/io/schemas/plasma.yml

@@ -70,6 +70,17 @@ properties:
         type: boolean
         default: false
         description: sets all beta_sobolevs to 1
+  continuum_interaction:
+    type: object
+    default: {}
+    additionalProperties: false
+    properties:
+      species:
+        type: array
+        default: []
+        description: Species that are requested to be treated with continuum
+            interactios (radiative/collisional ionization and recombination)
+            in the format ['Si II', 'Ca I', etc.]
   helium_treatment:
     type: string
     default: none

tardis/plasma/__init__.py

@@ -1,2 +1 @@
 from tardis.plasma.base import BasePlasma
-from tardis.plasma.standard_plasmas import LTEPlasma

tardis/plasma/properties/__init__.py

@@ -7,3 +7,4 @@
 from tardis.plasma.properties.radiative_properties import *
 from tardis.plasma.properties.nlte import *
 from tardis.plasma.properties.j_blues import *
+from tardis.plasma.properties.continuum_processes import *

tardis/plasma/properties/atomic.py

@@ -11,7 +11,8 @@
 logger = logging.getLogger(__name__)
 
 __all__ = ['Levels', 'Lines', 'LinesLowerLevelIndex', 'LinesUpperLevelIndex',
-           'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData']
+           'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData',
+           'PhotoIonizationData']
 
 class Levels(BaseAtomicDataProperty):
     """
@@ -65,6 +66,48 @@ def _set_index(self, lines):
         # lines.set_index('line_id', inplace=True)
         return lines, lines['nu'], lines['f_lu'], lines['wavelength_cm']
 
+
+class PhotoIonizationData(ProcessingPlasmaProperty):
+    """
+    Attributes
+    ----------
+    photo_ion_cross_sections: Pandas DataFrame (nu, x_sect,
+                                                index=['atomic_number',
+                                                       'ion_number',
+                                                       'level_number']),
+                                                dtype float)
+                              Table of photoionization cross sections as a
+                              function of frequency.
+    photo_ion_block_references: One-dimensional Numpy Array, dtype int
+                              Indices where the photoionization data for
+                              a given level starts. Needed for calculation
+                              of recombination rates.
+    photo_ion_index: Pandas MultiIndex, dtype int
+                              Atomic, ion and level numbers for which
+                              photoionization data exists.
+    """
+    outputs = ('photo_ion_cross_sections', 'photo_ion_block_references',
+               'photo_ion_index')
+    latex_name = ('\\xi_{\\textrm{i}}(\\nu)', '', '')
+
+    def calculate(self, atomic_data, continuum_interaction_species):
+        photoionization_data = atomic_data.photoionization_data.set_index(
+            ['atomic_number', 'ion_number', 'level_number']
+        )
+        selected_species_idx = pd.IndexSlice[
+            continuum_interaction_species.get_level_values('atomic_number'),
+            continuum_interaction_species.get_level_values('ion_number'),
+            slice(None)
+        ]
+        photoionization_data = photoionization_data.loc[selected_species_idx]
+        phot_nus = photoionization_data['nu']
+        block_references = np.hstack(
+            [[0], phot_nus.groupby(level=[0,1,2]).count().values.cumsum()]
+        )
+        photo_ion_index = photoionization_data.index.unique()
+        return photoionization_data, block_references, photo_ion_index
+
+
 class LinesLowerLevelIndex(HiddenPlasmaProperty):
     """
     Attributes:
@@ -92,18 +135,6 @@ def calculate(self, levels, lines):
         lines_index = lines.index.droplevel('level_number_lower')
         return np.array(levels_index.loc[lines_index])
 
-class IonCXData(BaseAtomicDataProperty):
-    outputs = ('ion_cx_data',)
-
-    def _filter_atomic_property(self, ion_cx_data, selected_atoms):
-        return ion_cx_data[ion_cx_data.atomic_number.isin([selected_atoms]
-                                                          if np.isscalar(
-            selected_atoms)
-                                                          else selected_atoms)]
-
-    def _set_index(self, ion_cx_data):
-        return ion_cx_data.set_index(['atomic_number', 'ion_number',
-                                      'level_number'])
 
 class AtomicMass(ProcessingPlasmaProperty):
     """

tardis/plasma/properties/continuum_processes.py

@@ -0,0 +1,92 @@
+import logging
+
+import numpy as np
+import pandas as pd
+
+from numba import prange, njit
+from astropy import constants as const
+
+from tardis.plasma.properties.base import ProcessingPlasmaProperty
+
+__all__ = ['SpontRecombRateCoeff']
+
+logger = logging.getLogger(__name__)
+
+njit_dict = {'fastmath': False, 'parallel': False}
+
+
+@njit(**njit_dict)
+def integrate_array_by_blocks(f, x, block_references):
+    """
+    Integrates a function f defined at locations x over blocks
+    given in block_references.
+
+    Parameters
+    ----------
+    f : Two-dimensional Numpy Array, dtype float
+    x : One-dimensional Numpy Array, dtype float
+    block_references : One-dimensional Numpy Array, dtype int
+
+    Returns
+    -------
+    integrated : Two-dimensional Numpy Array, dtype float
+
+    """
+    integrated = np.zeros((len(block_references) - 1, f.shape[1]))
+    for i in prange(f.shape[1]):  # columns
+        for j in prange(len(integrated)):  # rows
+            start = block_references[j]
+            stop = block_references[j + 1]
+            integrated[j, i] = np.trapz(f[start:stop, i], x[start:stop])
+    return integrated
+
+
+def get_ion_multi_index(multi_index_full, next_higher=True):
+    """
+    Integrates a function f defined at locations x over blocks
+    given in block_references.
+
+    Parameters
+    ----------
+    multi_index_full : Pandas MultiIndex (atomic_number, ion_number,
+                                          level_number)
+    next_higher : bool
+        If true use ion number of next higher ion, else use ion_number from
+        multi_index_full.
+
+    Returns
+    -------
+    multi_index : Pandas MultiIndex (atomic_number, ion_number)
+
+    """
+    atomic_number = multi_index_full.get_level_values(0)
+    ion_number = multi_index_full.get_level_values(1)
+    if next_higher is True:
+        ion_number += 1
+    return pd.MultiIndex.from_arrays([atomic_number, ion_number])
+
+
+class SpontRecombRateCoeff(ProcessingPlasmaProperty):
+    """
+    Attributes
+    ----------
+    alpha_sp : Pandas DataFrame, dtype float
+               The rate coefficient for spontaneous recombination.
+    """
+    outputs = ('alpha_sp',)
+    latex_name = ('\\alpha^{\\textrm{sp}}',)
+
+    def calculate(self, photo_ion_cross_sections, t_electrons,
+                  photo_ion_block_references, photo_ion_index, phi_ik):
+        x_sect = photo_ion_cross_sections['x_sect'].values
+        nu = photo_ion_cross_sections['nu'].values
+
+        alpha_sp = (8 * np.pi * x_sect * nu ** 2 / (const.c.cgs.value) ** 2)
+        alpha_sp = alpha_sp[:, np.newaxis]
+        boltzmann_factor = np.exp(-nu[np.newaxis].T / t_electrons *
+                                  (const.h.cgs.value / const.k_B.cgs.value))
+        alpha_sp = alpha_sp * boltzmann_factor
+        alpha_sp = integrate_array_by_blocks(alpha_sp, nu,
+                                             photo_ion_block_references)
+        alpha_sp = pd.DataFrame(alpha_sp, index=photo_ion_index)
+        return alpha_sp * phi_ik

tardis/plasma/properties/general.py

@@ -1,6 +1,7 @@
 import logging
 
 import numpy as np
+import pandas as pd
 from astropy import units as u
 from tardis import constants as const
 
@@ -10,7 +11,7 @@
 
 __all__ = ['BetaRadiation', 'GElectron', 'NumberDensity', 'SelectedAtoms',
            'ElectronTemperature', 'BetaElectron', 'LuminosityInner',
-           'TimeSimulation']
+           'TimeSimulation', 'ThermalGElectron']
 
 class BetaRadiation(ProcessingPlasmaProperty):
     """
@@ -44,6 +45,22 @@ def calculate(self, beta_rad):
         return ((2 * np.pi * const.m_e.cgs.value / beta_rad) /
                 (const.h.cgs.value ** 2)) ** 1.5
 
+
+class ThermalGElectron(GElectron):
+    """
+    Attributes
+    ----------
+    thermal_g_electron : Numpy Array, dtype float
+    """
+    outputs = ('thermal_g_electron',)
+    latex_name = ('g_{\\textrm{electron_thermal}}',)
+    latex_formula = ('\\Big(\\dfrac{2\\pi m_{e}/\
+                     \\beta_{\\textrm{electron}}}{h^2}\\Big)^{3/2}',)
+
+    def calculate(self, beta_electron):
+        return super(ThermalGElectron, self).calculate(beta_electron)
+
+
 class NumberDensity(ProcessingPlasmaProperty):
     """
     Attributes
@@ -63,7 +80,7 @@ class SelectedAtoms(ProcessingPlasmaProperty):
     """
     Attributes
     ----------
-    selected_atoms : Numpy Array, dtype int
+    selected_atoms : Pandas Int64Index, dtype int
                      Atomic numbers of elements required for particular simulation
     """
     outputs = ('selected_atoms',)