Use total recombination data if available to calculate recombination …

…rate (#262)

* use trparams data if available for recombination rate

* ignore pyarrow pandas warning

* fix deprectation warning

* add tests for total recombination rate priority, explicit tests for other rates

fiasco/conftest.py

@@ -72,7 +72,10 @@
     'ca_2.elvlc',
     'fe_2.elvlc',
     'fe_5.elvlc',
+    'fe_5.diparams',
     'fe_5.drparams',
+    'fe_5.rrparams',
+    'fe_5.trparams',
     'fe_5.easplups',
     'fe_5.scups',
     'fe_5.wgfa',

fiasco/ions.py

@@ -1181,6 +1181,17 @@ def dielectronic_recombination_rate(self) -> u.cm**3 / u.s:
         else:
             raise ValueError(f"Unrecognized fit type {self._drparams['fit_type']}")
 
+    @cached_property
+    @needs_dataset('trparams')
+    @u.quantity_input
+    def _total_recombination_rate(self) -> u.cm**3 / u.s:
+        temperature_data = self._trparams['temperature'].to_value('K')
+        rate_data = self._trparams['recombination_rate'].to_value('cm3 s-1')
+        f_interp = interp1d(temperature_data, rate_data, fill_value='extrapolate', kind='cubic')
+        f_interp = PchipInterpolator(np.log10(temperature_data), np.log10(rate_data), extrapolate=True)
+        rate_interp = 10**f_interp(np.log10(self.temperature.to_value('K')))
+        return u.Quantity(rate_interp, 'cm3 s-1')
+
     @cached_property
     @u.quantity_input
     def recombination_rate(self) -> u.cm**3 / u.s:
@@ -1194,18 +1205,43 @@ def recombination_rate(self) -> u.cm**3 / u.s:
 
             \alpha_{R} = \alpha_{RR} + \alpha_{DR}
 
+        .. warning::
+
+            For most ions, this total recombination rate is computed by summing the
+            outputs of the `radiative_recombination_rate` and `dielectronic_recombination_rate` methods.
+            However, for some ions, total recombination rate data is available in the
+            so-called ``.trparams`` files. For these ions, the output of this method
+            will *not* be equal to the sum of the `dielectronic_recombination_rate` and
+            `radiative_recombination_rate` method. As such, when computing the total
+            recombination rate, this method should always be used.
+
         See Also
         --------
         radiative_recombination_rate
         dielectronic_recombination_rate
         """
+        # NOTE: If the trparams data is available, then it is prioritized over the sum
+        # of the dielectronic and radiative recombination rates. This is also how the
+        # total recombination rates are computed in IDL. The reasoning here is that the
+        # total recombination rate data, if available, is more reliable than the sum of
+        # the radiative and dielectronic recombination rates. According to P. Young, there
+        # is some slight controversy over this within some communities, but CHIANTI has chosen
+        # to prioritize this data if it exists.
+        try:
+            tr_rate = self._total_recombination_rate
+        except MissingDatasetException:
+            self.log.debug(f'No total recombination data available for {self.ion_name}.')
+        else:
+            return tr_rate
         try:
             rr_rate = self.radiative_recombination_rate
         except MissingDatasetException:
+            self.log.debug(f'No radiative recombination data available for {self.ion_name}.')
             rr_rate = u.Quantity(np.zeros(self.temperature.shape), 'cm3 s-1')
         try:
             dr_rate = self.dielectronic_recombination_rate
         except MissingDatasetException:
+            self.log.debug(f'No dielectronic recombination data available for {self.ion_name}.')
             dr_rate = u.Quantity(np.zeros(self.temperature.shape), 'cm3 s-1')
         return rr_rate + dr_rate
 

fiasco/tests/test_ion.py

@@ -286,18 +286,28 @@ def test_intensity(ion, em):
     assert intens.shape == ion.temperature.shape + (1, ) + wave_shape
 
 
-def test_excitation_autoionization_rate(ion):
-    rate = ion.excitation_autoionization_rate
+@pytest.mark.parametrize(('rate_name','answer'), [
+    # NOTE: The expected values have not been tested for correctness and
+    # are only meant to indicate whether a value has changed or not.
+    ('direct_ionization_rate', 9.448935172152884e-13*u.cm**3 / u.s),
+    ('excitation_autoionization_rate', 1.14821255e-12 * u.cm**3 / u.s),
+    ('dielectronic_recombination_rate', 1.60593802e-11 * u.cm**3 / u.s),
+    ('radiative_recombination_rate', 1.6221634159408823e-12*u.cm**3 / u.s),
+])
+def test_rates(ion, rate_name, answer):
+    rate = getattr(ion, rate_name)
     assert rate.shape == ion.temperature.shape
-    # This value has not been tested for correctness
-    assert u.allclose(rate[0], 1.14821255e-12 * u.cm**3 / u.s)
+    assert u.allclose(rate[0], answer)
 
 
-def test_dielectronic_recombination_rate(ion):
-    rate = ion.dielectronic_recombination_rate
-    assert rate.shape == ion.temperature.shape
-    # This value has not been tested for correctness
-    assert u.allclose(rate[0], 1.60593802e-11 * u.cm**3 / u.s)
+def test_total_recombination_rate_priority(ion):
+    # This tests that in cases where the total recombination data in the trparams
+    # files is available that it is prioritized over the sum of the DR and RR rates.
+    # This is the case for this ion because Fe V has total recombination rate data
+    # available.
+    recomb_rate = ion.recombination_rate
+    tot_recomb_rate = ion._total_recombination_rate
+    assert np.all(recomb_rate == tot_recomb_rate)
 
 
 def test_free_free(ion):

pyproject.toml

@@ -96,6 +96,8 @@ filterwarnings =[
   "ignore:numpy.ndarray size changed",
   "ignore:The unit 'erg' has been deprecated in the VOUnit standard.*:astropy.units.core.UnitsWarning",
   "ignore:The unit 'Angstrom' has been deprecated in the VOUnit standard.*:astropy.units.core.UnitsWarning",
+  # Can be removed when pandas 3.0.0 is released
+  "ignore:\\nPyarrow will become a required dependency of pandas in the next major release of pandas:DeprecationWarning",
 ]
 
 [tool.coverage.run]