From 7933c3cdd4b1fc45fb62a99155e30bc55a0d3dfe Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 16 Sep 2024 11:17:00 -1000 Subject: [PATCH] move free-bound and klgfb to GauntFactor --- fiasco/base.py | 5 ---- fiasco/gaunt.py | 76 +++++++++++++++++++++++++++++++++++++++++++++++-- fiasco/ions.py | 14 ++------- 3 files changed, 76 insertions(+), 19 deletions(-) diff --git a/fiasco/base.py b/fiasco/base.py index 6aabe770..62c36dcf 100644 --- a/fiasco/base.py +++ b/fiasco/base.py @@ -92,11 +92,6 @@ class ContinuumBase(Base): .. note:: This is not meant to be instantiated directly by the user and primarily serves as a base class for `~fiasco.Ion`. """ - @property - def _klgfb(self): - data_path = '/'.join(['continuum', 'klgfb']) - return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path) - @property def _verner(self): data_path = '/'.join([self.atomic_symbol.lower(), self._ion_name, 'continuum', diff --git a/fiasco/gaunt.py b/fiasco/gaunt.py index 03d4835a..523afdb7 100644 --- a/fiasco/gaunt.py +++ b/fiasco/gaunt.py @@ -5,6 +5,7 @@ import astropy.units as u import numpy as np +from scipy.interpolate import splev, splrep from scipy.ndimage import map_coordinates from scipy.special import polygamma @@ -18,6 +19,10 @@ class GauntFactor: """ Class for calculating the Gaunt factor for various continuum processes. + + The Gaunt factor is defined as the ratio of the true cross-section to the + semi-classical Kramers cross-section, and thus is essentially a multiplicative + correction for quantum mechanical effects. It is a unitless quantity. """ def __init__(self, hdf5_dbase_root=None, *args, **kwargs): if hdf5_dbase_root is None: @@ -57,6 +62,11 @@ def _itohintnonrel(self): data_path = '/'.join(['continuum', 'itohintnonrel']) return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path) + @property + def _klgfb(self): + data_path = '/'.join(['continuum', 'klgfb']) + return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path) + @u.quantity_input def free_free(self, temperature: u.K, atomic_number, charge_state, wavelength: u.angstrom) -> u.dimensionless_unscaled: r""" @@ -77,6 +87,17 @@ def free_free(self, temperature: u.K, atomic_number, charge_state, wavelength: u g_{ff} = \sum_{i,j=0}^{10} a_{ij}t^iU^j, where :math:`t=(\log{T} - 7.25)/1.25` and :math:`U=(\log{u} + 1.5)/2.5`. + + Parameters + ---------- + temperature : `~astropy.units.Quantity` + The temperature(s) for which to calculate the Gaunt factor + atomic_number : `int` + The atomic number of the emitting element + charge_state : `int`, + The charge state of the emitting ion + wavelength : `~astropy.units.Quantity` + The wavelength(s) at which to calculate the Gaunt factor """ gf_itoh = self._free_free_itoh(temperature, atomic_number, wavelength) gf_sutherland = self._free_free_sutherland(temperature, charge_state, wavelength) @@ -138,6 +159,13 @@ def free_free_total(self, temperature: u.K, charge_state) -> u.dimensionless_uns """ The total (wavelength-averaged) free-free Gaunt factor, used for calculating the total radiative losses from free-free emission. + + Parameters + ---------- + temperature : `~astropy.units.Quantity` + The temperature(s) for which to calculate the Gaunt factor + charge_state : `int`, + The charge state of the ion """ if charge_state == 0: return u.Quantity(np.zeros(temperature.shape)) @@ -155,6 +183,13 @@ def _free_free_itoh_integrated_nonrelativistic(self, temperature: u.K, charge_st """ The total (wavelength-averaged) non-relativistic 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 """ Ry = const.h * const.c * const.Ryd gamma_squared = (charge_state**2) * Ry / (const.k_B * temperature) @@ -173,6 +208,13 @@ def _free_free_itoh_integrated_relativistic(self, temperature: u.K, charge_state """ The total (wavelength-averaged) relativistic 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 """ z = (charge_state - 14.5) / 13.5 t = (np.log10(temperature)-7.25)/1.25 @@ -185,6 +227,34 @@ def _free_free_itoh_integrated_relativistic(self, temperature: u.K, charge_state summation[j] += self._itohintrel['a_ik'][i][k] * z**i * t**k return summation + @needs_dataset('klgfb') + @u.quantity_input + def free_bound(self, E_scaled: u.dimensionless_unscaled, n, l) -> u.dimensionless_unscaled: + r""" + Free-bound Gaunt factor as a function of scaled energy. + + The empirical fits are taken from Table 1 of :cite:t:`karzas_electron_1961`. + + Parameters + ---------- + E_scaled : `~astropy.units.Quantity` + A scaled energy, the ratio of photon energy divided by ionization energy. + n : `int` + The principal quantum number + l : `int` + The azimuthal quantum number + """ + index_nl = np.where(np.logical_and(self._klgfb['n'] == n, self._klgfb['l'] == l))[0] + # If there is no Gaunt factor for n, l, set it to 1 + if index_nl.shape == (0,): + gf = 1 + else: + gf_interp = splrep(self._klgfb['log_pe'][index_nl, :].squeeze(), + self._klgfb['log_gaunt_factor'][index_nl, :].squeeze()) + gf = np.exp(splev(E_scaled, gf_interp)) + + return gf + @u.quantity_input def free_bound_total(self, temperature: u.K, atomic_number, charge_state, n_0, @@ -197,13 +267,13 @@ def free_bound_total(self, temperature: u.K, atomic_number, charge_state, n_0, Parameters ---------- - temperature : + temperature : `~astropy.units.Quantity` The temperature(s) for which to calculate the Gaunt factor atomic_number : `int` The atomic number of the element - charge_state : `int`, + charge_state : `int` The charge state of the ion - n_0 : + n_0 : `int` The principal quantum number n of the ground state of the recombined ion ionization_potential : The ionization potential of the recombined ion diff --git a/fiasco/ions.py b/fiasco/ions.py index 0e6bb639..a7e1fb88 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -6,7 +6,7 @@ import numpy as np from functools import cached_property -from scipy.interpolate import CubicSpline, interp1d, PchipInterpolator, splev, splrep +from scipy.interpolate import CubicSpline, interp1d, PchipInterpolator from fiasco import proton_electron_ratio from fiasco.base import ContinuumBase, IonBase @@ -1507,7 +1507,6 @@ def _verner_cross_section(self, energy: u.erg) -> u.cm**2: return np.where(energy < self._verner['E_thresh'], 0., F.decompose().value) * self._verner['sigma_0'] - @needs_dataset('klgfb') @u.quantity_input def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg, n, l) -> u.cm**2: """ @@ -1525,15 +1524,8 @@ def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg, Orbital angular momentum number """ prefactor = (2**4)*const.h*(const.e.gauss**2)/(3*np.sqrt(3)*const.m_e*const.c) - index_nl = np.where(np.logical_and(self._klgfb['n'] == n, self._klgfb['l'] == l))[0] - # If there is no Gaunt factor for n, l, set it to 1 - if index_nl.shape == (0,): - gaunt_factor = 1 - else: - E_scaled = np.log(photon_energy/ionization_energy) - gf_interp = splrep(self._klgfb['log_pe'][index_nl, :].squeeze(), - self._klgfb['log_gaunt_factor'][index_nl, :].squeeze()) - gaunt_factor = np.exp(splev(E_scaled, gf_interp)) + E_scaled = np.log(photon_energy/ionization_energy) + gaunt_factor = self.gaunt_factor.free_bound(E_scaled, n, l) cross_section = prefactor * ionization_energy**2 * photon_energy**(-3) * gaunt_factor / n cross_section[np.where(photon_energy < ionization_energy)] = 0.*cross_section.unit return cross_section