diff --git a/fiasco/collections.py b/fiasco/collections.py index f79699da..923222e2 100644 --- a/fiasco/collections.py +++ b/fiasco/collections.py @@ -303,7 +303,7 @@ def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength return wavelength, spectrum @u.quantity_input - def radiative_loss(self, density: u.cm**(-3), **kwargs) -> u.Unit('erg cm3 s-1'): + def radiative_loss(self, density: u.cm**(-3), use_itoh=False, **kwargs) -> u.Unit('erg cm3 s-1'): r""" Calculate the total wavelength-integrated radiative loss rate including the bound-bound, free-bound, and free-free emission contributions @@ -315,14 +315,8 @@ def radiative_loss(self, density: u.cm**(-3), **kwargs) -> u.Unit('erg cm3 s-1') ---------- density : `~astropy.units.Quantity` Electron number density - itoh : `bool`, optional - Specify whether to use the approximations specified by :cite:t:`itoh_radiative_2002` to - calculate the wavelength-integrated free-free Gaunt factor. - If true, use the forms by :cite:t:`itoh_radiative_2002`. If false (default), use the forms by - :cite:t:`sutherland_accurate_1998`. - relativistic : `bool`, optional - If using the :cite:t:`itoh_radiative_2002` approximations, use the relativistic form - instead of the non-relativistic form. + use_itoh : `bool`, optional + Whether to use the Itoh Gaunt Factors for free-free. Defaults to false. Returns ------- @@ -330,7 +324,7 @@ def radiative_loss(self, density: u.cm**(-3), **kwargs) -> u.Unit('erg cm3 s-1') The total bolometric radiative loss rate """ rad_loss_bound_bound = self.bound_bound_radiative_loss(density, **kwargs) - rad_loss_free_free = self.free_free_radiative_loss() + rad_loss_free_free = self.free_free_radiative_loss(use_itoh) rad_loss_free_bound = self.free_bound_radiative_loss() rad_loss_total = (rad_loss_bound_bound @@ -371,21 +365,15 @@ def bound_bound_radiative_loss(self, density, **kwargs) -> u.Unit('erg cm3 s-1') return rad_loss @u.quantity_input - def free_free_radiative_loss(self, itoh=False, relativistic=True) -> u.Unit('erg cm3 s-1'): + def free_free_radiative_loss(self, use_itoh=False) -> u.Unit('erg cm3 s-1'): r""" Calculate the radiative loss rate from free-free emission (bremsstrahlung) integrated over wavelength. Parameters ---------- - itoh : `bool`, optional - Specify whether to use the approximations specified by :cite:t:`itoh_radiative_2002` to - calculate the wavelength-integrated free-free Gaunt factor. - If true, use the forms by :cite:t:`itoh_radiative_2002`. If false (default), use the forms by - :cite:t:`sutherland_accurate_1998`. - relativistic : `bool`, optional - If using the :cite:t:`itoh_radiative_2002` approximations, use the relativistic form - instead of the non-relativistic form. + use_itoh : `bool`, optional + Whether to use the Itoh Gaunt Factors. Defaults to false. Returns ------- @@ -395,7 +383,7 @@ def free_free_radiative_loss(self, itoh=False, relativistic=True) -> u.Unit('erg free_free = u.Quantity(np.zeros(self.temperature.shape), 'erg cm^3 s^-1') for ion in self: try: - ff = ion.free_free_radiative_loss(itoh, relativistic) + ff = ion.free_free_radiative_loss(use_itoh) abundance = ion.abundance ioneq = ion.ioneq except MissingDatasetException as e: diff --git a/fiasco/gaunt.py b/fiasco/gaunt.py index 0e701a6a..b541bede 100644 --- a/fiasco/gaunt.py +++ b/fiasco/gaunt.py @@ -154,7 +154,7 @@ def _free_free_sutherland(self, temperature: u.K, charge_state, wavelength: u.an return u.Quantity(np.where(gf < 0., 0., gf)) @u.quantity_input - def free_free_integrated(self, temperature: u.K, charge_state, itoh=False, relativistic=True) -> u.dimensionless_unscaled: + def free_free_integrated(self, temperature: u.K, charge_state, use_itoh=False) -> u.dimensionless_unscaled: """ The wavelength-integrated free-free Gaunt factor, used for calculating the total radiative losses from free-free emission. @@ -165,24 +165,18 @@ def free_free_integrated(self, temperature: u.K, charge_state, itoh=False, relat The temperature(s) for which to calculate the Gaunt factor charge_state : `int`, The charge state of the ion - itoh : `bool`, optional - Specify whether to use the approximations specified by :cite:t:`itoh_radiative_2002`. - If true, use the forms by :cite:t:`itoh_radiative_2002`. If false (default), use the forms by - :cite:t:`sutherland_accurate_1998`. - relativistic : `bool`, optional - If using the :cite:t:`itoh_radiative_2002` approximations, use the relativistic form - instead of the non-relativistic form. + use_itoh : `bool`, optional + Whether to use the Itoh Gaunt Factors. Defaults to false. """ if charge_state == 0: return u.Quantity(np.zeros(temperature.shape)) else: - if itoh and relativistic: - return self._free_free_itoh_integrated_relativistic(temperature, charge_state) - elif itoh and not relativistic: - return self._free_free_itoh_integrated_nonrelativistic(temperature, charge_state) - else: - return self._free_free_sutherland_integrated(temperature, charge_state) - + gf = self._free_free_sutherland_integrated(temperature, charge_state) + if use_itoh: + gf_itoh = self._free_free_itoh_integrated(temperature, charge_state) + gf = np.where(np.isnan(gf_itoh), gf, gf_itoh) + return gf + @needs_dataset('gffint') @u.quantity_input def _free_free_sutherland_integrated(self, temperature: u.K, charge_state) -> u.dimensionless_unscaled: @@ -208,6 +202,28 @@ def _free_free_sutherland_integrated(self, temperature: u.K, charge_state) -> u. # The spline fit was pre-calculated by Sutherland 1998: return self._gffint['gaunt_factor'][index] + delta * (self._gffint['s1'][index] + delta * (self._gffint['s2'][index] + delta * self._gffint['s3'][index])) + @u.quantity_input + def _free_free_itoh_integrated(self, temperature: u.K, charge_state) -> u.dimensionless_unscaled: + r""" + The wavelength-integrated free-free Gaunt factor, as specified by :cite:t:`itoh_radiative_2002`. + + Parameters + ---------- + temperature : `~astropy.units.Quantity` + The temperature(s) for which to calculate the Gaunt factor + charge_state : `int`, + The charge state of the ion + """ + temperature = np.atleast_1d(temperature) + try: + gf_relativistic = self._free_free_itoh_integrated_relativistic(temperature, charge_state) + except MissingDatasetException: + gf_relativistic = np.full(len(temperature), np.nan) + try: + gf_nonrelativistic = self._free_free_itoh_integrated_nonrelativistic(temperature, charge_state) + except MissingDatasetException: + gf_nonrelativistic = np.full(len(temperature), np.nan) + return np.where(np.isnan(gf_nonrelativistic), gf_relativistic, gf_nonrelativistic) @needs_dataset('itoh_integrated_gaunt_nonrel') @u.quantity_input @@ -233,7 +249,7 @@ def _free_free_itoh_integrated_nonrelativistic(self, temperature: u.K, charge_st Gamma = (np.log10(gamma_squared) + 0.5) / 2.5 for j in range(len(summation)): if np.log10(gamma_squared[j]) < -3.0 or np.log10(gamma_squared[j]) > 2.0: - summation[j] = self._free_free_sutherland_integrated(temperature[j], charge_state) + summation[j] = np.nan else: b_array = self._itoh_integrated_gaunt_nonrel['b_i'] G_array = np.array([Gamma[j]**i for i in range(len(b_array))]) @@ -263,10 +279,8 @@ def _free_free_itoh_integrated_relativistic(self, temperature: u.K, charge_state t = (np.log10(temperature.data)-7.25)/1.25 summation = u.Quantity(np.zeros(temperature.shape)) for j in range(len(summation)): - if np.log10(temperature[j].data) < 6.0: - summation[j] = self._free_free_itoh_integrated_nonrelativistic(temperature[j], charge_state) - elif np.log10(temperature[j].data) > 8.5: - summation[j] = self._free_free_sutherland_integrated(temperature[j], charge_state) + if np.log10(temperature[j].data) < 6.0 or np.log10(temperature[j].data) > 8.5: + summation[j] = np.nan else: a_matrix = self._itoh_integrated_gaunt['a_ik'] z_array = np.array([z**i for i in range(len(a_matrix[:,0]))]) diff --git a/fiasco/ions.py b/fiasco/ions.py index af26d0c7..72056364 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -1318,7 +1318,7 @@ def free_free(self, wavelength: u.angstrom) -> u.erg * u.cm**3 / u.s / u.angstro / np.sqrt(self.temperature)[:, np.newaxis]) @u.quantity_input - def free_free_radiative_loss(self, itoh=False, relativistic=True) -> u.erg * u.cm**3 / u.s: + def free_free_radiative_loss(self, use_itoh=False) -> u.erg * u.cm**3 / u.s: r""" Free-free continuum radiative losses as a function of temperature. @@ -1339,20 +1339,14 @@ def free_free_radiative_loss(self, itoh=False, relativistic=True) -> u.erg * u.c wavelength-integrated free-free Gaunt factor. The prefactor :math:`F_{k}` is defined in Equation 19 of :cite:t:`sutherland_accurate_1998`, with a value of :math:`F_k\approx1.42555669\times10^{-27}\,\mathrm{cm}^{5}\,\mathrm{g}\,\mathrm{K}^{-1/2}\,\mathrm{s}^{3}`. - + Parameters ---------- - itoh : `bool`, optional - Specify whether to use the approximations specified by :cite:t:`itoh_radiative_2002` to - calculate the wavelength-integrated free-free Gaunt factor. - If true, use the forms by :cite:t:`itoh_radiative_2002`. If false (default), use the forms by - :cite:t:`sutherland_accurate_1998`. - relativistic : `bool`, optional - If using the :cite:t:`itoh_radiative_2002` approximations, use the relativistic form - instead of the non-relativistic form. + use_itoh : `bool`, optional + Whether to use the Itoh Gaunt Factors. Defaults to false. """ prefactor = (16./3**1.5) * np.sqrt(2. * np.pi * const.k_B/(const.hbar**2 * const.m_e**3)) * (const.e.esu**6 / const.c**3) - gf = self.gaunt_factor.free_free_integrated(self.temperature, self.charge_state, itoh, relativistic) + gf = self.gaunt_factor.free_free_integrated(self.temperature, self.charge_state, use_itoh) return (prefactor * self.charge_state**2 * gf * np.sqrt(self.temperature)) @needs_dataset('fblvl', 'ip') diff --git a/fiasco/tests/test_collections.py b/fiasco/tests/test_collections.py index 3344dcc0..f3e6ab07 100644 --- a/fiasco/tests/test_collections.py +++ b/fiasco/tests/test_collections.py @@ -142,19 +142,25 @@ def test_radiative_loss_bound_bound(collection, hdf5_dbase_root): 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): +@pytest.mark.parametrize(('index','expected'), + [(0, 2.72706455e-35), + (75, 5.59153955e-31),]) +def test_radiative_loss_free_free(collection, index, expected): 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) + u.isclose(rl[index], expected * u.erg * u.cm**3 / u.s) @pytest.mark.requires_dbase_version('>= 9.0.1') -@pytest.mark.parametrize(('itoh','relativistic','index','expected'), - [(True, True, 0, 2.71800458e-35), (True, False, 0, 2.71800458e-35), - (False, False, 0, 2.72706455e-35), (True, True, 75, 5.57512083e-31), - (True, False, 75, 5.57346003e-31), (False, False, 75, 5.59153955e-31)]) -def test_radiative_loss_free_free_itoh(collection, itoh, relativistic, index, expected): - rl = collection.free_free_radiative_loss(itoh, relativistic) +@pytest.mark.parametrize(('index','expected'), + [(0, 2.71800458e-35), + (75, 5.57346003e-31),]) +def test_radiative_loss_free_free_itoh(collection, index, expected): + """ + For database versions >= 9.0.1, the Itoh Gaunt factors give a different + result for the free-free radiative loss. + """ + rl = collection.free_free_radiative_loss(use_itoh=True) assert rl.shape == collection.temperature.shape # This value has not been checked for correctness u.isclose(rl[index], expected * u.erg * u.cm**3 / u.s) diff --git a/fiasco/tests/test_gaunt.py b/fiasco/tests/test_gaunt.py index b74426d0..d44a5a49 100644 --- a/fiasco/tests/test_gaunt.py +++ b/fiasco/tests/test_gaunt.py @@ -37,12 +37,13 @@ def test_gaunt_factor_free_free_integrated(gaunt_factor, charge_state, expected) assert u.allclose(gf[0], expected * u.dimensionless_unscaled) @pytest.mark.requires_dbase_version('>= 9.0.1') -@pytest.mark.parametrize(('charge_state','itoh','relativistic','index','expected'), - [(3, True, True, 0, 1.2673757), (3, True, False, 0, 1.2673757), - (3, False, False, 0, 1.27165875), (3, True, True, 75, 1.3466916), - (3, True, False, 75, 1.34662286), (3, False, False, 75, 1.35061768)]) -def test_gaunt_factor_free_free_integrated_itoh(gaunt_factor, charge_state, itoh, relativistic, index, expected): - gf = gaunt_factor.free_free_integrated(temperature, charge_state, itoh, relativistic) +@pytest.mark.parametrize(('charge_state','index','expected'), + [(1, 0, 1.40965757), (1, 75, 1.20459895), + (2, 0, 1.32203268), (2, 75, 1.28680539), + (3, 0, 1.2673757), (3, 75, 1.34662286), + ]) +def test_gaunt_factor_free_free_integrated_itoh(gaunt_factor, charge_state, index, expected): + gf = gaunt_factor.free_free_integrated(temperature, charge_state, use_itoh=True) assert gf.shape == temperature.shape # This value has not been tested for correctness assert u.allclose(gf[index], expected * u.dimensionless_unscaled) diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index f837b1a1..20a8637c 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -339,23 +339,28 @@ def test_free_free(ion): # This value has not been tested for correctness assert u.allclose(emission[0], 1.72804216e-29 * u.cm**3 * u.erg / u.Angstrom / u.s) - -def test_free_free_radiative_loss(h1, fe20): +@pytest.mark.parametrize(('index','expected'), + [(0, 1.79093013e-22), + (75, 3.06577525e-21), ]) +def test_free_free_radiative_loss(h1, fe20, index, expected): assert fe20.free_free_radiative_loss().shape == fe20.temperature.shape assert u.allclose(h1.free_free_radiative_loss(), 0.0 * u.erg * u.cm**3 / u.s) # This value has not been tested for correctness - assert u.isclose(fe20.free_free_radiative_loss()[0], 1.79093013e-22 * u.erg * u.cm**3 / u.s) + assert u.isclose(fe20.free_free_radiative_loss()[index], expected * u.erg * u.cm**3 / u.s) @pytest.mark.requires_dbase_version('>= 9.0.1') -@pytest.mark.parametrize(('itoh','relativistic','index','expected'), - [(True, True, 0, 1.79093013e-22), (True, False, 0, 1.79093013e-22), - (False, False, 0, 1.79093013e-22), (True, True, 75, 3.05586194e-21), - (True, False, 75, 3.05628894e-21), (False, False, 75, 3.06577525e-21)]) -def test_free_free_radiative_loss_itoh(h1, fe20, itoh, relativistic, index, expected): - assert fe20.free_free_radiative_loss(itoh, relativistic).shape == fe20.temperature.shape - assert u.allclose(h1.free_free_radiative_loss(), 0.0 * u.erg * u.cm**3 / u.s) +@pytest.mark.parametrize(('index','expected'), + [(0, 1.79093013e-22), + (75, 3.05628894e-21), ]) +def test_free_free_radiative_loss_itoh(h1, fe20, index, expected): + """ + For database versions >= 9.0.1, the Itoh Gaunt factors give a different + result for the free-free radiative loss at temperature > 1 MK. + """ + assert fe20.free_free_radiative_loss(use_itoh=True).shape == fe20.temperature.shape + assert u.allclose(h1.free_free_radiative_loss(use_itoh=True), 0.0 * u.erg * u.cm**3 / u.s) # This value has not been tested for correctness - assert u.isclose(fe20.free_free_radiative_loss(itoh, relativistic)[index], expected * u.erg * u.cm**3 / u.s) + assert u.isclose(fe20.free_free_radiative_loss(use_itoh=True)[index], expected * u.erg * u.cm**3 / u.s) def test_free_bound(ion): emission = ion.free_bound(200 * u.Angstrom)