Drop assumed conversion factor between H and electron density (#273)

examples/user_guide/aia_response.py

@@ -45,7 +45,7 @@
 ############################################################
 # Compute the contribution function,
 #
-# .. math:: G(n,T,\lambda) = 0.83\mathrm{Ab}\frac{hc}{\lambda}N_{\lambda}A_{\lambda}f\frac{1}{n}
+# .. math:: G(n,T,\lambda) = \mathrm{Ab}\frac{hc}{\lambda}N_{\lambda}A_{\lambda}f\frac{1}{n}
 #
 # for each transition of Fe 18 at constant pressure of :math:`10^{15}`
 # K :math:`\mathrm{cm}^{-3}`. Note that we use the ``couple_density_to_temperature``
@@ -73,13 +73,10 @@
 # .. math:: K_c(T) = \int_0^{\infty}\mathrm{d}\lambda\,G(\lambda,T)R_c(\lambda)
 #
 # We divide by :math:`hc/\lambda` in order to
-# convert from units of energy to photons. The factor of :math:`0.83` is a
-# relative scaling factor for the abundance of H and is not included in the
-# temperature responses computed by
-# `aia_get_response.pro <https://sohoftp.nascom.nasa.gov/solarsoft/sdo/aia/idl/response/aia_get_response.pro>`_.
+# convert from units of energy to photons.
 # For more more information on the AIA wavelength response calculation,
 # see :cite:t:`boerner_initial_2012`.
-K = (g / energy * response_transitions).sum(axis=2) / (4*np.pi*u.steradian) / 0.83
+K = (g / energy * response_transitions).sum(axis=2) / (4*np.pi*u.steradian)
 K = K.squeeze().to('cm5 ct pix-1 s-1')
 
 ############################################################

fiasco/collections.py

@@ -84,7 +84,7 @@ def free_free(self, wavelength: u.angstrom):
         r"""
         Compute combined free-free continuum emission (bremsstrahlung).
 
-        .. note:: Both abundance and ionization equilibrium are included here
+        .. note:: Both abundance and ionization fraction are included here
 
         The combined free-free continuum is given by,
 
@@ -126,7 +126,7 @@ def free_bound(self, wavelength: u.angstrom, **kwargs):
         r"""
         Compute combined free-bound continuum emission.
 
-        .. note:: Both abundance and ionization equilibrium are included here.
+        .. note:: Both abundance and ionization fraction are included here.
 
         The combined free-bound continuum is given by,
 

fiasco/ions.py

@@ -191,11 +191,11 @@ def ioneq(self):
         ionization equilibrium outside of this temperature range, it is better to use the ionization
         and recombination rates.
 
-        Note
-        ----
-        The cubic interpolation is performed in log-log spaceusing a Piecewise Cubic Hermite
-        Interpolating Polynomial with `~scipy.interpolate.PchipInterpolator`. This helps to
-        ensure smoothness while reducing oscillations in the interpolated ionization fractions.
+        .. note::
+
+            The cubic interpolation is performed in log-log space using a Piecewise Cubic Hermite
+            Interpolating Polynomial with `~scipy.interpolate.PchipInterpolator`. This helps to
+            ensure smoothness while reducing oscillations in the interpolated ionization fractions.
 
         See Also
         --------
@@ -709,7 +709,7 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er
 
         .. math::
 
-           G_{ij} = \frac{n_H}{n_e}\mathrm{Ab}(X)f_{X,k}N_jA_{ij}\Delta E_{ij}\frac{1}{n_e},
+           G_{ij} = mathrm{Ab}(X)f_{X,k}N_jA_{ij}\Delta E_{ij}\frac{1}{n_e},
 
         Note that the contribution function is often defined in differing ways by different authors.
         The contribution function is defined as above in :cite:t:`young_chianti_2016`.
@@ -720,6 +720,13 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er
 
            ion.transitions.wavelength[~ion.transitions.is_twophoton]
 
+        .. important::
+
+            The ratio :math:`n_H/n_e`, which is often approximated as :math:`n_H/n_e\approx0.83`, is
+            explicitly not included here. This means that when computing an intensity with the
+            result of this function, the accompanying emission measure is
+            :math:`\mathrm{EM}=\mathrm{d}hn_Hn_e` rather than :math:`n_e^2`.
+
         Parameters
         ----------
         density: `~astropy.units.Quantity`
@@ -753,7 +760,7 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er
         else:
             term = np.outer(self.ioneq, 1./density)
             term = term[:, :, np.newaxis]
-        term *= self.abundance * 0.83
+        term *= self.abundance
         # Exclude two-photon transitions
         upper_level = self.transitions.upper_level[~self.transitions.is_twophoton]
         wavelength = self.transitions.wavelength[~self.transitions.is_twophoton]
@@ -772,19 +779,24 @@ def emissivity(self, density: u.cm**(-3), **kwargs) -> u.erg * u.cm**(-3) / u.s:
 
         .. math::
 
-           \epsilon(n_e,T) = G(n_e,T)n_e^2
+           \epsilon(n_e,T) = G(n_e,T)n_Hn_e
 
-        where :math:`G` is the contribution function, :math:`n_e` is the electron
-        density, and :math:`T` is the temperature.
+        where :math:`G` is the contribution function, :math:`n_H` is the H (or proton) density,
+        :math:`n_e` is the electron density, and :math:`T` is the temperature.
         Note that, like the contribution function, emissivity is often defined in
         in differing ways by different authors.
         Here, we use the definition of the emissivity as given by Eq. 3 of
         :cite:t:`young_chianti_2016`.
 
+        .. note::
+
+            The H number density, :math:`n_H`, is computed using ``density`` combined with the
+            output of `~fiasco.proton_electron_ratio`.
+
         Parameters
         ----------
         density : `~astropy.units.Quantity`
-            Electron number density
+            Electron number density.
         couple_density_to_temperature: `bool`, optional
             If True, the density will vary along the same axis as temperature
             in the computed level populations. The number of densities must be the same as the number of temperatures. This is useful, for
@@ -808,14 +820,16 @@ def emissivity(self, density: u.cm**(-3), **kwargs) -> u.erg * u.cm**(-3) / u.s:
         contribution_function : Calculate contribution function, :math:`G(n,T)`
         """
         density = np.atleast_1d(density)
+        pe_ratio = proton_electron_ratio(self.temperature, **self._instance_kwargs)
+        pe_ratio = pe_ratio[:, np.newaxis, np.newaxis]
         g = self.contribution_function(density, **kwargs)
         density_squared = density**2
         couple_density_to_temperature = kwargs.get('couple_density_to_temperature', False)
         if couple_density_to_temperature:
             density_squared = density_squared[:, np.newaxis, np.newaxis]
         else:
             density_squared = density_squared[np.newaxis, :, np.newaxis]
-        return g * density_squared
+        return g * pe_ratio * density_squared
 
     @u.quantity_input
     def intensity(self,
@@ -828,7 +842,7 @@ def intensity(self,
 
         .. math::
 
-           I = \frac{1}{4\pi}\int\mathrm{d}T,G(n,T)n^2\frac{dh}{dT}
+           I = \frac{1}{4\pi}\int\mathrm{d}T,G(n,T)n_Hn_e\frac{dh}{dT}
 
         which, in the isothermal approximation, can be simplified to,
 
@@ -840,7 +854,7 @@ def intensity(self,
 
         .. math::
 
-           \mathrm{EM}(T) = \int\mathrm{d}h\,n^2
+           \mathrm{EM}(T) = \int\mathrm{d}h\,n_Hn_e
 
         is the column emission measure.
 
@@ -850,7 +864,8 @@ def intensity(self,
             Electron number density
         emission_measure : `~astropy.units.Quantity`
             Column emission measure. Must be either a scalar, an array of length 1, or
-            an array with the same length as ``temperature``.
+            an array with the same length as ``temperature``. Note that it is assumed
+            that the emission measure is the product of the H and electron density.
         couple_density_to_temperature: `bool`, optional
             If True, the density will vary along the same axis as temperature.
             The number of densities must be the same as the number of temperatures.
@@ -1216,7 +1231,7 @@ def recombination_rate(self) -> u.cm**3 / u.s:
 
             \alpha_{R} = \alpha_{RR} + \alpha_{DR}
 
-        .. warning::
+        .. important::
 
             For most ions, this total recombination rate is computed by summing the
             outputs of the `radiative_recombination_rate` and `dielectronic_recombination_rate` methods.
@@ -1261,6 +1276,10 @@ def free_free(self, wavelength: u.angstrom) -> u.erg * u.cm**3 / u.s / u.angstro
         r"""
         Free-free continuum emission as a function of temperature and wavelength.
 
+        .. note::
+
+            This does not include ionization fraction or abundance factors.
+
         Free-free emission, also known as *bremsstrahlung* (or “braking radiation”),
         is produced when an ion interacts with a free electron, reduces the momentum
         of the free electron, and, by conservation of energy and momentum, produces
@@ -1382,10 +1401,11 @@ def free_bound(self,
         r"""
         Free-bound continuum emission of the recombined ion.
 
-        .. note:: Does not include ionization equilibrium or abundance.
-                  Unlike the equivalent IDL routine, the output here is not
-                  expressed per steradian and as such the factor of
-                  :math:`1/4\pi` is not included.
+        .. note::
+
+            This does not include the ionization fraction or abundance factors.
+            Unlike the equivalent IDL routine, the output here is not
+            expressed per steradian and as such the factor of :math:`1/4\pi` is not included.
 
         When an electron is captured by an ion of charge :math:`z+1`
         (the recombining ion), it creates a an ion of charge :math:`z`

fiasco/tests/idl/test_idl_goft.py

@@ -76,7 +76,8 @@ def test_idl_compare_goft(idl_env, hdf5_dbase_root, dbase_version, ion_name, wav
                      ioneq_filename=idl_result['ioneq'])
     contribution_func = ion.contribution_function(idl_result['density'])
     idx = np.argmin(np.abs(ion.transitions.wavelength[~ion.transitions.is_twophoton] - idl_result['wavelength']))
-    goft_python = contribution_func[:, 0, idx]
+    # NOTE: Multiply by 0.83 because the fiasco calculation does not include the n_H/n_e ratio
+    goft_python = contribution_func[:, 0, idx] * 0.83
     # Find relevant range for comparison. The solutions may diverge many orders of magnitude below
     # the peak of the contribution function, but that is not relevant in assessing meaningful
     # differences between the IDL and Python approaches. Thus, we only assess the differences

fiasco/tests/test_collections.py

@@ -107,10 +107,10 @@ def test_free_bound(another_collection, wavelength):
     assert u.allclose(fb[index_t, index_w], 3.057781475607237e-36 * u.Unit('erg cm3 s-1 Angstrom-1'))
 
 
-def test_radiative_los(collection):
+def test_radiative_loss(collection):
     rl = collection.radiative_loss(1e9*u.cm**(-3))
     # This value has not been checked for correctness
-    assert u.allclose(rl[0,0], 3.2389535764824023e-24*u.Unit('erg cm3 s-1'))
+    assert u.allclose(rl[0,0], 3.90235371e-24*u.Unit('erg cm3 s-1'))
 
 
 def test_spectrum(hdf5_dbase_root):

fiasco/tests/test_ion.py

@@ -194,7 +194,7 @@ def test_contribution_function(ion):
     cont_func = ion.contribution_function(1e7 * u.cm**-3)
     assert cont_func.shape == ion.temperature.shape + (1, ) + ion._wgfa['wavelength'].shape
     # This value has not been tested for correctness
-    assert u.allclose(cont_func[0, 0, 0], 2.08668713e-30 * u.cm**3 * u.erg / u.s)
+    assert u.allclose(cont_func[0, 0, 0], 2.51408088e-30 * u.cm**3 * u.erg / u.s)
 
 @pytest.mark.parametrize(('density', 'use_coupling'), [
     ([1e9,] * u.cm**(-3), False),
@@ -264,7 +264,7 @@ def test_emissivity(ion):
     emm = ion.emissivity(1e7 * u.cm**-3)
     assert emm.shape == ion.temperature.shape + (1, ) + ion._wgfa['wavelength'].shape
     # This value has not been tested for correctness
-    assert u.allclose(emm[0, 0, 0], 2.08668713e-16 * u.erg / u.cm**3 / u.s)
+    assert u.allclose(emm[0, 0, 0], 2.18000422e-16 * u.erg / u.cm**3 / u.s)
 
 
 @pytest.mark.parametrize('em', [