From 8d75f18f3fc46098174485f9d8b8dd398eb314d8 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 2 Feb 2024 23:40:13 +0000 Subject: [PATCH] Backport PR #262 on branch 0.2 (Use total recombination data if available to calculate recombination rate (#263) * Backport PR #262: Use total recombination data if available to calculate recombination rate * had to move pyarrow ignore to setup.cfg * add hashes for rate files * unpin sphinx to make rtd happy --- fiasco/conftest.py | 3 +++ fiasco/ions.py | 36 ++++++++++++++++++++++++++++++++++++ fiasco/tests/test_ion.py | 28 +++++++++++++++++++--------- setup.cfg | 3 ++- 4 files changed, 60 insertions(+), 10 deletions(-) diff --git a/fiasco/conftest.py b/fiasco/conftest.py index 365be2dc..d59a42fa 100644 --- a/fiasco/conftest.py +++ b/fiasco/conftest.py @@ -67,6 +67,9 @@ 'fe_5.elvlc': 'f7dd75c6bcea77584c5e00ddf2684d9c', 'fe_5.fblvl': 'e13a9becff489ea044fbd3cb1e02af13', 'fe_5.drparams': '0f47e2fbf4be9eb5564a8adc19531999', + 'fe_5.rrparams': '27b2f883abdc67e23eff8c544d82c7f3', + 'fe_5.diparams': 'c36dcc6d14af93f08db16428bc2957a2', + 'fe_5.trparams': '5dbfea52730e08ef9ee94cd477d83dd7', 'fe_5.easplups': 'ec5443d429722ccdd9903743ca4718ab', 'fe_5.scups': '0a863212d59c5eda5936e0664b513fb9', 'fe_5.wgfa': '6f8e4d41760d5a0540008e15aca1038c', diff --git a/fiasco/ions.py b/fiasco/ions.py index df9c9360..a6387e4f 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -1159,6 +1159,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: @@ -1172,18 +1183,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 diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 37fe4450..ce3987d6 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -287,18 +287,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): diff --git a/setup.cfg b/setup.cfg index 6c6e9e16..6b2b8cfb 100644 --- a/setup.cfg +++ b/setup.cfg @@ -47,7 +47,7 @@ test = pytest-cov matplotlib docs = - sphinx>=4.2,<5 + sphinx>=4.2 sphinx-automodapi sphinx-gallery towncrier @@ -67,6 +67,7 @@ remote_data_strict = False filterwarnings = error ignore:numpy.ndarray size changed + ignore:\nPyarrow will become a required dependency of pandas in the next major release of pandas:DeprecationWarning [coverage:run] omit =