update the logic of Itoh integrated f-f gf

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,22 +315,16 @@ 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
         -------
         rad_loss_total : `~astropy.units.Quantity`
             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:

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]))])

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')

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)

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)

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)