From c0cd5bf7bf34e881efe8a6d6ba47a30114c9cba5 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Thu, 29 Aug 2024 15:52:41 -0400 Subject: [PATCH] respect T-n coupling in radiative loss calculation --- fiasco/collections.py | 6 +++++- fiasco/tests/test_collections.py | 18 ++++++++++++++++-- 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/fiasco/collections.py b/fiasco/collections.py index 9bc879ef..65698318 100644 --- a/fiasco/collections.py +++ b/fiasco/collections.py @@ -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) diff --git a/fiasco/tests/test_collections.py b/fiasco/tests/test_collections.py index a5a877b3..3a822ecb 100644 --- a/fiasco/tests/test_collections.py +++ b/fiasco/tests/test_collections.py @@ -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) @@ -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): @@ -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 @@ -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 @@ -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)