Added line_ratio and example with EIS data

examples/user_guide/plot_jet_density_figure.py

@@ -0,0 +1,135 @@
+"""
+========================================
+Fe XII Density Diagnostics with EIS Data
+========================================
+
+This example shows how to compute a density diagnostic from Hinode/EIS data
+using `fiasco`.
+
+The goal is to reproduce the main elements of Fig. 7 from
+`Temperature and density structure of a recurring active region jet <https://www.aanda.org/articles/aa/full_html/2017/02/aa28796-16/aa28796-16.html>`__
+by Mulay et al. (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.utils.data import download_file
+from astropy.visualization import quantity_support
+
+import fiasco
+
+quantity_support()
+
+###############################################################################
+# First download 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/>`__.
+data_base_url = 'https://media.githubusercontent.com/media/nabobalis/data/main/fiasco'
+observed_path = download_file(
+        f'{data_base_url}/jet_footpoint_fe12_observed.fits',
+        cache=True,
+        pkgname='fiasco',
+        show_progress=False,
+)
+
+intensity_path = download_file(
+        f'{data_base_url}/jet_footpoint_fe12_intensities.fits',
+        cache=True,
+        pkgname='fiasco',
+        show_progress=False,
+)
+
+
+###############################################################################
+# Load the observed Fe XII 195 map and the fitted Fe XII 186 and 195 intensity
+# maps from the downloaded FITS files.
+with fits.open(observed_path) as observed_hdul:
+    observed_195 = observed_hdul['OBSERVED_195'].data
+
+with fits.open(intensity_path) 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.
+temperature = 10**6.2 * u.K
+temperature_grid = np.array([temperature.value]) * temperature.unit
+density = np.logspace(8, 12, 80) * u.cm**-3
+
+fe12 = fiasco.Ion('Fe XII', temperature_grid)
+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,
+) * u.dimensionless_unscaled
+is_finite = np.isfinite(ratio_curve.value)
+sort_index = np.argsort(ratio_curve.value[is_finite])
+density_map = u.Quantity(
+    np.interp(
+        observed_ratio.to_value(u.dimensionless_unscaled),
+        ratio_curve.value[is_finite][sort_index],
+        density.to_value('cm-3')[is_finite][sort_index],
+        left=np.nan,
+        right=np.nan,
+    ),
+    'cm-3',
+)
+
+###############################################################################
+# 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', lw=2)
+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')
+
+observed_image = axes['observed'].imshow(observed_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,76 @@ 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 density.
+
+    Parameters
+    ----------
+    ion : `~fiasco.Ion`
+        Ion used to compute the contribution function.
+    numerator, denominator : `~astropy.units.Quantity` or array-like of `str`
+        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 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`.
+    """
+    bound_mask = ion.transitions.is_bound_bound
+    bound_wavelength = ion.transitions.wavelength[bound_mask]
+    bound_label = ion.transitions.label[bound_mask]
+
+    def transition_indices(transitions):
+        if isinstance(transitions, u.Quantity):
+            transitions = np.atleast_1d(transitions.to(bound_wavelength.unit))
+            if transitions.size == 0:
+                raise ValueError('At least one transition must be provided.')
+            return np.array(
+                [np.argmin(np.abs(bound_wavelength - transition)) for transition in transitions],
+                dtype=int,
+            )
+
+        transitions = np.atleast_1d(transitions)
+        if transitions.size == 0:
+            raise ValueError('At least one transition must be provided.')
+        if not all(isinstance(transition, str) for transition in transitions):
+            raise TypeError('Transitions must be provided as wavelengths or transition labels.')
+
+        indices = []
+        for transition in transitions:
+            index = np.flatnonzero(bound_label == transition)
+            if index.size == 0:
+                raise ValueError(f'No bound-bound transition found for label "{transition}".')
+            indices.append(index[0])
+        return np.array(indices, dtype=int)
+
+    numerator_idx = transition_indices(numerator)
+    denominator_idx = transition_indices(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)

fiasco/tests/data/test_file_list.json

@@ -775,8 +775,12 @@
     "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",
@@ -2015,4 +2019,4 @@
     "zn_30.rrparams",
     "zn_31.rrparams"
   ]
-}
+}

fiasco/tests/idl/conftest.py

@@ -5,20 +5,12 @@
 
 @pytest.fixture(scope="session")
 def idl_env(ascii_dbase_root, request):
-    # NOTE: The reason that we return None here rather than calling pytest.skip
-    # is that the IDL tests can still be run if there is a file containing the IDL
-    # results available. Thus, we have to wait until later to decide if we want to
-    # skip the test.
     idl_executable = request.config.getoption('--idl-executable', default=None)
     idl_codebase_root = request.config.getoption('--idl-codebase-root', default=None)
-    if idl_executable is None or idl_codebase_root is None:
-        return None
     return setup_idl_environment(ascii_dbase_root, idl_codebase_root, idl_executable)
 
 
 @pytest.fixture(scope="session")
 def chianti_idl_version(request):
     idl_codebase_root = request.config.getoption('--idl-codebase-root', default=None)
-    if idl_codebase_root is None:
-        return None
     return get_chianti_idl_version(idl_codebase_root)