From 66105ed1948c606ee11c674e622f148a406204ed Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Wed, 1 Apr 2026 10:50:59 -0700 Subject: [PATCH 01/12] Add ratio curve helpers and density-mapping utilities --- fiasco/__init__.py | 2 + fiasco/fiasco.py | 84 ++++++++++++++++++++++++++ fiasco/ions.py | 69 +++++++++++++++++++++ fiasco/tests/data/eis_fe_xiii.readme | 6 ++ fiasco/tests/data/eis_fe_xiii.sav | Bin 0 -> 2440 bytes fiasco/tests/data/test_file_list.json | 6 +- fiasco/tests/test_fiasco.py | 33 ++++++++++ fiasco/tests/test_ion.py | 63 +++++++++++++++++++ 8 files changed, 262 insertions(+), 1 deletion(-) create mode 100644 fiasco/tests/data/eis_fe_xiii.readme create mode 100644 fiasco/tests/data/eis_fe_xiii.sav diff --git a/fiasco/__init__.py b/fiasco/__init__.py index f5b42c57..5f2ef79e 100644 --- a/fiasco/__init__.py +++ b/fiasco/__init__.py @@ -7,6 +7,7 @@ get_isoelectronic_sequence, list_elements, list_ions, + map_ratio_to_quantity, proton_electron_ratio, ) from fiasco.gaunt import GauntFactor @@ -54,6 +55,7 @@ def _get_bibtex(): "list_elements", "list_ions", "get_isoelectronic_sequence", + "map_ratio_to_quantity", "proton_electron_ratio", "Ion", "Levels", diff --git a/fiasco/fiasco.py b/fiasco/fiasco.py index e109608d..a7008637 100644 --- a/fiasco/fiasco.py +++ b/fiasco/fiasco.py @@ -18,6 +18,7 @@ 'list_ions', 'proton_electron_ratio', 'get_isoelectronic_sequence', + 'map_ratio_to_quantity', ] @@ -165,3 +166,86 @@ def proton_electron_ratio(temperature: u.K, **kwargs): fill_value=(ratio[0], ratio[-1])) return u.Quantity(f_interp(temperature.value)) + + +def map_ratio_to_quantity(observed_ratio, + quantity, + theoretical_ratio, + *, + bounds_error=False, + fill_value=np.nan): + """ + Map an observed line ratio to the associated physical quantity using a theoretical curve. + + Parameters + ---------- + observed_ratio : array-like or `~astropy.units.Quantity` + Observed line ratio values. Must be dimensionless. + quantity : `~astropy.units.Quantity` + Quantity grid associated with ``theoretical_ratio``, e.g. an electron density grid. + theoretical_ratio : array-like or `~astropy.units.Quantity` + Theoretical line ratio sampled on ``quantity``. Must be dimensionless and monotonic. + bounds_error : `bool`, optional + If True, raise an exception when ``observed_ratio`` falls outside of the + theoretical ratio range. By default, points outside the range are set by + ``fill_value``. + fill_value : scalar, tuple, or `~astropy.units.Quantity`, optional + Value used outside of the interpolation interval. If a quantity is given, + it must be convertible to the units of ``quantity``. + + Returns + ------- + `~astropy.units.Quantity` + Interpolated quantity with the same shape as ``observed_ratio``. + """ + observed_ratio = u.Quantity(observed_ratio, u.dimensionless_unscaled) + if not isinstance(quantity, u.Quantity): + raise TypeError('quantity must be an astropy Quantity.') + quantity = np.ravel(quantity) + theoretical_ratio = np.ravel( + u.Quantity(theoretical_ratio, u.dimensionless_unscaled).value + ) + + if quantity.shape != theoretical_ratio.shape: + raise ValueError('quantity and theoretical_ratio must have the same shape.') + + is_finite = np.isfinite(quantity.value) & np.isfinite(theoretical_ratio) + quantity = quantity[is_finite] + theoretical_ratio = theoretical_ratio[is_finite] + + if quantity.size < 2: + raise ValueError('At least two finite samples are required to interpolate a ratio curve.') + + diff = np.diff(theoretical_ratio) + nonzero = diff != 0.0 + if not np.any(nonzero): + raise ValueError('theoretical_ratio must vary over the provided quantity grid.') + if np.all(diff[nonzero] < 0.0): + quantity = quantity[::-1] + theoretical_ratio = theoretical_ratio[::-1] + elif not np.all(diff[nonzero] > 0.0): + raise ValueError( + 'theoretical_ratio must be monotonic to map ratios back to a quantity. ' + 'Restrict the grid to a monotonic interval first.' + ) + + theoretical_ratio, unique_index = np.unique(theoretical_ratio, return_index=True) + quantity = quantity[unique_index] + if quantity.size < 2: + raise ValueError('Theoretical ratio must contain at least two unique samples.') + + if isinstance(fill_value, tuple): + fill_value = tuple( + value.to_value(quantity.unit) if isinstance(value, u.Quantity) else value + for value in fill_value + ) + elif isinstance(fill_value, u.Quantity): + fill_value = fill_value.to_value(quantity.unit) + + interp = interp1d( + theoretical_ratio, + quantity.value, + bounds_error=bounds_error, + fill_value=fill_value, + ) + return u.Quantity(interp(observed_ratio.to_value(u.dimensionless_unscaled)), quantity.unit) diff --git a/fiasco/ions.py b/fiasco/ions.py index 810078d7..6943a99d 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -236,6 +236,28 @@ def transitions(self): "A `~fiasco.Transitions` object holding the information about transitions for this ion." return Transitions(self.levels, self._wgfa, n_levels=self.n_levels) + def _transition_indices(self, transitions): + """ + Resolve one or more transition wavelengths to indices in the bound-bound transition array. + + Parameters + ---------- + transitions : `~astropy.units.Quantity` + One or more bound-bound transition wavelengths. + + Returns + ------- + `numpy.ndarray` + Array of integer indices into + ``self.transitions.wavelength[self.transitions.is_bound_bound]``. + """ + bound_mask = self.transitions.is_bound_bound + wavelengths = self.transitions.wavelength[bound_mask] + transitions = np.atleast_1d(transitions.to(wavelengths.unit)) + if transitions.size == 0: + raise ValueError('At least one transition must be provided.') + return np.array([np.argmin(np.abs(wavelengths - transition)) for transition in transitions], dtype=int) + @property @u.quantity_input def ionization_fraction(self) -> u.dimensionless_unscaled: @@ -1154,6 +1176,53 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er g = term * populations[:, :, i_upper] * (A * energy) return g + @u.quantity_input + @u.quantity_input + def line_ratio(self, + numerator: u.angstrom, + denominator: u.angstrom, + density: u.cm**(-3), + temperature: u.K = None, + **kwargs): + r""" + Theoretical line ratio for one or more bound-bound transitions as a function of density. + + Parameters + ---------- + numerator, denominator : `~astropy.units.Quantity` + Transition wavelength(s) for the numerator and denominator. If multiple + transitions are provided, their contribution functions are summed before + taking the ratio. + density : `~astropy.units.Quantity` + Electron number density. + temperature : `~astropy.units.Quantity` or `None`, optional + Temperature grid on which to evaluate the ratio. If `None`, the + temperature grid of the current ion instance is used. + **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`. + """ + ion = self if temperature is None else self._new_instance(temperature=temperature) + numerator_idx = ion._transition_indices(numerator) + denominator_idx = ion._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) + @u.quantity_input def emissivity(self, density: u.cm**(-3), **kwargs) -> u.erg * u.cm**(-3) / u.s: r""" diff --git a/fiasco/tests/data/eis_fe_xiii.readme b/fiasco/tests/data/eis_fe_xiii.readme new file mode 100644 index 00000000..e373dd8b --- /dev/null +++ b/fiasco/tests/data/eis_fe_xiii.readme @@ -0,0 +1,6 @@ +Generated on 04/01/2026 with SSWIDL on V11 of the database + +``` +eis_fe_xiii=eis_chianti_dens_ratio(1) +save,filename='eis_fe_xiii.sav',eis_fe_xiii +``` diff --git a/fiasco/tests/data/eis_fe_xiii.sav b/fiasco/tests/data/eis_fe_xiii.sav new file mode 100644 index 0000000000000000000000000000000000000000..56d116da5bc25625e114794cbab4d97bc278c97c GIT binary patch literal 2440 zcmeHHYe*DP6u#>^yZNA;GTESV=xxfjJDaO+?(7UX;gGw!?)s>xVN2Hsq8VA1Wkpm( zWc?_x7ZE{;K`+Q&fBX~LLrNnkv$E0<3aJe=g0ORDb|Mz~*G~_8_k8C&_nv#_aOe68 zNi-oOmXPRN@RK&#fyoX`cED6pUd9NUw=j&Av05A!o5N~lSPN^1)=;&9ja4;jj|{=` zvbrq+C723t$XyZ!#0WAm@p6xv`=e>ONFRgq!LS~lO*4G|5&@+EO65vvBXHPH2Ot_} z!W^6qkFU*t2^{Ny=>}lX1E8jE=r2A39P5DTe!vi}8HMM89w-)_BYJ$ms6xRn<*S;} zQ{Uy_T(Ld?bb zCATR1{DQIy_VpG?D_3Ji?Gt{JA21W3G=&g58_sE3WS_MdzYEa@(B=YV7#L93G;MyQ z{~zn|{8$J3i}k^312FupMeeJ`I0x4wG!3X?nA5a>d7}W61ij5B7=l3|=G18+{z;Dj zw@H`*S`uipK${I31KNDh7J-%mS{i6@y966(cF?jw%La`n#MMd5y^g1LS9h}a`m$R! zuD&Dnmzm~cu7LsR+puqz>(i#>jQcGOuCG52cmC?%>H2Z0WzX=;T|9YFVtU<= 8') +def test_line_ratio(ion): + density = [1e9, 1e10] * u.cm**(-3) + wavelengths = ion.transitions.wavelength[ion.transitions.is_bound_bound] + ratio = ion.line_ratio(wavelengths[0], wavelengths[1], density) + assert ratio.shape == ion.temperature.shape + density.shape + + grouped_ratio = ion.line_ratio(wavelengths[:2], wavelengths[2], density) + cont_func = ion.contribution_function(density) + with np.errstate(divide='ignore', invalid='ignore'): + expected_grouped_ratio = cont_func[..., [0, 1]].sum(axis=-1) / cont_func[..., [2]].sum(axis=-1) + assert np.allclose(grouped_ratio.value, expected_grouped_ratio.to_value(grouped_ratio.unit), equal_nan=True) + + density_coupled = 1e15 * u.K * u.cm**(-3) / ion.temperature + ratio_coupled = ion.line_ratio(wavelengths[0], wavelengths[1], density_coupled, couple_density_to_temperature=True) + assert ratio_coupled.shape == ion.temperature.shape + (1,) + + ratio_formation = ion.line_ratio( + wavelengths[:2], + wavelengths[2], + density, + temperature=ion.formation_temperature, + ) + ion_formation = ion._new_instance(temperature=ion.formation_temperature) + cont_func_formation = ion_formation.contribution_function(density) + with np.errstate(divide='ignore', invalid='ignore'): + expected_ratio_formation = cont_func_formation[..., [0, 1]].sum(axis=-1) / cont_func_formation[..., [2]].sum(axis=-1) + assert ratio_formation.shape == (1,) + density.shape + assert u.allclose(ratio_formation, expected_ratio_formation) + assert np.all(np.isfinite(ratio_formation)) + + +@pytest.mark.requires_dbase_version('>= 10.1') +def test_eis_fe_xiii_density_ratio_from_idl(hdf5_dbase_root): + sav_path = Path(__file__).resolve().parents[1] / 'tests/data/eis_fe_xiii.sav' + idl_ratio = readsav(sav_path)['eis_fe_xiii'] + + log_density = np.asarray(idl_ratio['DENS'][0], dtype=float) + density = (10.0**log_density) * u.cm**(-3) + expected_ratio = np.asarray(idl_ratio['RATIO'][0], dtype=float) + temperature = (10.0**float(idl_ratio['TEMP'][0])) * u.K + + numerator = np.array([ + float(wavelength) + for wavelength in idl_ratio['NUM_WVL'][0].decode().split('+') + ]) * u.angstrom + denominator = float(idl_ratio['DEN_WVL'][0].decode()) * u.angstrom + ion_name = idl_ratio['ION'][0].decode() + + ion = fiasco.Ion(ion_name, np.array([temperature.to_value('K')]) * u.K, hdf5_dbase_root=hdf5_dbase_root) + ratio = ion.line_ratio( + numerator, + denominator, + density, + temperature=temperature, + use_two_ion_model=False, + ) + + assert np.allclose(ratio.value.squeeze(), expected_ratio, rtol=2e-3) + + @pytest.mark.requires_dbase_version('>= 8') @pytest.mark.parametrize(('density', 'use_coupling'), [ ([1e9,] * u.cm**(-3), False), From 346ebb4e8521040cb574df89935e5f6580f2a177 Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Wed, 1 Apr 2026 17:40:44 -0700 Subject: [PATCH 02/12] Added an EIS example --- .../user_guide/plot_jet_density_figure.py | 131 ++++++++++++++++++ fiasco/ions.py | 1 - 2 files changed, 131 insertions(+), 1 deletion(-) create mode 100644 examples/user_guide/plot_jet_density_figure.py 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..69260059 --- /dev/null +++ b/examples/user_guide/plot_jet_density_figure.py @@ -0,0 +1,131 @@ +""" +======================================== +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 `__ +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 `__. +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 = fe12.line_ratio( + [186.854, 186.887] * u.angstrom, + [195.119, 195.179] * u.angstrom, + density, + temperature=temperature, + 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, +) +ratio_min = np.nanmin(ratio_curve.value) +ratio_max = np.nanmax(ratio_curve.value) +observed_ratio = np.where( + (observed_ratio >= ratio_min) & (observed_ratio <= ratio_max), + observed_ratio, + np.nan, +) * u.dimensionless_unscaled +density_map = fiasco.map_ratio_to_quantity(observed_ratio, density, ratio_curve) + +############################################################################### +# 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() diff --git a/fiasco/ions.py b/fiasco/ions.py index 6943a99d..8df7b077 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -1176,7 +1176,6 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er g = term * populations[:, :, i_upper] * (A * energy) return g - @u.quantity_input @u.quantity_input def line_ratio(self, numerator: u.angstrom, From 04c3d61b1091cc05b756e3096efb91ae325b7a93 Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Wed, 1 Apr 2026 18:08:18 -0700 Subject: [PATCH 03/12] revamp some of the tests --- fiasco/tests/test_fiasco.py | 14 +++---- fiasco/tests/test_ion.py | 79 +++++++++++++++++++++++++++++-------- 2 files changed, 68 insertions(+), 25 deletions(-) diff --git a/fiasco/tests/test_fiasco.py b/fiasco/tests/test_fiasco.py index b4e95e3e..3d44498e 100644 --- a/fiasco/tests/test_fiasco.py +++ b/fiasco/tests/test_fiasco.py @@ -9,15 +9,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): @@ -32,11 +30,10 @@ def test_get_isoelectronic_sequence(hdf5_dbase_root): def test_proton_electron_ratio(hdf5_dbase_root): t = np.logspace(4, 9, 100) * u.K - # NOTE: this number will not be accurate as we are using only a subset of - # the database 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.all(np.isfinite(pe_ratio)) + assert np.all(pe_ratio > 0) + assert np.ptp(pe_ratio.value) > 0 def test_map_ratio_to_quantity(): @@ -46,7 +43,6 @@ def test_map_ratio_to_quantity(): mapped_density = fiasco.map_ratio_to_quantity(observed_ratio, density, theoretical_ratio) - assert type(mapped_density) is u.Quantity assert mapped_density.shape == (3,) assert mapped_density.unit == density.unit assert u.allclose(mapped_density[:2], [1e8, 5.5e9] * u.cm**(-3)) diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index a2a4c40b..16c877da 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -280,32 +280,79 @@ def test_contribution_function(ion): @pytest.mark.requires_dbase_version('>= 8') def test_line_ratio(ion): density = [1e9, 1e10] * u.cm**(-3) - wavelengths = ion.transitions.wavelength[ion.transitions.is_bound_bound] - ratio = ion.line_ratio(wavelengths[0], wavelengths[1], density) - assert ratio.shape == ion.temperature.shape + density.shape + numerator = 703729.77 * u.angstrom + denominator = 363372.09 * u.angstrom + grouped_numerator = [703729.77, 363372.09] * u.angstrom + grouped_denominator = 151285.93 * u.angstrom - grouped_ratio = ion.line_ratio(wavelengths[:2], wavelengths[2], density) - cont_func = ion.contribution_function(density) - with np.errstate(divide='ignore', invalid='ignore'): - expected_grouped_ratio = cont_func[..., [0, 1]].sum(axis=-1) / cont_func[..., [2]].sum(axis=-1) - assert np.allclose(grouped_ratio.value, expected_grouped_ratio.to_value(grouped_ratio.unit), equal_nan=True) + ratio = ion.line_ratio(numerator, denominator, density, use_two_ion_model=False) + assert ratio.shape == ion.temperature.shape + density.shape + assert np.isclose(ratio[0, 0].value, 0.040713073461089135) + assert np.isclose(ratio[0, 1].value, 0.04085437059595315) + assert np.isclose(ratio[10, 0].value, 0.04053179876351098) + assert np.isclose(ratio[10, 1].value, 0.040757339413215445) + assert np.isclose(ratio[50, 0].value, 0.0394011731373532) + assert np.isclose(ratio[50, 1].value, 0.04054035992828684) + assert np.isnan(ratio[-1, :]).all() + + grouped_ratio = ion.line_ratio( + grouped_numerator, + grouped_denominator, + density, + use_two_ion_model=False, + ) + assert np.isclose(grouped_ratio[0, 0].value, 367497.60772844875) + assert np.isclose(grouped_ratio[0, 1].value, 367496.8140678737) + assert np.isclose(grouped_ratio[10, 1].value, 366355.21439281706) + assert np.isclose(grouped_ratio[50, 0].value, 364554.8767819584) + assert np.isnan(grouped_ratio[-1, :]).all() density_coupled = 1e15 * u.K * u.cm**(-3) / ion.temperature - ratio_coupled = ion.line_ratio(wavelengths[0], wavelengths[1], density_coupled, couple_density_to_temperature=True) + ratio_coupled = ion.line_ratio( + numerator, + denominator, + density_coupled, + couple_density_to_temperature=True, + use_two_ion_model=False, + ) assert ratio_coupled.shape == ion.temperature.shape + (1,) + assert np.allclose( + ratio_coupled[[0, 10, 50], 0].value, + [0.04085437059595315, 0.040731805023382285, 0.03821161915021005], + ) + assert np.isnan(ratio_coupled[-1, 0].value) + + +@pytest.mark.requires_dbase_version('>= 8') +def test_line_ratio_formation_temperature_values(hdf5_dbase_root): + temperature = np.logspace(5.2, 8, 100) * u.K + ion = fiasco.Ion('Fe 5', temperature, hdf5_dbase_root=hdf5_dbase_root) + density = [1e9, 1e10] * u.cm**(-3) ratio_formation = ion.line_ratio( - wavelengths[:2], - wavelengths[2], + [703729.77, 363372.09] * u.angstrom, + 151285.93 * u.angstrom, density, temperature=ion.formation_temperature, + use_two_ion_model=False, + ) + ion_formation = fiasco.Ion( + 'Fe 5', + np.array([ion.formation_temperature.to_value('K')]) * u.K, + hdf5_dbase_root=hdf5_dbase_root, + ) + direct_ratio = ion_formation.line_ratio( + [703729.77, 363372.09] * u.angstrom, + 151285.93 * u.angstrom, + density, + use_two_ion_model=False, ) - ion_formation = ion._new_instance(temperature=ion.formation_temperature) - cont_func_formation = ion_formation.contribution_function(density) - with np.errstate(divide='ignore', invalid='ignore'): - expected_ratio_formation = cont_func_formation[..., [0, 1]].sum(axis=-1) / cont_func_formation[..., [2]].sum(axis=-1) assert ratio_formation.shape == (1,) + density.shape - assert u.allclose(ratio_formation, expected_ratio_formation) + assert np.allclose( + ratio_formation[0, :].value, + [366710.87156523, 366655.15338345], + ) + assert u.allclose(ratio_formation, direct_ratio) assert np.all(np.isfinite(ratio_formation)) From 4c6fee97091e17042c0e0ec692013c1d2d4deac4 Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Mon, 6 Apr 2026 12:07:11 -0700 Subject: [PATCH 04/12] Addressed overall review comments --- .../user_guide/plot_jet_density_figure.py | 24 ++-- fiasco/__init__.py | 4 +- fiasco/conftest.py | 17 +-- fiasco/fiasco.py | 134 ++++++++---------- fiasco/ions.py | 68 --------- fiasco/tests/data/eis_fe_xiii.readme | 6 - fiasco/tests/data/eis_fe_xiii.sav | Bin 2440 -> 0 bytes fiasco/tests/idl/conftest.py | 8 -- .../eis_fe_xiii_density_ratio_v11.0.2.asdf | Bin 0 -> 3374 bytes fiasco/tests/idl/helpers.py | 61 ++++++-- fiasco/tests/idl/test_idl_line_ratio.py | 68 +++++++++ fiasco/tests/test_fiasco.py | 40 +----- fiasco/tests/test_ion.py | 94 ++++-------- 13 files changed, 238 insertions(+), 286 deletions(-) delete mode 100644 fiasco/tests/data/eis_fe_xiii.readme delete mode 100644 fiasco/tests/data/eis_fe_xiii.sav create mode 100644 fiasco/tests/idl/data/eis_fe_xiii_density_ratio_v11.0.2.asdf create mode 100644 fiasco/tests/idl/test_idl_line_ratio.py diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index 69260059..0581e8df 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -64,11 +64,11 @@ density = np.logspace(8, 12, 80) * u.cm**-3 fe12 = fiasco.Ion('Fe XII', temperature_grid) -ratio_curve = fe12.line_ratio( +ratio_curve = fiasco.line_ratio( + fe12, [186.854, 186.887] * u.angstrom, [195.119, 195.179] * u.angstrom, density, - temperature=temperature, use_two_ion_model=False, ).squeeze() @@ -80,15 +80,19 @@ intensity_195, out=np.full_like(intensity_186, np.nan), where=intensity_195 > 0, -) -ratio_min = np.nanmin(ratio_curve.value) -ratio_max = np.nanmax(ratio_curve.value) -observed_ratio = np.where( - (observed_ratio >= ratio_min) & (observed_ratio <= ratio_max), - observed_ratio, - np.nan, ) * u.dimensionless_unscaled -density_map = fiasco.map_ratio_to_quantity(observed_ratio, density, ratio_curve) +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 diff --git a/fiasco/__init__.py b/fiasco/__init__.py index 5f2ef79e..be850670 100644 --- a/fiasco/__init__.py +++ b/fiasco/__init__.py @@ -5,9 +5,9 @@ from fiasco.elements import Element from fiasco.fiasco import ( get_isoelectronic_sequence, + line_ratio, list_elements, list_ions, - map_ratio_to_quantity, proton_electron_ratio, ) from fiasco.gaunt import GauntFactor @@ -55,7 +55,7 @@ def _get_bibtex(): "list_elements", "list_ions", "get_isoelectronic_sequence", - "map_ratio_to_quantity", + "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 a7008637..e88deb06 100644 --- a/fiasco/fiasco.py +++ b/fiasco/fiasco.py @@ -18,7 +18,7 @@ 'list_ions', 'proton_electron_ratio', 'get_isoelectronic_sequence', - 'map_ratio_to_quantity', + 'line_ratio', ] @@ -168,84 +168,74 @@ def proton_electron_ratio(temperature: u.K, **kwargs): return u.Quantity(f_interp(temperature.value)) -def map_ratio_to_quantity(observed_ratio, - quantity, - theoretical_ratio, - *, - bounds_error=False, - fill_value=np.nan): +@u.quantity_input(density=u.cm**(-3)) +def line_ratio(ion, + numerator, + denominator, + density: u.cm**(-3), + **kwargs): """ - Map an observed line ratio to the associated physical quantity using a theoretical curve. + Theoretical line ratio for one or more bound-bound transitions as a function of density. Parameters ---------- - observed_ratio : array-like or `~astropy.units.Quantity` - Observed line ratio values. Must be dimensionless. - quantity : `~astropy.units.Quantity` - Quantity grid associated with ``theoretical_ratio``, e.g. an electron density grid. - theoretical_ratio : array-like or `~astropy.units.Quantity` - Theoretical line ratio sampled on ``quantity``. Must be dimensionless and monotonic. - bounds_error : `bool`, optional - If True, raise an exception when ``observed_ratio`` falls outside of the - theoretical ratio range. By default, points outside the range are set by - ``fill_value``. - fill_value : scalar, tuple, or `~astropy.units.Quantity`, optional - Value used outside of the interpolation interval. If a quantity is given, - it must be convertible to the units of ``quantity``. + 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` - Interpolated quantity with the same shape as ``observed_ratio``. + 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`. """ - observed_ratio = u.Quantity(observed_ratio, u.dimensionless_unscaled) - if not isinstance(quantity, u.Quantity): - raise TypeError('quantity must be an astropy Quantity.') - quantity = np.ravel(quantity) - theoretical_ratio = np.ravel( - u.Quantity(theoretical_ratio, u.dimensionless_unscaled).value + 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, ) - - if quantity.shape != theoretical_ratio.shape: - raise ValueError('quantity and theoretical_ratio must have the same shape.') - - is_finite = np.isfinite(quantity.value) & np.isfinite(theoretical_ratio) - quantity = quantity[is_finite] - theoretical_ratio = theoretical_ratio[is_finite] - - if quantity.size < 2: - raise ValueError('At least two finite samples are required to interpolate a ratio curve.') - - diff = np.diff(theoretical_ratio) - nonzero = diff != 0.0 - if not np.any(nonzero): - raise ValueError('theoretical_ratio must vary over the provided quantity grid.') - if np.all(diff[nonzero] < 0.0): - quantity = quantity[::-1] - theoretical_ratio = theoretical_ratio[::-1] - elif not np.all(diff[nonzero] > 0.0): - raise ValueError( - 'theoretical_ratio must be monotonic to map ratios back to a quantity. ' - 'Restrict the grid to a monotonic interval first.' - ) - - theoretical_ratio, unique_index = np.unique(theoretical_ratio, return_index=True) - quantity = quantity[unique_index] - if quantity.size < 2: - raise ValueError('Theoretical ratio must contain at least two unique samples.') - - if isinstance(fill_value, tuple): - fill_value = tuple( - value.to_value(quantity.unit) if isinstance(value, u.Quantity) else value - for value in fill_value - ) - elif isinstance(fill_value, u.Quantity): - fill_value = fill_value.to_value(quantity.unit) - - interp = interp1d( - theoretical_ratio, - quantity.value, - bounds_error=bounds_error, - fill_value=fill_value, - ) - return u.Quantity(interp(observed_ratio.to_value(u.dimensionless_unscaled)), quantity.unit) + return u.Quantity(ratio, numerator_term.unit / denominator_term.unit) diff --git a/fiasco/ions.py b/fiasco/ions.py index 8df7b077..810078d7 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -236,28 +236,6 @@ def transitions(self): "A `~fiasco.Transitions` object holding the information about transitions for this ion." return Transitions(self.levels, self._wgfa, n_levels=self.n_levels) - def _transition_indices(self, transitions): - """ - Resolve one or more transition wavelengths to indices in the bound-bound transition array. - - Parameters - ---------- - transitions : `~astropy.units.Quantity` - One or more bound-bound transition wavelengths. - - Returns - ------- - `numpy.ndarray` - Array of integer indices into - ``self.transitions.wavelength[self.transitions.is_bound_bound]``. - """ - bound_mask = self.transitions.is_bound_bound - wavelengths = self.transitions.wavelength[bound_mask] - transitions = np.atleast_1d(transitions.to(wavelengths.unit)) - if transitions.size == 0: - raise ValueError('At least one transition must be provided.') - return np.array([np.argmin(np.abs(wavelengths - transition)) for transition in transitions], dtype=int) - @property @u.quantity_input def ionization_fraction(self) -> u.dimensionless_unscaled: @@ -1176,52 +1154,6 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er g = term * populations[:, :, i_upper] * (A * energy) return g - @u.quantity_input - def line_ratio(self, - numerator: u.angstrom, - denominator: u.angstrom, - density: u.cm**(-3), - temperature: u.K = None, - **kwargs): - r""" - Theoretical line ratio for one or more bound-bound transitions as a function of density. - - Parameters - ---------- - numerator, denominator : `~astropy.units.Quantity` - Transition wavelength(s) for the numerator and denominator. If multiple - transitions are provided, their contribution functions are summed before - taking the ratio. - density : `~astropy.units.Quantity` - Electron number density. - temperature : `~astropy.units.Quantity` or `None`, optional - Temperature grid on which to evaluate the ratio. If `None`, the - temperature grid of the current ion instance is used. - **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`. - """ - ion = self if temperature is None else self._new_instance(temperature=temperature) - numerator_idx = ion._transition_indices(numerator) - denominator_idx = ion._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) - @u.quantity_input def emissivity(self, density: u.cm**(-3), **kwargs) -> u.erg * u.cm**(-3) / u.s: r""" diff --git a/fiasco/tests/data/eis_fe_xiii.readme b/fiasco/tests/data/eis_fe_xiii.readme deleted file mode 100644 index e373dd8b..00000000 --- a/fiasco/tests/data/eis_fe_xiii.readme +++ /dev/null @@ -1,6 +0,0 @@ -Generated on 04/01/2026 with SSWIDL on V11 of the database - -``` -eis_fe_xiii=eis_chianti_dens_ratio(1) -save,filename='eis_fe_xiii.sav',eis_fe_xiii -``` diff --git a/fiasco/tests/data/eis_fe_xiii.sav b/fiasco/tests/data/eis_fe_xiii.sav deleted file mode 100644 index 56d116da5bc25625e114794cbab4d97bc278c97c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2440 zcmeHHYe*DP6u#>^yZNA;GTESV=xxfjJDaO+?(7UX;gGw!?)s>xVN2Hsq8VA1Wkpm( zWc?_x7ZE{;K`+Q&fBX~LLrNnkv$E0<3aJe=g0ORDb|Mz~*G~_8_k8C&_nv#_aOe68 zNi-oOmXPRN@RK&#fyoX`cED6pUd9NUw=j&Av05A!o5N~lSPN^1)=;&9ja4;jj|{=` zvbrq+C723t$XyZ!#0WAm@p6xv`=e>ONFRgq!LS~lO*4G|5&@+EO65vvBXHPH2Ot_} z!W^6qkFU*t2^{Ny=>}lX1E8jE=r2A39P5DTe!vi}8HMM89w-)_BYJ$ms6xRn<*S;} zQ{Uy_T(Ld?bb zCATR1{DQIy_VpG?D_3Ji?Gt{JA21W3G=&g58_sE3WS_MdzYEa@(B=YV7#L93G;MyQ z{~zn|{8$J3i}k^312FupMeeJ`I0x4wG!3X?nA5a>d7}W61ij5B7=l3|=G18+{z;Dj zw@H`*S`uipK${I31KNDh7J-%mS{i6@y966(cF?jw%La`n#MMd5y^g1LS9h}a`m$R! zuD&Dnmzm~cu7LsR+puqz>(i#>jQcGOuCG52cmC?%>H2Z0WzX=;T|9YFVtU<T?L(iSQG*>+p}NomHrJ9l@7&d)M) zr$1JF&|oBpFD7by(1b)Z7~o0ePrx@KL=%yih=PDoVj?xsM8SxjJ2N{wEk)o7`_Q@f z%sJor&Uf!UcXoN}K>KPG>D+XI1N7abp zLc*X%PU47YCa@vq!@jz@I#iw0bP{IsI`D=_HWXwjt7Cm4;k-7qJRQS20T;-KCb3D2 znw~&|MS?6Q?PQe5T8ZdJ2o*JjlyHG0P@qU@DG?4ABw94HTuxK$M0rhDFtx^kkPjgh zE8xu%htMcsNSc~JagOIAn|wvdpc*g%H!@C%%3d05R;KP+1%*&7Vv22+;Nl4paxyjy zI|=88aDswWDNhW_^*Do$a*`SGrY_mhP>-dST?D7;g|H=|3u;(%Ai>swktpX7$}(O{ z)`|{-L9r@gUGxYN<>E{bhIr#Zo_I&lTdW+_wbF#88aAwP*dv^&O4I=5Vvce*UT2+y z=BYQrMPa)#!Tv{E#iFTd%7nLU*9fJOXT>tYpXyFId_sAEPJ zI9SVLXr`Wn)b{Ku_Lfg}f)Y&^2^1|$l*+`8EMivX5FZ`(F%^xRE|s9rW$RUTwC(>m zYeG&nM+C#nijr=i7Bm!at$`2(+Clmv#aiWBmZm8%7c;=K%3H3 zEDO=*W?sk>S%k(U!0F>GT)@nv(}O~LXQrik+#ID^DTD^rtzvUGahXoddKHRU4}2Y+ zfnQ!y#R5@-e6$ADq6o*Mm8gk}*z}-llZ_TJ@T|~swi9Q)w(f#!atPQykeSIOlDsit zjD`Z$7B>qy5UPx4dR6iVLY1*w1tkhxUYMFX0e-onKovY54CILrje#&>Su2E)6~}sX zO9jE9Llg)jiewl<4$HY9-@vwKG#0Y-X>mjfA(p4blgt6N=zCB@s8We7zU9F%;p9k! zLR!1R5y(4MBsvKag_GGW97P@dD2m`~00`M8o|*Va2&p5keC~-9VA;qxAh)2Ap@_4a zDzMXWW3%Aa*PwI;&Ef}`(;38TJQP$@VTxTcBZM|FwfS(Of#!6Pi?!NFq*#RCy8M_?9w&gOPo8V?UPr0u_ zJx*J@Uhm<-6;8s;&S_aay(cRbER#`gJi|Uq;F~0PZ+z?JRa)1}6T-@TGN;$K`m`}3vb(T7gIy8n(-$uGAq`eTjpUh>;hhcZhJ zFG~KH-TT6}b90lwUT)hXeL6S!=Z>alFUAU~IVXcJKl;;w)I2Ibx#sSjsk`1d(_AEN zsrhFv_kRA*=G483%sgSqcxr)oxZ~@`&ZQQd-SGSQIa^ZxmFFM+Z0l<&{}acLgg)Ak z@~g)$`&Txn{F75}M+bJL{4d^fVE54Bl>eov(}Ncmr~F&0CVv?IW!CEMbgFZo8;`Et z{5;>^*}RMY&D{4)sjchS=C-aAru9Pyt?EviVlvPQ! literal 0 HcmV?d00001 diff --git a/fiasco/tests/idl/helpers.py b/fiasco/tests/idl/helpers.py index 7a4d4773..183e0fb5 100644 --- a/fiasco/tests/idl/helpers.py +++ b/fiasco/tests/idl/helpers.py @@ -23,6 +23,35 @@ ] +def get_hissw_defaults(): + """ + Read default hissw configuration, if available. + """ + try: + import hissw + except ImportError: + return {} + return hissw.defaults + + +def find_chianti_idl_root(idl_codebase_root=None): + """ + Locate the CHIANTI IDL code root. + """ + if idl_codebase_root is not None: + idl_codebase_root = pathlib.Path(idl_codebase_root) + if idl_codebase_root.name != 'idl' and (idl_codebase_root / 'idl').is_dir(): + idl_codebase_root = idl_codebase_root / 'idl' + return idl_codebase_root + + ssw_home = get_hissw_defaults().get('ssw_home') + if ssw_home is None: + return None + + chianti_idl_root = pathlib.Path(ssw_home) / 'packages' / 'chianti' / 'idl' + return chianti_idl_root if chianti_idl_root.is_dir() else None + + def get_idl_test_output_filepath(name, version): if not isinstance(version, Version): version = Version(version) @@ -79,7 +108,12 @@ def setup_idl_environment(ascii_dbase_root, import hissw except ImportError: return None - extra_paths = [d for d, _, _ in os.walk(idl_codebase_root)] + + hissw_defaults = get_hissw_defaults() + chianti_root = find_chianti_idl_root(idl_codebase_root) + extra_paths = [] + if chianti_root is not None: + extra_paths.extend(d for d, _, _ in os.walk(chianti_root)) header = f''' defsysv,'!xuvtop','{ascii_dbase_root}' defsysv,'!abund_file','' @@ -89,16 +123,19 @@ def setup_idl_environment(ascii_dbase_root, extra_paths += [get_pkg_data_path('ssw_gen_functions', package='fiasco.tests.idl')] else: extra_paths += [d for d, _, _ in os.walk(ssw_gen_root)] + idl_only = chianti_root is not None env = hissw.Environment( - idl_only=True, - idl_home=idl_executable, + ssw_packages=[] if idl_only else ['chianti'], + ssw_home=hissw_defaults.get('ssw_home'), + idl_home=idl_executable or hissw_defaults.get('idl_home'), + idl_only=idl_only, header=header, extra_paths=extra_paths, filters={'version_check': version_check}, ) try: _ = env.run('print,!xuvtop', verbose=False) - except (hissw.util.SSWIDLError, hissw.util.IDLLicenseError): + except (ValueError, hissw.util.SSWIDLError, hissw.util.IDLLicenseError): return None else: return env @@ -151,10 +188,11 @@ def run_idl_script(idl_env, # Import here so that this can be used without a hard pytest dependency import pytest pytest.skip("""To run the IDL comparison tests, you must: - 1. Specify a path to a working IDL executable, - 2. Specify a path to the CHIANTI IDL code, - 3. Install the hissw package (pip install hissw) - Without the following, you will not be able to generate new IDL comparison test results.""") + 1. Install the hissw package (pip install hissw), + 2. Configure hissw defaults in ~/.hissw/hisswrc or pass + --idl-executable/--idl-codebase-root explicitly, + 3. Have a working IDL and CHIANTI IDL setup. + Without that, you will not be able to generate new IDL comparison test results.""") # Add versions so they are accessible from within the script if needed input_args = { **input_args, @@ -182,7 +220,7 @@ def run_idl_script(idl_env, return read_idl_test_output(file_name, dbase_version) -def get_chianti_idl_version(idl_codebase_root): +def get_chianti_idl_version(idl_codebase_root=None): """ Get the version of the CHIANTI IDL code being used. @@ -193,6 +231,11 @@ def get_chianti_idl_version(idl_codebase_root): read it from the ``VERSION`` file. If that fails, raise a warning and return None. """ + idl_codebase_root = find_chianti_idl_root(idl_codebase_root) + if idl_codebase_root is None: + fiasco.log.warning('Cannot determine CHIANTI IDL version because no CHIANTI IDL root was found.') + return None + # First try to get the version from the current git tag try: import git 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..1d227c96 --- /dev/null +++ b/fiasco/tests/idl/test_idl_line_ratio.py @@ -0,0 +1,68 @@ +""" +IDL comparison tests for line ratio calculations. +""" +import astropy.units as u +import numpy as np +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), + 'ratio': lambda x: x * u.dimensionless_unscaled, + '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', + np.atleast_1d(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 np.allclose(ratio.value.squeeze(), idl_result['ratio'].value, rtol=2e-3) diff --git a/fiasco/tests/test_fiasco.py b/fiasco/tests/test_fiasco.py index 3d44498e..4180bdad 100644 --- a/fiasco/tests/test_fiasco.py +++ b/fiasco/tests/test_fiasco.py @@ -3,7 +3,6 @@ """ import astropy.units as u import numpy as np -import pytest import fiasco @@ -30,38 +29,9 @@ def test_get_isoelectronic_sequence(hdf5_dbase_root): def test_proton_electron_ratio(hdf5_dbase_root): t = np.logspace(4, 9, 100) * u.K + # NOTE: this number will not be accurate as we are using only a subset of + # the database pe_ratio = fiasco.proton_electron_ratio(t, hdf5_dbase_root=hdf5_dbase_root) - assert np.all(np.isfinite(pe_ratio)) - assert np.all(pe_ratio > 0) - assert np.ptp(pe_ratio.value) > 0 - - -def test_map_ratio_to_quantity(): - density = [1e8, 1e9, 1e10] * u.cm**(-3) - theoretical_ratio = [0.2, 0.5, 0.8] * u.dimensionless_unscaled - observed_ratio = [0.2, 0.65, 1.0] - - mapped_density = fiasco.map_ratio_to_quantity(observed_ratio, density, theoretical_ratio) - - assert mapped_density.shape == (3,) - assert mapped_density.unit == density.unit - assert u.allclose(mapped_density[:2], [1e8, 5.5e9] * u.cm**(-3)) - assert np.isnan(mapped_density[-1].value) - - mapped_density = fiasco.map_ratio_to_quantity(0.35, density, theoretical_ratio[::-1]) - assert u.allclose(mapped_density, 5.5e9 * u.cm**(-3)) - - -def test_map_ratio_to_quantity_non_monotonic_curve(): - density = [1e8, 1e9, 1e10] * u.cm**(-3) - theoretical_ratio = [0.2, 0.8, 0.5] * u.dimensionless_unscaled - - with pytest.raises(ValueError, match='must be monotonic'): - fiasco.map_ratio_to_quantity(0.35, density, theoretical_ratio) - - -def test_map_ratio_to_quantity_requires_quantity_axis(): - theoretical_ratio = [0.2, 0.5, 0.8] * u.dimensionless_unscaled - - with pytest.raises(TypeError, match='quantity must be an astropy Quantity'): - fiasco.map_ratio_to_quantity(0.35, [1e8, 1e9, 1e10], theoretical_ratio) + assert type(pe_ratio) is u.Quantity + assert pe_ratio.shape == t.shape + assert np.isfinite(pe_ratio).all() diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 16c877da..b5f340a5 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -6,8 +6,6 @@ import pytest from fractions import Fraction -from pathlib import Path -from scipy.io import readsav import fiasco @@ -284,8 +282,11 @@ def test_line_ratio(ion): denominator = 363372.09 * u.angstrom grouped_numerator = [703729.77, 363372.09] * u.angstrom grouped_denominator = 151285.93 * u.angstrom + bound_mask = ion.transitions.is_bound_bound + bound_wavelength = ion.transitions.wavelength[bound_mask] + bound_label = ion.transitions.label[bound_mask] - ratio = ion.line_ratio(numerator, denominator, density, use_two_ion_model=False) + ratio = fiasco.line_ratio(ion, numerator, denominator, density, use_two_ion_model=False) assert ratio.shape == ion.temperature.shape + density.shape assert np.isclose(ratio[0, 0].value, 0.040713073461089135) assert np.isclose(ratio[0, 1].value, 0.04085437059595315) @@ -295,7 +296,8 @@ def test_line_ratio(ion): assert np.isclose(ratio[50, 1].value, 0.04054035992828684) assert np.isnan(ratio[-1, :]).all() - grouped_ratio = ion.line_ratio( + grouped_ratio = fiasco.line_ratio( + ion, grouped_numerator, grouped_denominator, density, @@ -307,8 +309,27 @@ def test_line_ratio(ion): assert np.isclose(grouped_ratio[50, 0].value, 364554.8767819584) assert np.isnan(grouped_ratio[-1, :]).all() + label_ratio = fiasco.line_ratio( + ion, + bound_label[np.argmin(np.abs(bound_wavelength - numerator))], + bound_label[np.argmin(np.abs(bound_wavelength - denominator))], + density, + use_two_ion_model=False, + ) + assert u.allclose(label_ratio, ratio, equal_nan=True) + + grouped_label_ratio = fiasco.line_ratio( + ion, + bound_label[[np.argmin(np.abs(bound_wavelength - w)) for w in grouped_numerator]], + bound_label[np.argmin(np.abs(bound_wavelength - grouped_denominator))], + density, + use_two_ion_model=False, + ) + assert u.allclose(grouped_label_ratio, grouped_ratio, equal_nan=True) + density_coupled = 1e15 * u.K * u.cm**(-3) / ion.temperature - ratio_coupled = ion.line_ratio( + ratio_coupled = fiasco.line_ratio( + ion, numerator, denominator, density_coupled, @@ -322,69 +343,6 @@ def test_line_ratio(ion): ) assert np.isnan(ratio_coupled[-1, 0].value) - -@pytest.mark.requires_dbase_version('>= 8') -def test_line_ratio_formation_temperature_values(hdf5_dbase_root): - temperature = np.logspace(5.2, 8, 100) * u.K - ion = fiasco.Ion('Fe 5', temperature, hdf5_dbase_root=hdf5_dbase_root) - density = [1e9, 1e10] * u.cm**(-3) - - ratio_formation = ion.line_ratio( - [703729.77, 363372.09] * u.angstrom, - 151285.93 * u.angstrom, - density, - temperature=ion.formation_temperature, - use_two_ion_model=False, - ) - ion_formation = fiasco.Ion( - 'Fe 5', - np.array([ion.formation_temperature.to_value('K')]) * u.K, - hdf5_dbase_root=hdf5_dbase_root, - ) - direct_ratio = ion_formation.line_ratio( - [703729.77, 363372.09] * u.angstrom, - 151285.93 * u.angstrom, - density, - use_two_ion_model=False, - ) - assert ratio_formation.shape == (1,) + density.shape - assert np.allclose( - ratio_formation[0, :].value, - [366710.87156523, 366655.15338345], - ) - assert u.allclose(ratio_formation, direct_ratio) - assert np.all(np.isfinite(ratio_formation)) - - -@pytest.mark.requires_dbase_version('>= 10.1') -def test_eis_fe_xiii_density_ratio_from_idl(hdf5_dbase_root): - sav_path = Path(__file__).resolve().parents[1] / 'tests/data/eis_fe_xiii.sav' - idl_ratio = readsav(sav_path)['eis_fe_xiii'] - - log_density = np.asarray(idl_ratio['DENS'][0], dtype=float) - density = (10.0**log_density) * u.cm**(-3) - expected_ratio = np.asarray(idl_ratio['RATIO'][0], dtype=float) - temperature = (10.0**float(idl_ratio['TEMP'][0])) * u.K - - numerator = np.array([ - float(wavelength) - for wavelength in idl_ratio['NUM_WVL'][0].decode().split('+') - ]) * u.angstrom - denominator = float(idl_ratio['DEN_WVL'][0].decode()) * u.angstrom - ion_name = idl_ratio['ION'][0].decode() - - ion = fiasco.Ion(ion_name, np.array([temperature.to_value('K')]) * u.K, hdf5_dbase_root=hdf5_dbase_root) - ratio = ion.line_ratio( - numerator, - denominator, - density, - temperature=temperature, - use_two_ion_model=False, - ) - - assert np.allclose(ratio.value.squeeze(), expected_ratio, rtol=2e-3) - - @pytest.mark.requires_dbase_version('>= 8') @pytest.mark.parametrize(('density', 'use_coupling'), [ ([1e9,] * u.cm**(-3), False), From 69caeed9582f16287fb362ce12215ea05fcb3b7f Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Mon, 27 Apr 2026 10:50:58 -0700 Subject: [PATCH 05/12] Rename to line_ratio_density --- examples/user_guide/plot_jet_density_figure.py | 2 +- fiasco/__init__.py | 4 ++-- fiasco/fiasco.py | 4 ++-- fiasco/tests/idl/test_idl_line_ratio.py | 2 +- fiasco/tests/test_ion.py | 12 ++++++------ 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index 0581e8df..6c1698ed 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -64,7 +64,7 @@ density = np.logspace(8, 12, 80) * u.cm**-3 fe12 = fiasco.Ion('Fe XII', temperature_grid) -ratio_curve = fiasco.line_ratio( +ratio_curve = fiasco.line_ratio_density( fe12, [186.854, 186.887] * u.angstrom, [195.119, 195.179] * u.angstrom, diff --git a/fiasco/__init__.py b/fiasco/__init__.py index be850670..baa56dc6 100644 --- a/fiasco/__init__.py +++ b/fiasco/__init__.py @@ -5,7 +5,7 @@ from fiasco.elements import Element from fiasco.fiasco import ( get_isoelectronic_sequence, - line_ratio, + line_ratio_density, list_elements, list_ions, proton_electron_ratio, @@ -55,7 +55,7 @@ def _get_bibtex(): "list_elements", "list_ions", "get_isoelectronic_sequence", - "line_ratio", + "line_ratio_density", "proton_electron_ratio", "Ion", "Levels", diff --git a/fiasco/fiasco.py b/fiasco/fiasco.py index e88deb06..0f10e37e 100644 --- a/fiasco/fiasco.py +++ b/fiasco/fiasco.py @@ -18,7 +18,7 @@ 'list_ions', 'proton_electron_ratio', 'get_isoelectronic_sequence', - 'line_ratio', + 'line_ratio_density', ] @@ -169,7 +169,7 @@ def proton_electron_ratio(temperature: u.K, **kwargs): @u.quantity_input(density=u.cm**(-3)) -def line_ratio(ion, +def line_ratio_density(ion, numerator, denominator, density: u.cm**(-3), diff --git a/fiasco/tests/idl/test_idl_line_ratio.py b/fiasco/tests/idl/test_idl_line_ratio.py index 1d227c96..01248dbc 100644 --- a/fiasco/tests/idl/test_idl_line_ratio.py +++ b/fiasco/tests/idl/test_idl_line_ratio.py @@ -57,7 +57,7 @@ def test_eis_fe_xiii_density_ratio_from_idl(idl_env, dbase_version, chianti_idl_ np.atleast_1d(idl_result['temperature']), hdf5_dbase_root=hdf5_dbase_root, ) - ratio = fiasco.line_ratio( + ratio = fiasco.line_ratio_density( ion, idl_result['numerator'], idl_result['denominator'], diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index b5f340a5..d44195fc 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -276,7 +276,7 @@ def test_contribution_function(ion): @pytest.mark.requires_dbase_version('>= 8') -def test_line_ratio(ion): +def test_line_ratio_density(ion): density = [1e9, 1e10] * u.cm**(-3) numerator = 703729.77 * u.angstrom denominator = 363372.09 * u.angstrom @@ -286,7 +286,7 @@ def test_line_ratio(ion): bound_wavelength = ion.transitions.wavelength[bound_mask] bound_label = ion.transitions.label[bound_mask] - ratio = fiasco.line_ratio(ion, numerator, denominator, density, use_two_ion_model=False) + ratio = fiasco.line_ratio_density(ion, numerator, denominator, density, use_two_ion_model=False) assert ratio.shape == ion.temperature.shape + density.shape assert np.isclose(ratio[0, 0].value, 0.040713073461089135) assert np.isclose(ratio[0, 1].value, 0.04085437059595315) @@ -296,7 +296,7 @@ def test_line_ratio(ion): assert np.isclose(ratio[50, 1].value, 0.04054035992828684) assert np.isnan(ratio[-1, :]).all() - grouped_ratio = fiasco.line_ratio( + grouped_ratio = fiasco.line_ratio_density( ion, grouped_numerator, grouped_denominator, @@ -309,7 +309,7 @@ def test_line_ratio(ion): assert np.isclose(grouped_ratio[50, 0].value, 364554.8767819584) assert np.isnan(grouped_ratio[-1, :]).all() - label_ratio = fiasco.line_ratio( + label_ratio = fiasco.line_ratio_density( ion, bound_label[np.argmin(np.abs(bound_wavelength - numerator))], bound_label[np.argmin(np.abs(bound_wavelength - denominator))], @@ -318,7 +318,7 @@ def test_line_ratio(ion): ) assert u.allclose(label_ratio, ratio, equal_nan=True) - grouped_label_ratio = fiasco.line_ratio( + grouped_label_ratio = fiasco.line_ratio_density( ion, bound_label[[np.argmin(np.abs(bound_wavelength - w)) for w in grouped_numerator]], bound_label[np.argmin(np.abs(bound_wavelength - grouped_denominator))], @@ -328,7 +328,7 @@ def test_line_ratio(ion): assert u.allclose(grouped_label_ratio, grouped_ratio, equal_nan=True) density_coupled = 1e15 * u.K * u.cm**(-3) / ion.temperature - ratio_coupled = fiasco.line_ratio( + ratio_coupled = fiasco.line_ratio_density( ion, numerator, denominator, From 6b40d8c82211c522f6d96a9c1da288a184bfcabc Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 1 May 2026 17:58:46 -0400 Subject: [PATCH 06/12] simplified and renamed ratio function --- fiasco/__init__.py | 2 +- fiasco/fiasco.py | 67 +++++++++++++++--------------------- fiasco/tests/idl/conftest.py | 8 +++++ fiasco/tests/idl/helpers.py | 61 +++++--------------------------- 4 files changed, 45 insertions(+), 93 deletions(-) diff --git a/fiasco/__init__.py b/fiasco/__init__.py index baa56dc6..2bc9068d 100644 --- a/fiasco/__init__.py +++ b/fiasco/__init__.py @@ -5,7 +5,7 @@ from fiasco.elements import Element from fiasco.fiasco import ( get_isoelectronic_sequence, - line_ratio_density, + line_ratio, list_elements, list_ions, proton_electron_ratio, diff --git a/fiasco/fiasco.py b/fiasco/fiasco.py index 0f10e37e..d00e2f2f 100644 --- a/fiasco/fiasco.py +++ b/fiasco/fiasco.py @@ -18,7 +18,7 @@ 'list_ions', 'proton_electron_ratio', 'get_isoelectronic_sequence', - 'line_ratio_density', + 'line_ratio', ] @@ -169,22 +169,18 @@ def proton_electron_ratio(temperature: u.K, **kwargs): @u.quantity_input(density=u.cm**(-3)) -def line_ratio_density(ion, - numerator, - denominator, - density: u.cm**(-3), - **kwargs): +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. + 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 : `~astropy.units.Quantity` or array-like of `str` + 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 match + 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. @@ -198,36 +194,8 @@ def line_ratio_density(ion, 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) + 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) @@ -236,6 +204,25 @@ def transition_indices(transitions): numerator_term.value, denominator_term.value, out=ratio, - where=denominator_term.value != 0.0, + 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/idl/conftest.py b/fiasco/tests/idl/conftest.py index a8cfabb2..baabcd27 100644 --- a/fiasco/tests/idl/conftest.py +++ b/fiasco/tests/idl/conftest.py @@ -5,12 +5,20 @@ @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) diff --git a/fiasco/tests/idl/helpers.py b/fiasco/tests/idl/helpers.py index 183e0fb5..7a4d4773 100644 --- a/fiasco/tests/idl/helpers.py +++ b/fiasco/tests/idl/helpers.py @@ -23,35 +23,6 @@ ] -def get_hissw_defaults(): - """ - Read default hissw configuration, if available. - """ - try: - import hissw - except ImportError: - return {} - return hissw.defaults - - -def find_chianti_idl_root(idl_codebase_root=None): - """ - Locate the CHIANTI IDL code root. - """ - if idl_codebase_root is not None: - idl_codebase_root = pathlib.Path(idl_codebase_root) - if idl_codebase_root.name != 'idl' and (idl_codebase_root / 'idl').is_dir(): - idl_codebase_root = idl_codebase_root / 'idl' - return idl_codebase_root - - ssw_home = get_hissw_defaults().get('ssw_home') - if ssw_home is None: - return None - - chianti_idl_root = pathlib.Path(ssw_home) / 'packages' / 'chianti' / 'idl' - return chianti_idl_root if chianti_idl_root.is_dir() else None - - def get_idl_test_output_filepath(name, version): if not isinstance(version, Version): version = Version(version) @@ -108,12 +79,7 @@ def setup_idl_environment(ascii_dbase_root, import hissw except ImportError: return None - - hissw_defaults = get_hissw_defaults() - chianti_root = find_chianti_idl_root(idl_codebase_root) - extra_paths = [] - if chianti_root is not None: - extra_paths.extend(d for d, _, _ in os.walk(chianti_root)) + extra_paths = [d for d, _, _ in os.walk(idl_codebase_root)] header = f''' defsysv,'!xuvtop','{ascii_dbase_root}' defsysv,'!abund_file','' @@ -123,19 +89,16 @@ def setup_idl_environment(ascii_dbase_root, extra_paths += [get_pkg_data_path('ssw_gen_functions', package='fiasco.tests.idl')] else: extra_paths += [d for d, _, _ in os.walk(ssw_gen_root)] - idl_only = chianti_root is not None env = hissw.Environment( - ssw_packages=[] if idl_only else ['chianti'], - ssw_home=hissw_defaults.get('ssw_home'), - idl_home=idl_executable or hissw_defaults.get('idl_home'), - idl_only=idl_only, + idl_only=True, + idl_home=idl_executable, header=header, extra_paths=extra_paths, filters={'version_check': version_check}, ) try: _ = env.run('print,!xuvtop', verbose=False) - except (ValueError, hissw.util.SSWIDLError, hissw.util.IDLLicenseError): + except (hissw.util.SSWIDLError, hissw.util.IDLLicenseError): return None else: return env @@ -188,11 +151,10 @@ def run_idl_script(idl_env, # Import here so that this can be used without a hard pytest dependency import pytest pytest.skip("""To run the IDL comparison tests, you must: - 1. Install the hissw package (pip install hissw), - 2. Configure hissw defaults in ~/.hissw/hisswrc or pass - --idl-executable/--idl-codebase-root explicitly, - 3. Have a working IDL and CHIANTI IDL setup. - Without that, you will not be able to generate new IDL comparison test results.""") + 1. Specify a path to a working IDL executable, + 2. Specify a path to the CHIANTI IDL code, + 3. Install the hissw package (pip install hissw) + Without the following, you will not be able to generate new IDL comparison test results.""") # Add versions so they are accessible from within the script if needed input_args = { **input_args, @@ -220,7 +182,7 @@ def run_idl_script(idl_env, return read_idl_test_output(file_name, dbase_version) -def get_chianti_idl_version(idl_codebase_root=None): +def get_chianti_idl_version(idl_codebase_root): """ Get the version of the CHIANTI IDL code being used. @@ -231,11 +193,6 @@ def get_chianti_idl_version(idl_codebase_root=None): read it from the ``VERSION`` file. If that fails, raise a warning and return None. """ - idl_codebase_root = find_chianti_idl_root(idl_codebase_root) - if idl_codebase_root is None: - fiasco.log.warning('Cannot determine CHIANTI IDL version because no CHIANTI IDL root was found.') - return None - # First try to get the version from the current git tag try: import git From efda5ca5defb5fa2518193c5473ecb99c0866852 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 1 May 2026 18:00:01 -0400 Subject: [PATCH 07/12] simplify and fix up tests --- fiasco/tests/idl/test_idl_line_ratio.py | 13 ++--- fiasco/tests/test_fiasco.py | 33 ++++++++++++ fiasco/tests/test_ion.py | 68 ------------------------- 3 files changed, 36 insertions(+), 78 deletions(-) diff --git a/fiasco/tests/idl/test_idl_line_ratio.py b/fiasco/tests/idl/test_idl_line_ratio.py index 01248dbc..da2ea8cf 100644 --- a/fiasco/tests/idl/test_idl_line_ratio.py +++ b/fiasco/tests/idl/test_idl_line_ratio.py @@ -2,7 +2,6 @@ IDL comparison tests for line ratio calculations. """ import astropy.units as u -import numpy as np import pytest import fiasco @@ -37,7 +36,6 @@ def test_eis_fe_xiii_density_ratio_from_idl(idl_env, dbase_version, chianti_idl_ """ formatters = { 'density': lambda x: (10.0**x) * u.cm**(-3), - 'ratio': lambda x: x * u.dimensionless_unscaled, 'temperature': lambda x: (10.0**x) * u.K, 'numerator': lambda x: x * u.angstrom, 'denominator': lambda x: x * u.angstrom, @@ -52,17 +50,12 @@ def test_eis_fe_xiii_density_ratio_from_idl(idl_env, dbase_version, chianti_idl_ chianti_idl_version, format_func=formatters, ) - ion = fiasco.Ion( - 'Fe XIII', - np.atleast_1d(idl_result['temperature']), - hdf5_dbase_root=hdf5_dbase_root, - ) - ratio = fiasco.line_ratio_density( + 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 np.allclose(ratio.value.squeeze(), idl_result['ratio'].value, rtol=2e-3) + 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 4180bdad..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 @@ -35,3 +37,34 @@ def test_proton_electron_ratio(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]) diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index d44195fc..f585a687 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -275,74 +275,6 @@ def test_contribution_function(ion): assert u.allclose(cont_func[0, 0, 0], 2.51408088e-30 * u.cm**3 * u.erg / u.s) -@pytest.mark.requires_dbase_version('>= 8') -def test_line_ratio_density(ion): - density = [1e9, 1e10] * u.cm**(-3) - numerator = 703729.77 * u.angstrom - denominator = 363372.09 * u.angstrom - grouped_numerator = [703729.77, 363372.09] * u.angstrom - grouped_denominator = 151285.93 * u.angstrom - bound_mask = ion.transitions.is_bound_bound - bound_wavelength = ion.transitions.wavelength[bound_mask] - bound_label = ion.transitions.label[bound_mask] - - ratio = fiasco.line_ratio_density(ion, numerator, denominator, density, use_two_ion_model=False) - assert ratio.shape == ion.temperature.shape + density.shape - assert np.isclose(ratio[0, 0].value, 0.040713073461089135) - assert np.isclose(ratio[0, 1].value, 0.04085437059595315) - assert np.isclose(ratio[10, 0].value, 0.04053179876351098) - assert np.isclose(ratio[10, 1].value, 0.040757339413215445) - assert np.isclose(ratio[50, 0].value, 0.0394011731373532) - assert np.isclose(ratio[50, 1].value, 0.04054035992828684) - assert np.isnan(ratio[-1, :]).all() - - grouped_ratio = fiasco.line_ratio_density( - ion, - grouped_numerator, - grouped_denominator, - density, - use_two_ion_model=False, - ) - assert np.isclose(grouped_ratio[0, 0].value, 367497.60772844875) - assert np.isclose(grouped_ratio[0, 1].value, 367496.8140678737) - assert np.isclose(grouped_ratio[10, 1].value, 366355.21439281706) - assert np.isclose(grouped_ratio[50, 0].value, 364554.8767819584) - assert np.isnan(grouped_ratio[-1, :]).all() - - label_ratio = fiasco.line_ratio_density( - ion, - bound_label[np.argmin(np.abs(bound_wavelength - numerator))], - bound_label[np.argmin(np.abs(bound_wavelength - denominator))], - density, - use_two_ion_model=False, - ) - assert u.allclose(label_ratio, ratio, equal_nan=True) - - grouped_label_ratio = fiasco.line_ratio_density( - ion, - bound_label[[np.argmin(np.abs(bound_wavelength - w)) for w in grouped_numerator]], - bound_label[np.argmin(np.abs(bound_wavelength - grouped_denominator))], - density, - use_two_ion_model=False, - ) - assert u.allclose(grouped_label_ratio, grouped_ratio, equal_nan=True) - - density_coupled = 1e15 * u.K * u.cm**(-3) / ion.temperature - ratio_coupled = fiasco.line_ratio_density( - ion, - numerator, - denominator, - density_coupled, - couple_density_to_temperature=True, - use_two_ion_model=False, - ) - assert ratio_coupled.shape == ion.temperature.shape + (1,) - assert np.allclose( - ratio_coupled[[0, 10, 50], 0].value, - [0.04085437059595315, 0.040731805023382285, 0.03821161915021005], - ) - assert np.isnan(ratio_coupled[-1, 0].value) - @pytest.mark.requires_dbase_version('>= 8') @pytest.mark.parametrize(('density', 'use_coupling'), [ ([1e9,] * u.cm**(-3), False), From 1a55f1fca82268d95c91e3570e6c349983d5c2ca Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 1 May 2026 18:00:19 -0400 Subject: [PATCH 08/12] clean up example --- docs/references.bib | 16 +++++ .../user_guide/plot_jet_density_figure.py | 68 +++++-------------- 2 files changed, 33 insertions(+), 51 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 24217dd7..f8aa93b7 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -620,3 +620,19 @@ @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} +} diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index 6c1698ed..343a11dd 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -4,11 +4,7 @@ ======================================== 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 `__ -by Mulay et al. (2017). +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 @@ -19,7 +15,6 @@ import numpy as np from astropy.io import fits -from astropy.utils.data import download_file from astropy.visualization import quantity_support import fiasco @@ -27,44 +22,24 @@ quantity_support() ############################################################################### -# First download the precomputed EIS products. These two FITS files are hosted +# 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 `__. +# 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/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: +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. -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_density( +# 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, @@ -80,19 +55,8 @@ 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', ) +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 @@ -104,13 +68,15 @@ layout='constrained', ) -axes['curve'].plot(density, ratio_curve, color='black', lw=2) +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') +axes['curve'].legend() -observed_image = axes['observed'].imshow(observed_195, origin='lower', cmap='magma', aspect="auto") +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]') From 53b207ebaac50057381abe09a43d4f2feec7efa5 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 1 May 2026 18:05:59 -0400 Subject: [PATCH 09/12] point data url to sunpy --- examples/user_guide/plot_jet_density_figure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index 343a11dd..b9ce65c2 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -26,7 +26,7 @@ # 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/nabobalis/data/main/fiasco' +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: From 1f19a94a44b92269fcb31f8f3b5914f053081491 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Fri, 1 May 2026 18:17:17 -0400 Subject: [PATCH 10/12] fix import; add fe_12 test data --- fiasco/__init__.py | 2 +- fiasco/tests/data/test_file_list.json | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/fiasco/__init__.py b/fiasco/__init__.py index 2bc9068d..be850670 100644 --- a/fiasco/__init__.py +++ b/fiasco/__init__.py @@ -55,7 +55,7 @@ def _get_bibtex(): "list_elements", "list_ions", "get_isoelectronic_sequence", - "line_ratio_density", + "line_ratio", "proton_electron_ratio", "Ion", "Levels", diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 7fc2c83d..47fa18ab 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -771,7 +771,10 @@ "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", @@ -2019,4 +2022,4 @@ "zn_30.rrparams", "zn_31.rrparams" ] -} +} \ No newline at end of file From 7e2d3382d95c9a3084c1b6aae86e4a4a5ef804dd Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Fri, 1 May 2026 22:40:17 -0700 Subject: [PATCH 11/12] Apply suggestion from @nabobalis --- examples/user_guide/plot_jet_density_figure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index b9ce65c2..e67a4bce 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -24,7 +24,7 @@ ############################################################################### # 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/)_. +# 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: From 200897074380880e3a76032c704bfc4ffbea9358 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 4 May 2026 11:32:24 -0400 Subject: [PATCH 12/12] fix citation; add histogram --- docs/references.bib | 6 ++++-- examples/user_guide/plot_jet_density_figure.py | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index f8aa93b7..b87e08af 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -624,7 +624,7 @@ @article{nikolic_suppression_2018 @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, + year = {2017}, month = feb, journal = {Astronomy \& Astrophysics}, volume = {598}, @@ -634,5 +634,7 @@ @article{mulay_Temperature_2017 doi = {10.1051/0004-6361/201628796}, urldate = {2026-05-01}, copyright = {\copyright{} ESO, 2017}, - langid = {english} + langid = {english}, + keywords = {cited}, } + diff --git a/examples/user_guide/plot_jet_density_figure.py b/examples/user_guide/plot_jet_density_figure.py index b9ce65c2..6fdd9901 100644 --- a/examples/user_guide/plot_jet_density_figure.py +++ b/examples/user_guide/plot_jet_density_figure.py @@ -74,6 +74,9 @@ 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")