Use the cached proton-electron ratio in Ion.emissivity and make it settable

changelog/471.bugfix.rst

@@ -0,0 +1 @@
+`fiasco.Ion.emissivity` now uses the cached `fiasco.Ion.proton_electron_ratio` property rather than recomputing the proton-to-electron ratio, which walks the entire database, on every call. `fiasco.Ion.proton_electron_ratio` can now also be set directly so that a ratio computed once with `fiasco.proton_electron_ratio` can be shared between ions with the same temperature.

fiasco/ions.py

@@ -70,6 +70,7 @@ def __init__(self,
         super().__init__(ion_name, *args, **kwargs)
         self.temperature = np.atleast_1d(temperature)
         self._dset_names = {}
+        self._proton_electron_ratio = None
         self.abundance = abundance
         self.ionization_fraction = ionization_fraction
         self.ionization_potential = ionization_potential
@@ -203,10 +204,26 @@ def thermal_energy(self) -> u.erg:
         """
         return self.temperature.to('erg', equivalencies=u.equivalencies.temperature_energy())
 
-    @cached_property
+    @property
     @u.quantity_input
     def proton_electron_ratio(self) -> u.dimensionless_unscaled:
-        return proton_electron_ratio(self.temperature, **self._instance_kwargs)
+        """
+        Ratio of proton to electron number density as a function of temperature.
+
+        Computed with `fiasco.proton_electron_ratio` on first access and cached afterward.
+        The ratio depends only on the temperature and the abundance and ionization fraction
+        datasets, so it can also be set directly to avoid recomputing it for every ion that
+        shares the same temperature, e.g.
+        ``ion.proton_electron_ratio = fiasco.proton_electron_ratio(ion.temperature)``.
+        """
+        if self._proton_electron_ratio is None:
+            self._proton_electron_ratio = proton_electron_ratio(self.temperature, **self._instance_kwargs)
+        return self._proton_electron_ratio
+
+    @proton_electron_ratio.setter
+    def proton_electron_ratio(self, value):
+        # Multiplying by np.ones allows for passing in scalar values
+        self._proton_electron_ratio = np.atleast_1d(value) * np.ones(self.temperature.shape)
 
     def next_ion(self):
         """
@@ -1214,8 +1231,7 @@ 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]
+        pe_ratio = self.proton_electron_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)

fiasco/tests/test_ion.py

@@ -582,6 +582,25 @@ def test_ionization_fraction_setter(ion, ioneq_input, ioneq_output):
         assert u.allclose(ion._instance_kwargs['ionization_fraction'], ioneq_input)
 
 
+@pytest.mark.parametrize('value', [
+    0.83,
+    0.83 * np.ones(temperature.shape),
+])
+def test_proton_electron_ratio_setter(ion, value):
+    ion.proton_electron_ratio = value
+    assert ion.proton_electron_ratio.shape == ion.temperature.shape
+    assert u.allclose(ion.proton_electron_ratio, 0.83)
+
+
+@pytest.mark.requires_dbase_version('>= 8')
+def test_emissivity_uses_proton_electron_ratio(ion):
+    # Setting the ratio to 0 should zero the emissivity, which is only the case
+    # if emissivity uses the (cached) property rather than recomputing the ratio.
+    ion.proton_electron_ratio = 0.0
+    emm = ion.emissivity(1e7 * u.cm**-3)
+    assert u.allclose(emm, 0 * u.erg / u.cm**3 / u.s)
+
+
 def test_ionization_fraction_setter_exception(ion):
     # This should fail because the input has len>1 but is not the same
     # shape as the temperature array