diff --git a/docs/references.bib b/docs/references.bib index 24217dd7..b87e08af 100644 --- a/docs/references.bib +++ b/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}, +} + diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py new file mode 100644 index 00000000..c4fa7b9e --- /dev/null +++ b/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 `__. +# 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() diff --git a/fiasco/__init__.py b/fiasco/__init__.py index f5b42c57..be850670 100644 --- a/fiasco/__init__.py +++ b/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", diff --git a/fiasco/conftest.py b/fiasco/conftest.py index 4600cd47..1397daab 100644 --- a/fiasco/conftest.py +++ b/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.", - ) + ) diff --git a/fiasco/fiasco.py b/fiasco/fiasco.py index e109608d..d00e2f2f 100644 --- a/fiasco/fiasco.py +++ b/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) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 9836d0b4..47fa18ab 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/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", diff --git a/fiasco/tests/idl/data/eis_fe_xiii_density_ratio_v11.0.2.asdf b/fiasco/tests/idl/data/eis_fe_xiii_density_ratio_v11.0.2.asdf new file mode 100644 index 00000000..b5d8b213 Binary files /dev/null and b/fiasco/tests/idl/data/eis_fe_xiii_density_ratio_v11.0.2.asdf differ diff --git a/fiasco/tests/idl/test_idl_line_ratio.py b/fiasco/tests/idl/test_idl_line_ratio.py new file mode 100644 index 00000000..da2ea8cf --- /dev/null +++ b/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) diff --git a/fiasco/tests/test_fiasco.py b/fiasco/tests/test_fiasco.py index d71bb2be..aae0e391 100644 --- a/fiasco/tests/test_fiasco.py +++ b/fiasco/tests/test_fiasco.py @@ -1,6 +1,8 @@ """ Test package level functions """ +import pytest + import astropy.units as u import numpy as np @@ -8,15 +10,13 @@ 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])