add 2p method for IonCollection

Author: Jeffrey Reep

fiasco/collections.py

@@ -167,6 +167,26 @@ def free_bound(self, wavelength: u.angstrom, **kwargs):
                 free_bound += fb * abundance * ioneq[:, np.newaxis]
         return free_bound
 
+    @u.quantity_input
+    def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**-3, **kwargs):
+        """
+        """
+        wavelength = np.atleast_1d(wavelength)
+        two_photon = u.Quantity(np.zeros(self.temperature.shape + wavelength.shape),
+                                'erg cm^3 s^-1 Angstrom^-1')
+        for ion in self:
+            if ion.hydrogenic or ion.helium_like:
+                try:
+                    tp = ion.two_photon(wavelength, electron_density, **kwargs)
+                    abundance = ion.abundance
+                    ioneq = ion.ioneq
+                except MissingDatasetException as e:
+                    self.log.warning(f'{ion.ion_name} not included in two-photon emission. {e}')
+                    continue
+                else:
+                    two_photon += tp * abundance * ioneq[:, np.newaxis]
+        return two_photon
+
     @u.quantity_input
     def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength_range=None,
                  bin_width=None, kernel=None, **kwargs):

fiasco/ions.py

@@ -1471,7 +1471,10 @@ def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg,
     def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**(-3)) -> u.erg * u.cm**3 / u.s / u.angstrom:
         """
         """
+        wavelength = np.atleast_1d(wavelength)
+
         prefactor = (const.h * const.c) / (4*np.pi)
+        EM = 1 * u.cm**-5  # CHIANTI assumes a column EM with value of 1 if not specified
 
         if self.hydrogenic:
             A_ji = self._hseq['A']
@@ -1497,8 +1500,8 @@ def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**(-3)) -> u.
         psi_interp = cubic_spline(y_interp)
         energy_dist = (A_ji * y_interp * psi_interp) / (A_sum * wavelength[indices])
 
-        pe_ratio = proton_electron_ratio(self.temperature, **self._instance_kwargs)
-        level_population = self.level_populations(electron_density, **self._instance_kwargs)
+        pe_ratio = proton_electron_ratio(self.temperature)
+        level_population = self.level_populations(electron_density)
 
         # N_j(X+m) = N_j(X+m)/N(X+m) * N(X+m)/N(X) * N(X)/N(H) * N(H)/Ne * Ne
         level_density = level_population[:,0,1] * self.ioneq * self.abundance * pe_ratio * electron_density
@@ -1507,4 +1510,4 @@ def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**(-3)) -> u.
 
         #prefactor * matrix = u.erg / (u.s * u.angstrom * u.cm**2)
 
-        return 0.0 * u.erg * u.cm**3 / u.s / u.angstrom
+        return (prefactor * matrix / EM)