Skip to content

Backport PR #480 on branch 0.8 (Fix bug in allowing arrays of density into two-photon continuum calculation) #481

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
wtbarnes merged 1 commit into from
Oct 23, 2026
Merged

Backport PR #480 on branch 0.8 (Fix bug in allowing arrays of density into two-photon continuum calculation) #481

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions fiasco/ions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2010,7 +2010,6 @@ def two_photon(self,
"""
wavelength = np.atleast_1d(wavelength)
electron_density = np.atleast_1d(electron_density)

final_shape = self.temperature.shape + electron_density.shape + wavelength.shape
couple_density_to_temperature = kwargs.setdefault('couple_density_to_temperature', False)
if couple_density_to_temperature:
Expand All @@ -2028,27 +2027,28 @@ def two_photon(self,
else:
return u.Quantity(np.zeros(final_shape), 'erg cm^3 s^-1 Angstrom^-1')
level_index = np.where(self.levels.label==label)
if level_index[0].shape != (1,):
raise ValueError(f'Multiple levels {self.levels.level[level_index]} for label {label}.')
rest_wavelength = self.levels.energy[level_index].to('AA', equivalencies=u.equivalencies.spectral())

# NOTE: Explicitly setting the boundary condition type here to match the behavior of the
# IDL spline interpolation functions. See https://github.com/wtbarnes/fiasco/pull/297 for
# additional details.
cubic_spline = CubicSpline(x_interp, y_interp, bc_type='natural')
x_new = (rest_wavelength / wavelength).decompose().to_value(u.dimensionless_unscaled)
psi_interp = cubic_spline(x_new)
psi_interp = np.where(x_new>1.0, 0.0, psi_interp)

energy_dist = (A_ji * rest_wavelength * psi_interp) / (psi_norm * wavelength**3)
# NOTE: There are known issues when including the proton rates here for some ions so these
# are excluded by default. See https://github.com/wtbarnes/fiasco/pull/260#issuecomment-1955770878
# for more details.
kwargs.setdefault('include_protons', False)
level_population = self.level_populations(electron_density, **kwargs)
level_population = level_population[..., self.levels[level_index].level-1]

# NOTE: This is necessary as the level population grid is not necessarily the same as self.levels.level
# NOTE: Ensure idx is a scalar so that level_population is always 2D
idx = np.where(np.arange(1, self.n_levels+1, 1)==self.levels.level[level_index])[0][0]
level_population = level_population[..., idx]
if couple_density_to_temperature:
electron_density = electron_density[:, np.newaxis]
level_density = level_population / electron_density
matrix = np.outer(level_density, energy_dist).reshape(final_shape)

return const.h * const.c * matrix
21 changes: 17 additions & 4 deletions fiasco/tests/test_ion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -467,13 +467,14 @@ def test_free_bound_radiative_loss(ion, h1):
assert u.isclose(ion.free_bound_radiative_loss()[20], 1.977311794139327e-23 * u.erg * u.cm**3 / u.s)


# The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl
@pytest.mark.requires_dbase_version('>= 8','<= 8.0.7')
@pytest.mark.parametrize(('ion_fixture_name', 'include_protons', 'expected'), [
('c4', False, 0.0 * u.cm**3 * u.erg / u.Angstrom / u.s),
('c5', False, 4.08382025e-25 * u.cm**3 * u.erg / u.Angstrom / u.s),
# The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl
pytest.param(
'c5', False, 4.08382025e-25 * u.cm**3 * u.erg / u.Angstrom / u.s,
marks=pytest.mark.requires_dbase_version('>= 8', '<= 8.0.7')
),
('c6', False, 8.23471268e-26 * u.cm**3 * u.erg / u.Angstrom / u.s),
('c6', True, 6.77829033e-29 * u.cm**3 * u.erg / u.Angstrom / u.s),
])
def test_two_photon(request, ion_fixture_name, include_protons, expected):
ion = request.getfixturevalue(ion_fixture_name)
Expand All @@ -482,6 +483,18 @@ def test_two_photon(request, ion_fixture_name, include_protons, expected):
assert u.allclose(emission[30, 0, 0], expected)


@pytest.mark.parametrize(('density', 'couple_density_to_temperature', 'expected_shape'), [
(1e9*u.cm**(-3), False, (1,)),
([1e9,]*u.cm**(-3), False, (1,)),
([1e8, 1e9, 1e10]*u.cm**(-3), False, (3,)),
(1e9*u.cm**(-3)*np.ones((100,)), True, (1,)),
])
def test_two_photon_density_shape(c5, density, couple_density_to_temperature, expected_shape):
wavelength = np.linspace(20, 100, 5) * u.AA
emission = c5.two_photon(wavelength, density, couple_density_to_temperature=couple_density_to_temperature)
assert emission.shape == c5.temperature.shape+expected_shape+wavelength.shape


def test_free_bound_no_recombining(h1):
# This is test the case where there is no data available for the recombining
# ion (H 2)
Expand Down
Loading