Skip to content
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.

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

Respect temperature-density coupling in radiative loss calculation #295

Merged
merged 1 commit into from
Aug 30, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
6 changes: 5 additions & 1 deletion fiasco/collections.py
Original file line number Diff line number Diff line change
Expand Up @@ -348,7 +348,11 @@ def bound_bound_radiative_loss(self, density, **kwargs) -> u.Unit('erg cm3 s-1')
The bolometric bound-bound radiative loss rate per unit emission measure
"""
density = np.atleast_1d(density)
rad_loss = u.Quantity(np.zeros(self.temperature.shape + density.shape), 'erg cm^3 s^-1')
if kwargs.get('couple_density_to_temperature', False):
shape = self.temperature.shape + (1,)
else:
shape = self.temperature.shape + density.shape
rad_loss = u.Quantity(np.zeros(shape), 'erg cm^3 s^-1')
for ion in self:
try:
g = ion.contribution_function(density, **kwargs)
Expand Down
18 changes: 16 additions & 2 deletions fiasco/tests/test_collections.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ def test_contains(collection, hdf5_dbase_root):
def test_length(collection):
assert len(collection) == len(collection._ion_list)


@pytest.mark.parametrize('wavelength', [wavelength, wavelength[50]])
def test_free_free(another_collection, wavelength):
ff = another_collection.free_free(wavelength)
Expand All @@ -106,6 +107,7 @@ def test_free_bound(another_collection, wavelength):
index_t = 24 # This is approximately where the ioneq for Fe V peaks
assert u.allclose(fb[index_t, index_w], 3.057781475607237e-36 * u.Unit('erg cm3 s-1 Angstrom-1'))


@pytest.mark.requires_dbase_version('>= 8')
@pytest.mark.parametrize('wavelength', [wavelength, wavelength[450]])
def test_two_photon(collection, wavelength, hdf5_dbase_root):
Expand All @@ -121,6 +123,7 @@ def test_two_photon(collection, wavelength, hdf5_dbase_root):
# This value has not been checked for correctness
assert u.allclose(tp[index_w, index_t, 0], 3.48580645e-27 * u.Unit('erg cm3 s-1 Angstrom-1'))


@pytest.mark.requires_dbase_version('>= 8')
def test_radiative_loss(collection, hdf5_dbase_root):
# add Li III to the test to include an ion that throws a MissingDatasetException
Expand All @@ -131,6 +134,7 @@ def test_radiative_loss(collection, hdf5_dbase_root):
# These values have not been checked for correctness
u.allclose(rl[0], [3.90235371e-24, 4.06540902e-24, 4.08411295e-24] * u.erg * u.cm**3 / u.s)


@pytest.mark.requires_dbase_version('>= 8')
def test_radiative_loss_bound_bound(collection, hdf5_dbase_root):
# add Li III to the test to include an ion that throws a MissingDatasetException
Expand All @@ -141,20 +145,30 @@ def test_radiative_loss_bound_bound(collection, hdf5_dbase_root):
# These values have not been checked for correctness
u.allclose(rl[0], [3.90235371e-24, 4.06540902e-24, 4.08411295e-24] * u.erg * u.cm**3 / u.s)


@pytest.mark.requires_dbase_version('>= 8')
def test_radiative_loss_free_free(collection, hdf5_dbase_root):
def test_radiative_loss_free_free(collection):
rl = collection.free_free_radiative_loss()
assert rl.shape == collection.temperature.shape
# This value has not been checked for correctness
u.isclose(rl[0], 2.72706455e-35 * u.erg * u.cm**3 / u.s)


@pytest.mark.requires_dbase_version('>= 8')
def test_radiative_loss_free_bound(collection, hdf5_dbase_root):
def test_radiative_loss_free_bound(collection):
rl = collection.free_bound_radiative_loss()
assert rl.shape == collection.temperature.shape
# This value has not been checked for correctness
u.isclose(rl[0], 1.13808317e-33 * u.erg * u.cm**3 / u.s)


def test_radiative_loss_temperature_density_coupling(another_collection):
p0 = 1e15 * u.K * u.cm**(-3)
density = p0 / another_collection.temperature
rl = another_collection.radiative_loss(density, couple_density_to_temperature=True)
assert rl.shape == another_collection.temperature.shape + (1,)


@pytest.mark.requires_dbase_version('>= 8')
def test_spectrum(hdf5_dbase_root):
i1 = fiasco.Ion('H 1', 1 * u.MK, hdf5_dbase_root=hdf5_dbase_root)
Expand Down