Added line_ratio and example with EIS data

docs/references.bib

@@ -620,3 +620,21 @@ @article{nikolic_suppression_2018
   urldate = {2025-08-13},
   langid = {english}
 }
+
+@article{mulay_Temperature_2017,
+  title = {Temperature and Density Structure of a Recurring Active Region Jet},
+  author = {Mulay, Sargam M. and Zanna, Giulio Del and Mason, Helen},
+  year = {2017},
+  month = feb,
+  journal = {Astronomy \& Astrophysics},
+  volume = {598},
+  pages = {A11},
+  publisher = {EDP Sciences},
+  issn = {0004-6361, 1432-0746},
+  doi = {10.1051/0004-6361/201628796},
+  urldate = {2026-05-01},
+  copyright = {\copyright{} ESO, 2017},
+  langid = {english},
+  keywords = {cited},
+}
+

examples/user_guide/plot_jet_density_figure.py

@@ -0,0 +1,104 @@
+"""
+========================================
+Fe XII Density Diagnostics with EIS Data
+========================================
+
+This example shows how to compute a density diagnostic from Hinode/EIS data
+using `fiasco` by reproducing some components of Fig. 7 from :cite:t:`mulay_Temperature_2017`.
+
+The exact density values differ from the paper because this example uses the
+current CHIANTI database through `fiasco`, together with precomputed EIS
+products generated from reprocessed data.
+"""
+import astropy.units as u
+import matplotlib.pyplot as plt
+import numpy as np
+
+from astropy.io import fits
+from astropy.visualization import quantity_support
+
+import fiasco
+
+quantity_support()
+
+###############################################################################
+# First download and read in the precomputed EIS products. These two FITS files are hosted
+# in a separate data repository so that this example does not need to run the
+# EIS fitting step with `eispac <https://eispac.readthedocs.io/en/latest/>`__.
+# The intensity maps include the wavelength-integrated Fe XII 195 Å intensity and the fitted Fe XII 186 Å and 195 Å peak intensity.
+data_base_url = 'https://media.githubusercontent.com/media/sunpy/data/main/fiasco'
+with fits.open(f'{data_base_url}/jet_footpoint_fe12_observed.fits') as observed_hdul:
+    integrated_195 = observed_hdul['OBSERVED_195'].data
+with fits.open(f'{data_base_url}/jet_footpoint_fe12_intensities.fits') as intensity_hdul:
+    intensity_186 = intensity_hdul['INTENSITY_186'].data
+    intensity_195 = intensity_hdul['INTENSITY_195'].data
+    intensity_header = intensity_hdul['INTENSITY_195'].header
+
+###############################################################################
+# Compute the theoretical Fe XII density-sensitive ratio curve at the temperature
+# at which the ionization fraction of Fe XII is highest.
+density = 10**np.arange(7, 12, 0.1) * u.cm**-3
+fe12 = fiasco.Ion('Fe XII', 1.43 * u.MK)
+ratio_curve = fiasco.line_ratio(
+    fe12,
+    [186.854, 186.887] * u.angstrom,
+    [195.119, 195.179] * u.angstrom,
+    density,
+    use_two_ion_model=False,
+).squeeze()
+
+###############################################################################
+# Map the observed Fe XII intensity ratio onto the theoretical curve to derive
+# the density map.
+observed_ratio = np.divide(
+    intensity_186,
+    intensity_195,
+    out=np.full_like(intensity_186, np.nan),
+    where=intensity_195 > 0,
+)
+density_map = np.interp(observed_ratio, ratio_curve, density, left=np.nan, right=np.nan)
+
+###############################################################################
+# Plot the theoretical ratio curve, the observed Fe XII map, and the derived
+# density map for the jet footpoint.
+fig, axes = plt.subplot_mosaic(
+    [['curve', 'curve'],
+     ['observed', 'density']],
+    figsize=(10, 8),
+    layout='constrained',
+)
+
+axes['curve'].plot(density, ratio_curve, color='black', label='Theoretical')
+axes['curve'].plot(density_map.flatten(), observed_ratio.flatten(), marker='.', ls='', label='Observed')
+axes['curve'].set_title('Fe XII $(186.854 + 186.887) / (195.119 + 195.179)$')
+axes['curve'].set_xlabel(r'Electron density [$\mathrm{cm^{-3}}$]')
+axes['curve'].set_ylabel('Ratio')
+axes['curve'].set_xscale('log')
+ax_hist = axes['curve'].twinx()
+ax_hist.hist(density_map.flatten(), bins=np.logspace(7,12,100), histtype='step', log=True)
+ax_hist.set_ylabel('Number of Pixels')
+axes['curve'].legend()
+
+observed_image = axes['observed'].imshow(integrated_195, origin='lower', cmap='magma', aspect="auto")
+axes['observed'].set_title('Observed Fe XII 195')
+axes['observed'].set_xlabel('X [Pixels]')
+axes['observed'].set_ylabel('Y [Pixels]')
+
+density_image = axes['density'].imshow(np.log10(density_map.to_value('cm-3')), origin='lower', cmap='viridis', aspect="auto")
+axes['density'].set_title('Derived Electron Density')
+axes['density'].set_xlabel('X [Pixels]')
+axes['density'].set_ylabel('Y [Pixels]')
+
+fig.colorbar(
+    observed_image,
+    cax=axes['observed'].inset_axes([1.02, 0.0, 0.04, 1.0]),
+    label='Wavelength-integrated Fe XII 195',
+)
+
+fig.colorbar(
+    density_image,
+    cax=axes['density'].inset_axes([1.02, 0.0, 0.04, 1.0]),
+    label=r'$\log_{10}(n_e / \mathrm{cm^{-3}})$',
+)
+
+plt.show()

fiasco/__init__.py

@@ -5,6 +5,7 @@
 from fiasco.elements import Element
 from fiasco.fiasco import (
     get_isoelectronic_sequence,
+    line_ratio,
     list_elements,
     list_ions,
     proton_electron_ratio,
@@ -54,6 +55,7 @@ def _get_bibtex():
     "list_elements",
     "list_ions",
     "get_isoelectronic_sequence",
+    "line_ratio",
     "proton_electron_ratio",
     "Ion",
     "Levels",

fiasco/conftest.py

@@ -1,3 +1,4 @@
+import contextlib
 import numpy as np
 import pathlib
 import pytest
@@ -8,12 +9,12 @@
 from fiasco.util import check_database, read_chianti_version
 
 # Force MPL to use non-gui backends for testing.
-try:
-    import matplotlib
-except ImportError:
-    pass
-else:
-    matplotlib.use('Agg')
+with contextlib.suppress(ImportError):
+    import matplotlib as mpl
+
+    # Force MPL to use non-gui backends for testing.
+    mpl.use("Agg")
+
 
 @pytest.fixture(scope='session')
 def ascii_dbase_tree(tmpdir_factory, request):
@@ -110,6 +111,6 @@ def _evaluate_condtion(condition_string):
 
 
 def pytest_configure(config):
-  config.addinivalue_line(
+    config.addinivalue_line(
         "markers", "requires_dbase_version(dbase_version): Skip tests based on CHIANTI database version requirements.",
-  )
+    )

fiasco/fiasco.py

@@ -18,6 +18,7 @@
     'list_ions',
     'proton_electron_ratio',
     'get_isoelectronic_sequence',
+    'line_ratio',
 ]
 
 
@@ -165,3 +166,63 @@ def proton_electron_ratio(temperature: u.K, **kwargs):
                         fill_value=(ratio[0], ratio[-1]))
 
     return u.Quantity(f_interp(temperature.value))
+
+
+@u.quantity_input(density=u.cm**(-3))
+def line_ratio(ion, numerator, denominator, density: u.cm**(-3), **kwargs):
+    """
+    Theoretical line ratio for one or more bound-bound transitions as a function of temperature and density.
+
+    Parameters
+    ----------
+    ion : `~fiasco.Ion`
+        Ion used to compute the contribution function.
+    numerator, denominator : array-like of `~astropy.units.Quantity`, `str` or mix of both
+        Bound-bound transition wavelengths or transition labels for the numerator
+        and denominator. If multiple transitions are provided, their contribution
+        functions are summed before taking the ratio. Transition labels must exactly match
+        entries in ``ion.transitions.label[ion.transitions.is_bound_bound]``.
+    density : `~astropy.units.Quantity`
+        Electron number density.
+    **kwargs
+        Passed to :meth:`~fiasco.Ion.contribution_function`.
+
+    Returns
+    -------
+    `~astropy.units.Quantity`
+        The line ratio with shape ``(n_temperature, n_density)`` for independent
+        density arrays or ``(n_temperature, 1)`` for ``couple_density_to_temperature=True``.
+        Points where the denominator vanishes are returned as `~numpy.nan`.
+    """
+    numerator_idx = _transition_indices(ion, numerator)
+    denominator_idx = _transition_indices(ion, denominator)
+    contribution = ion.contribution_function(density, **kwargs)
+    numerator_term = contribution[..., numerator_idx].sum(axis=-1)
+    denominator_term = contribution[..., denominator_idx].sum(axis=-1)
+    ratio = np.full(numerator_term.shape, np.nan)
+    np.divide(
+        numerator_term.value,
+        denominator_term.value,
+        out=ratio,
+        where=denominator_term.value!=0.0,
+    )
+    return u.Quantity(ratio, numerator_term.unit / denominator_term.unit)
+
+
+def _transition_indices(ion, transitions):
+    # NOTE: This conditional avoids implicit type conversion of Quantity arrays to strings
+    if not isinstance(transitions, list):
+        transitions = np.atleast_1d(transitions)
+    wavelengths = ion.transitions.wavelength[ion.transitions.is_bound_bound]
+    labels = ion.transitions.label[ion.transitions.is_bound_bound]
+    indices = []
+    for transition in transitions:
+        if isinstance(transition, u.Quantity):
+            index = np.argmin(np.abs((wavelengths - transition).to('Angstrom')))
+        elif isinstance(transition, str):
+            index = np.flatnonzero(labels==transition)
+            if index.size == 0:
+                raise ValueError(f'No bound-bound transition found for label "{transition}".')
+            index = index[0]
+        indices.append(index)
+    return np.array(indices, dtype=int)

fiasco/tests/data/test_file_list.json

@@ -771,12 +771,19 @@
     "fe_12.easplups",
     "fe_12.elvlc",
     "fe_12.fblvl",
+    "fe_12.psplups",
     "fe_12.rrparams",
+    "fe_12.scups",
+    "fe_12.wgfa",
     "fe_13.diparams",
     "fe_13.drparams",
     "fe_13.easplups",
+    "fe_13.elvlc",
     "fe_13.fblvl",
+    "fe_13.psplups",
     "fe_13.rrparams",
+    "fe_13.scups",
+    "fe_13.wgfa",
     "fe_14.diparams",
     "fe_14.drparams",
     "fe_14.easplups",

fiasco/tests/idl/test_idl_line_ratio.py

@@ -0,0 +1,61 @@
+"""
+IDL comparison tests for line ratio calculations.
+"""
+import astropy.units as u
+import pytest
+
+import fiasco
+
+from fiasco.tests.idl.helpers import run_idl_script
+
+
+@pytest.mark.requires_dbase_version('>= 10.1')
+def test_eis_fe_xiii_density_ratio_from_idl(idl_env, dbase_version, chianti_idl_version, hdf5_dbase_root):
+    script = """
+    {% if database_version | version_check('>=', '10.1') %}
+    abundance_subdirs = ['abundance', 'archive']
+    {% else %}
+    abundance_subdirs = 'abundance'
+    {% endif %}
+    abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
+    density = findgen(21) * 0.2 + 8.0
+    ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')
+    defsysv,'!abund_file',abund_file
+    defsysv,'!ioneq_file',ioneq_file
+    temperature = ch_tmax('fe_13', /log, ioneqname=ioneq_file)
+    em = emiss_calc(26, 13, dens=density, temp=temperature, ioneq_file=ioneq_file, abund_file=abund_file, /quiet)
+
+    k = where(em.level1 EQ 1 AND em.level2 EQ 20, nk)
+    denominator_index = k[0]
+    k = where((em.level1 EQ 3 AND em.level2 EQ 24) OR (em.level1 EQ 3 AND em.level2 EQ 25), nk)
+    numerator_index = k
+
+    ratio = total(reform(em[numerator_index].em), 2) / reform(em[denominator_index].em)
+    numerator = em[numerator_index].lambda
+    denominator = em[denominator_index].lambda
+    """
+    formatters = {
+        'density': lambda x: (10.0**x) * u.cm**(-3),
+        'temperature': lambda x: (10.0**x) * u.K,
+        'numerator': lambda x: x * u.angstrom,
+        'denominator': lambda x: x * u.angstrom,
+    }
+    idl_result = run_idl_script(
+        idl_env,
+        script,
+        {},
+        ['density', 'ratio', 'numerator', 'denominator', 'temperature'],
+        'eis_fe_xiii_density_ratio',
+        dbase_version,
+        chianti_idl_version,
+        format_func=formatters,
+    )
+    ion = fiasco.Ion('Fe XIII', idl_result['temperature'], hdf5_dbase_root=hdf5_dbase_root)
+    ratio = fiasco.line_ratio(
+        ion,
+        idl_result['numerator'],
+        idl_result['denominator'],
+        idl_result['density'],
+        use_two_ion_model=False,
+    )
+    assert u.allclose(ratio.squeeze(), idl_result['ratio'], rtol=2e-3)

fiasco/tests/test_fiasco.py

@@ -1,22 +1,22 @@
 """
 Test package level functions
 """
+import pytest
+
 import astropy.units as u
 import numpy as np
 
 import fiasco
 
 
 def test_list_elements(hdf5_dbase_root):
-    # FIXME: actually test the expected elements
     elements = fiasco.list_elements(hdf5_dbase_root)
-    assert isinstance(elements, list)
+    assert elements[:5] == ['H', 'He', 'Li', 'Be', 'B']
 
 
 def test_list_ions(hdf5_dbase_root):
-    # FIXME: actually test the expected ions
     ions = fiasco.list_ions(hdf5_dbase_root)
-    assert isinstance(ions, list)
+    assert ions[:5] == ['H 1', 'H 2', 'He 1', 'He 2', 'He 3']
 
 
 def test_get_isoelectronic_sequence(hdf5_dbase_root):
@@ -36,3 +36,35 @@ def test_proton_electron_ratio(hdf5_dbase_root):
     pe_ratio = fiasco.proton_electron_ratio(t, hdf5_dbase_root=hdf5_dbase_root)
     assert type(pe_ratio) is u.Quantity
     assert pe_ratio.shape == t.shape
+    assert np.isfinite(pe_ratio).all()
+
+
+@pytest.fixture
+def fe_13(hdf5_dbase_root):
+    return fiasco.Ion('Fe XIII', 10**np.arange(5.5, 6.5, 0.05)*u.K, hdf5_dbase_root=hdf5_dbase_root)
+
+
+@pytest.mark.parametrize(('numerator', 'denominator', 'result'),[
+    (203.795*u.AA, 202.044*u.AA, [0.03940454, 0.30392223, 0.95424066]),
+    ('3s2 3p 3d 3D2 -- 3s2 3p2 3P2', '3s2 3p 3d 3P1 -- 3s2 3p2 3P0', [0.03940454, 0.30392223, 0.95424066]),
+    ([203.795, 203.826]*u.AA, 202.044*u.AA, [0.12535491, 1.05989413, 3.60720143]),
+    (['3s2 3p 3d 3D2 -- 3s2 3p2 3P2', '3s2 3p 3d 3D3 -- 3s2 3p2 3P2'], '3s2 3p 3d 3P1 -- 3s2 3p2 3P0', [0.12535491, 1.05989413, 3.60720143]),
+    ([203.795*u.AA, '3s2 3p 3d 3D3 -- 3s2 3p2 3P2'], '3s2 3p 3d 3P1 -- 3s2 3p2 3P0', [0.12535491, 1.05989413, 3.60720143]),
+])
+def test_line_ratio_wavelengths(fe_13, numerator, denominator, result):
+    density = [1e8, 1e9, 1e10] * u.cm**(-3)
+    ratio = fiasco.line_ratio(fe_13, numerator, denominator, density, use_two_ion_model=False)
+    assert ratio.shape == fe_13.temperature.shape + density.shape
+    assert u.allclose(ratio[np.argmax(fe_13.ionization_fraction), :], result)
+
+
+def test_line_ratio_temperature(fe_13):
+    density = 1e15 * u.K * u.cm**(-3) / fe_13.temperature
+    ratio = fiasco.line_ratio(fe_13,
+                              [203.795, 203.826]*u.AA,
+                              202.044*u.AA,
+                              density,
+                              couple_density_to_temperature=True,
+                              use_two_ion_model=False)
+    assert ratio.shape == fe_13.temperature.shape + (1,)
+    assert u.allclose(ratio.squeeze()[[0, 10, 19]], [3.04999452, 1.22415167, 0.3683284])