Source code for rabacus.rad_src.plane

""" An plane parallel radiation source class. """ 


import numpy as np
from rabacus.utils import utils
from rabacus.constants import physical

from hm12 import HM12Spectrum
from thermal import ThermalSpectrum
from powerlaw import PowerlawSpectrum
from monochromatic import MonochromaticSpectrum
from source import Source

import species


__all__ = ['PlaneSource']



[docs]class PlaneSource(Source): r""" A plane parallel radiation source class. Args: `q_min` (float): minimum photon energy / Rydbergs [dimensionless] `q_max` (float): maximum photon energy / Rydbergs [dimensionless] `spectrum_type` (str): the spectral shape {``monochromatic``, ``hm12``, ``thermal``, ``powerlaw``, ``user``} Kwargs: `Nnu` (int): number of spectral samples, log spaced in energy `segmented` (bool): if ``True``, forces spectral samples at H and He ionization thresholds `px_fit_type` (str): source to use for photoionization cross section fits {``verner``} `verbose` (bool): verbose screen output? `z` (float): redshift, need if `spectrum_type` = ``hm12`` `T_eff` (float): effective temperature, need if `spectrum_type` = ``thermal`` `alpha` (float): powerlaw index, need if `spectrum_type` = ``powerlaw`` `user_E` (array): energy samples for user defined spectrum. should have units of energy. `user_shape` (array): shape of user defined spectrum. should be dimensionless Attributes: `U` (:class:`~rabacus.constants.units.Units`) `PC` (:class:`~rabacus.constants.physical.PhysicalConstants`) `H` (:class:`~rabacus.atomic.hydrogen.Hydrogen`) `He` (:class:`~rabacus.atomic.helium.Helium`) `PX` (:class:`~rabacus.atomic.photo_xsection.PhotoXsections`) `E` (array): Spectrum energy samples `lam` (array): Spectrum wavelength samples `nu` (array): Spectrum frequency samples `q` (array): Spectrum energy samples / Rydbergs `dFu_over_dnu` (array): Energy flux density |dFu_dnu| `monochromatic` (bool): ``True`` for monochromatic spectra otherwise ``False`` `grey` (object): Grey photoionization cross sections and energies. Only created if `segmented` = ``True``. `log` (object): Convenience object storing the log of other attributes `sigma` (object): photoionization cross sections and ratios `th` (object): all quantities evaluated at ionization thresholds `thin` (object): all optically thin quantities Notes: The fundamental characterization of the radiation field for this class is the energy flux density |dFu_dnu| with units of ``erg/(s cm^2 Hz)`` or energy flux :math:`F_{u}` with units of ``erg/(s cm^2)`` for monochromatic spectra. Integrating over frequency yields several other useful characterizations of the radiation field which we summarize below. .. math:: F_{u} = \int \frac{dF_{u}}{d\nu} d\nu .. math:: u = \frac{1}{c} \int \frac{dF_{u}}{d\nu} d\nu \rightarrow \frac{du}{d\nu} = \frac{1}{c} \frac{dF_{u}}{d\nu} .. math:: F_{n} = \int \frac{1}{h \nu} \frac{dF_{u}}{d\nu} d\nu \rightarrow \frac{dF_{n}}{d\nu} = \frac{1}{h \nu} \frac{dF_{u}}{d\nu} .. math:: n = \frac{1}{c} \int \frac{1}{h \nu} \frac{dF_{u}}{d\nu} d\nu \rightarrow \frac{dn}{d\nu} = \frac{1}{c} \frac{1}{h \nu} \frac{dF_{u}}{d\nu} +-----------+---------------------+----------------------+-----------+ | Variable | Description | Units | Symbol | +===========+=====================+======================+===========+ | dFu_dnu | energy flux density | [erg/(s cm^2 Hz)] | |dFu_dnu| | +-----------+---------------------+----------------------+-----------+ | thin.Fu | energy flux | [erg/(s cm^2)] | |Fu| | +-----------+---------------------+----------------------+-----------+ | thin.u | energy density | [erg/(cm^3)] | |u| | +-----------+---------------------+----------------------+-----------+ | thin.Fn | photon flux | [1/(s cm^2)] | |Fn| | +-----------+---------------------+----------------------+-----------+ | thin.n | photon density | [1/(cm^3)] | |n| | +-----------+---------------------+----------------------+-----------+ The photoionization and photoheating rates are also integrals over the energy flux density. .. math:: \Gamma_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{h \nu} \frac{dF_{u}}{d\nu} d\nu \rightarrow \frac{ d\Gamma_{\rm _{X}} }{d\nu} = \frac{ \sigma_{\rm _{X}} }{h \nu} \frac{dF_{u}}{d\nu} .. math:: H_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{h \nu} \frac{dF_{u}}{d\nu} h(\nu - \nu_{\rm _{X}}) d\nu \rightarrow \frac{ dH_{\rm _{X}} }{d\nu} = \frac{ \sigma_{\rm _{X}} }{h \nu} \frac{dF_{u}}{d\nu} h(\nu - \nu_{\rm _{X}}) where :math:`\scriptstyle{\rm X}` is one of :math:`\{ {\scriptstyle{\rm HI, \, HeI, \, HeII }} \}`. Note that the :math:`\sigma_{\rm _{X}}` are frequency dependent photoionization cross-sections and the :math:`\nu_{\rm _{X}}` are (scalar) frequencies that correspond to the ionization energies for a given species. In many applications, we prefer to use energy as a variable instead of frequency. The above equations can easily be recast with a change of variable :math:`E=h\nu`. .. math:: \Gamma_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} dE \rightarrow \frac{ d\Gamma_{\rm _{X}} }{dE} = \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} .. math:: H_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} (E - E_{\rm _{X}}) dE \rightarrow \frac{ dH_{\rm _{X}} }{dE} = \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} (E - E_{\rm _{X}}) where the :math:`E_{\rm _{X}}` are (scalar) ionization energies. +-----------+---------------------+----------------------+--------+ | Variable | Description | Units | Symbol | +===========+=====================+======================+========+ | thin.H1i | H1 photo ion. rate | [1/s] | |H1i| | +-----------+---------------------+----------------------+--------+ | thin.He1i | He1 photo ion. rate | [1/s] | |He1i| | +-----------+---------------------+----------------------+--------+ | thin.He2i | He2 photo ion. rate | [1/s] | |He2i| | +-----------+---------------------+----------------------+--------+ | thin.H1h | H1 photo heat rate | [erg/s] | |H1h| | +-----------+---------------------+----------------------+--------+ | thin.He1h | He1 photo heat rate | [erg/s] | |He1h| | +-----------+---------------------+----------------------+--------+ | thin.He2h | He2 photo heat rate | [erg/s] | |He2h| | +-----------+---------------------+----------------------+--------+ In addition, the class provides methods for calculating the shielded photoionization/heating rates given an optical depth. .. math:: \Gamma_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} e^{-\tau} dE \rightarrow \frac{ d\Gamma_{\rm _{X}} }{dE} = \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE}e^{-\tau} .. math:: H_{\rm _{X}} = \int \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} (E - E_{\rm _{X}}) e^{-\tau} dE \rightarrow \frac{ dH_{\rm _{X}} }{dE} = \frac{ \sigma_{\rm _{X}} }{E} \frac{dF_{u}}{dE} (E - E_{\rm _{X}}) e^{-\tau} where :math:`\tau = \sum_{\rm _X} \tau_{\rm _X}`, :math:`\tau_{\rm _X} = N_{\rm _X} \sigma_{\rm _X}`, :math:`N_{\rm _X}` are column densities, and :math:`\sigma_{\rm _X}` are energy dependent photoionization cross sections. .. note:: The :math:`\tau_{\rm _X}` are energy dependent optical depths while the parameters to the shielding functions (:func:`shld_H1i`, :func:`shld_H1h`, :func:`shld_He1i`, :func:`shld_He1h`, :func:`shld_He2i`, :func:`shld_He2h`) are the :math:`\tau_{\rm _X}` evaluated at the species threhold energy (i.e. scalars). +-----------+---------------------+----------------------+--------+ | Function | Description | Units | Symbol | +===========+=====================+======================+========+ | shld_H1i | H1 photo ion. rate | [1/s] | |H1i| | +-----------+---------------------+----------------------+--------+ | shld_He1i | He1 photo ion. rate | [1/s] | |He1i| | +-----------+---------------------+----------------------+--------+ | shld_He2i | He2 photo ion. rate | [1/s] | |He2i| | +-----------+---------------------+----------------------+--------+ | shld_H1h | H1 photo heat rate | [erg/s] | |H1h| | +-----------+---------------------+----------------------+--------+ | shld_He1h | He1 photo heat rate | [erg/s] | |He1h| | +-----------+---------------------+----------------------+--------+ | shld_He2h | He2 photo heat rate | [erg/s] | |He2h| | +-----------+---------------------+----------------------+--------+ .. seealso:: :class:`~rabacus.rad_src.background.BackgroundSource`, :class:`~rabacus.rad_src.point.PointSource` """ def __init__( self, q_min, q_max, spectrum_type, Nnu = 128, segmented = True, # puts spectrum points at E_th px_fit_type = 'verner', verbose = False, z = None, # only for spectrum_type = 'hm12' T_eff = None, # only for spectrum_type = 'thermal' alpha = None, # only for spectrum_type = 'powerlaw' user_E = None, # photon E array for spectrum_type='user' user_shape = None, # shape for spectrum_type='user' ): # attach specific input background #----------------------------------------------------------------- self.source_type = 'plane' # do generic initialization #----------------------------------------------------------------- self.initialize( q_min, q_max, spectrum_type, Nnu, segmented, px_fit_type, verbose, z, T_eff, alpha, user_E, user_shape ) # set spectral shape (from yval) #----------------------------------------------------------------- if spectrum_type == 'thermal': spec = ThermalSpectrum( self, self.T_eff ) elif spectrum_type == 'powerlaw': spec = PowerlawSpectrum( self, self.alpha ) elif spectrum_type == 'hm12': spec = HM12Spectrum( self, self.z ) elif spectrum_type == 'monochromatic': spec = MonochromaticSpectrum( self ) # set primary spectrum shape from spec.yvals #----------------------------------------------------------------- num = self.U.erg if self.monochromatic: den = self.U.s * self.U.cm**2 units = num / den self.Fu = np.ones( self.Nnu ) * spec.yvals * units self.Fu.__doc__ = 'energy flux' self.log.Fu = np.log10( self.Fu.magnitude ) else: den = self.U.Hz * self.U.s * self.U.cm**2 units = num / den if spectrum_type == 'user': self.dFu_over_dnu = np.ones( self.Nnu ) * user_shape * units else: self.dFu_over_dnu = np.ones( self.Nnu ) * spec.yvals * units if spectrum_type == 'hm12': self.dFu_over_dnu *= 4 * np.pi self.dFu_over_dnu.__doc__ = 'energy flux per unit frequency' self.log.dFu_over_dnu = np.log10( self.dFu_over_dnu.magnitude ) # set integrals #----------------------------------------------------------------- self.set_integrands() self.thin = OpticallyThinQuantities( self ) # calculate grey X-sections and energies # (note segmented = False if monochromatic = True) #----------------------------------------------------------------- if self.segmented: self.grey = GreyQuantities( self )
[docs] def set_integrands( self ): """ Defines tabulated functions of photon energy/frequency which represent integrands for optically thin quantities. All attributes defined in this function are integrands of the fundamental quantity, energy flux density or dFu_over_dnu which has units [erg/(s cm^2 Hz)]. """ # dont need to set these for monochromatic spectra #-------------------------------------------------------- if self.monochromatic: return # first we define differentials wrt frequency. these are # just for convenience. #-------------------------------------------------------- self._du_over_dnu = self.dFu_over_dnu / self.PC.c self._du_over_dnu.units = 'erg/cm^3/Hz' self._dFn_over_dnu = self.dFu_over_dnu / ( self.PC.h * self.nu ) self._dFn_over_dnu.units = '1/s/cm^2/Hz' self._dn_over_dnu = self._dFn_over_dnu / self.PC.c self._dn_over_dnu.units = '1/cm^3/Hz' #-------------------------------------------------------- self._dH1i_over_dnu = self._dFn_over_dnu * self.sigma.H1 self._dH1i_over_dnu.units = '1/s/Hz' self._dH1h_over_dnu = self._dH1i_over_dnu * (self.E - self.th.E_H1) self._dH1h_over_dnu.units = 'erg/s/Hz' self._dHe1i_over_dnu = self._dFn_over_dnu * self.sigma.He1 self._dHe1i_over_dnu.units = '1/s/Hz' self._dHe1h_over_dnu = self._dHe1i_over_dnu * (self.E - self.th.E_He1) self._dHe1h_over_dnu.units = 'erg/s/Hz' self._dHe2i_over_dnu = self._dFn_over_dnu * self.sigma.He2 self._dHe2i_over_dnu.units = '1/s/Hz' self._dHe2h_over_dnu = self._dHe2i_over_dnu * (self.E - self.th.E_He2) self._dHe2h_over_dnu.units = 'erg/s/Hz' # now we define differentials wrt energy #-------------------------------------------------------- self._dFu_over_dE = self.dFu_over_dnu / self.PC.h self._dFu_over_dE.units = 'erg/s/cm^2/eV' self._du_over_dE = self._dFu_over_dE / self.PC.c self._du_over_dE.units = 'erg/cm^3/eV' self._dFn_over_dE = self.dFu_over_dnu / ( self.PC.h * self.E ) self._dFn_over_dE.units = '1/s/cm^2/eV' self._dn_over_dE = self._dFn_over_dE / self.PC.c self._dn_over_dE.units = '1/cm^3/eV' #-------------------------------------------------------- self._dH1i_over_dE = self._dFn_over_dE * self.sigma.H1 self._dH1i_over_dE.units = '1/s/eV' self._dH1h_over_dE = self._dH1i_over_dE * (self.E - self.th.E_H1) self._dH1h_over_dE.units = 'erg/s/eV' self._dHe1i_over_dE = self._dFn_over_dE * self.sigma.He1 self._dHe1i_over_dE.units = '1/s/eV' self._dHe1h_over_dE = self._dHe1i_over_dE * (self.E - self.th.E_He1) self._dHe1h_over_dE.units = 'erg/s/eV' self._dHe2i_over_dE = self._dFn_over_dE * self.sigma.He2 self._dHe2i_over_dE.units = '1/s/eV' self._dHe2h_over_dE = self._dHe2i_over_dE * (self.E - self.th.E_He2) self._dHe2h_over_dE.units = 'erg/s/eV' # functions that calculate shielded integrals #====================================================================
[docs] def return_attenuation( self, tauH1_th, tauHe1_th, tauHe2_th ): """ Calculate the attenuation array. First the energy dependence of tau for each absorbing species (H1, He1, He2) is calculated using the input scalar optical depths at the ionization thresholds. Next an energy dependent total tau is calculated and exponentiated to arrive at the attenuation array atten = np.exp(-tau). Args: `tauH1_th` (float): H1 optical depth at the H1 ionizing threshold `tauHe1_th` (float): He1 optical depth at the He1 ionizing threshold `tauHe2_th` (float): He2 optical depth at the He2 ionizing threshold Returns: `atten` (array): attenuation as a function of energy """ assert isinstance(tauH1_th, float) assert isinstance(tauHe1_th, float) assert isinstance(tauHe2_th, float) tauH1 = tauH1_th * self.sigma.H1_ratio tauHe1 = tauHe1_th * self.sigma.He1_ratio tauHe2 = tauHe2_th * self.sigma.He2_ratio tau = tauH1 + tauHe1 + tauHe2 atten = np.exp( -tau ) assert atten.size == self.E.size return atten
[docs] def shld_H1i(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HI photoionization rate after passing through a column with optical depth tau Args: `tauH1_th` (float): H1 optical depth at the H1 ionizing threshold `tauHe1_th` (float): He1 optical depth at the He1 ionizing threshold `tauHe2_th` (float): He2 optical depth at the He2 ionizing threshold Returns: `H1i` (float): attenuated H1 photoionization rate """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: H1i = self.thin.H1i * atten else: H1i = utils.trap( self.E, self._dH1i_over_dE * atten ) H1i.units = '1/s' return H1i
[docs] def shld_H1h(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HI photoheating rate after passing through a column with optical depth tau Args: `tauH1_th` (float): H1 optical depth at the H1 ionizing threshold `tauHe1_th` (float): He1 optical depth at the He1 ionizing threshold `tauHe2_th` (float): He2 optical depth at the He2 ionizing threshold Returns: `H1h` (float): attenuated H1 photoheating rate """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: H1h = self.thin.H1h * atten else: H1h = utils.trap( self.E, self._dH1h_over_dE * atten ) H1h.units = 'erg/s' return H1h
[docs] def shld_He1i(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HeI photoionization rate after passing through a column with optical depth tau. .. seealso:: :func:`shld_H1i` """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: He1i = self.thin.He1i * atten else: He1i = utils.trap( self.E, self._dHe1i_over_dE * atten ) He1i.units = '1/s' return He1i
[docs] def shld_He1h(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HeI photoheating rate after passing through a column with optical depth tau. .. seealso:: :func:`shld_H1h` """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: He1h = self.thin.He1h * atten else: He1h = utils.trap( self.E, self._dHe1h_over_dE * atten ) He1h.units = 'erg/s' return He1h
[docs] def shld_He2i(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HeII photoionization rate after passing through a column with optical depth tau. .. seealso:: :func:`shld_H1i` """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: He2i = self.thin.He2i * atten else: He2i = utils.trap( self.E, self._dHe2i_over_dE * atten ) He2i.units = '1/s' return He2i
[docs] def shld_He2h(self, tauH1_th, tauHe1_th, tauHe2_th): """ Integrates spectrum to calculate the HeII photoheating rate after passing through a column with optical depth tau. .. seealso:: :func:`shld_H1h` """ atten = self.return_attenuation( tauH1_th, tauHe1_th, tauHe2_th ) if self.monochromatic: He2h = self.thin.He2h * atten else: He2h = utils.trap( self.E, self._dHe2h_over_dE * atten ) He2h.units = 'erg/s' return He2h # functions that scale and normalize the spectrum #====================================================================
[docs] def scale_spectrum( self, fac ): """ Changes normalization of spectrum and recalculates spectral integrands. Args: `fac` (float): multiplicative factor """ if self.monochromatic: self.Fu *= fac else: self.dFu_over_dnu *= fac self.set_integrands() self.thin = OpticallyThinQuantities( self )
[docs] def normalize_n( self, n ): """ Normalize the spectrum such that the number density of photons between `q_min` and `q_max` is `n`. Args: `n` (float): target photon number density """ if not hasattr(n,'units'): raise utils.NeedUnitsError, '\n n must have units, e.g. cm^-3 \n' else: n.units = 'cm^-3' fac = n / self.thin.n self.scale_spectrum( fac )
[docs] def normalize_Fn( self, Fn ): """ Normalize the spectrum such that the photon flux between `q_min` and `q_max` is `Fn`. Args: `Fn` (float): target photon flux """ if not hasattr(Fn,'units'): msg = '\n Fn must have units, e.g. s^-1 cm^-2 \n' raise utils.NeedUnitsError, msg else: Fn.units = 's^-1 * cm^-2' fac = Fn / self.thin.Fn self.scale_spectrum( fac )
[docs] def normalize_H1i( self, H1i ): """ Normalize the spectrum such that the H1 photoionization rate integral yields `H1i`. Args: `H1i` (float): target photoionization rate """ if not hasattr(H1i,'units'): raise utils.NeedUnitsError, '\n H1i must have units, e.g. s^-1 \n' else: H1i.units = 's^-1' fac = H1i / self.thin.H1i self.scale_spectrum( fac ) # functions for "grey" quantities #====================================================================
class GreyQuantities: """ This class stores grey photoionization cross-sections and energies for each of the absorbing species (H1, He1, He2). Note that this requires a segmented spectrum. """ def __init__( self, rad_src ): # make sure we have what we need #-------------------------------------------- assert rad_src.segmented == True # setup containers #-------------------------------------------- self.sigma = species.AbsorbingSpecies() self.E = species.AbsorbingSpecies() # integrate between H1 and He1 thresholds #-------------------------------------------- ii = rad_src.th.i_H1 ff = rad_src.th.i_He1 xx = rad_src.E[ii:ff] yy = rad_src._dH1i_over_dE[ii:ff] H1i_thin = utils.trap( xx, yy ) yy = rad_src._dFn_over_dE[ii:ff] Fn_thin = utils.trap( xx, yy ) self.sigma.H1 = H1i_thin / Fn_thin self.sigma.H1.__doc__ = 'grey H1 photoionization cross section' log_sig = np.log10( rad_src.sigma.H1[ii:ff].magnitude ) log_gry = np.log10( self.sigma.H1.magnitude ) log_E = np.log10( rad_src.E[ii:ff].magnitude ) i_lo = np.argmin( np.abs( log_sig - log_gry ) ) if log_sig[i_lo] < log_gry: i_hi = i_lo + 1 else: i_hi = i_lo i_lo = i_lo - 1 dlogsig = log_sig[i_hi] - log_sig[i_lo] frac = (log_gry - log_sig[i_lo]) / dlogsig dlogE = log_E[i_hi] - log_E[i_lo] grey_E = 10.0**(log_E[i_lo] + frac * dlogE) self.E.H1 = grey_E * rad_src.U.eV self.E.H1.__doc__ = 'grey H1 energy' # integrate between He1 and He2 thresholds #-------------------------------------------- ii = rad_src.th.i_He1 ff = rad_src.th.i_He2 xx = rad_src.E[ii:ff] yy = rad_src._dHe1i_over_dE[ii:ff] He1i_thin = utils.trap( xx, yy ) yy = rad_src._dFn_over_dE[ii:ff] Fn_thin = utils.trap( xx, yy ) self.sigma.He1 = He1i_thin / Fn_thin self.sigma.He1.__doc__ = 'grey He1 photoionization cross section' log_sig = np.log10( rad_src.sigma.He1[ii:ff].magnitude ) log_gry = np.log10( self.sigma.He1.magnitude ) log_E = np.log10( rad_src.E[ii:ff].magnitude ) i_lo = np.argmin( np.abs( log_sig - log_gry ) ) if log_sig[i_lo] < log_gry: i_hi = i_lo + 1 else: i_hi = i_lo i_lo = i_lo - 1 dlogsig = log_sig[i_hi] - log_sig[i_lo] frac = (log_gry - log_sig[i_lo]) / dlogsig dlogE = log_E[i_hi] - log_E[i_lo] grey_E = 10.0**(log_E[i_lo] + frac * dlogE) self.E.He1 = grey_E * rad_src.U.eV self.E.He1.__doc__ = 'grey He1 energy' # integrate above He2 threshold #-------------------------------------------- ii = rad_src.th.i_He2 xx = rad_src.E[ii:] yy = rad_src._dHe2i_over_dE[ii:] He2i_thin = utils.trap( xx, yy ) yy = rad_src._dFn_over_dE[ii:] Fn_thin = utils.trap( xx, yy ) self.sigma.He2 = He2i_thin / Fn_thin self.sigma.He2.__doc__ = 'grey He2 photoionization cross section' log_sig = np.log10( rad_src.sigma.He2[ii:].magnitude ) log_gry = np.log10( self.sigma.He2.magnitude ) log_E = np.log10( rad_src.E[ii:].magnitude ) i_lo = np.argmin( np.abs( log_sig - log_gry ) ) if log_sig[i_lo] < log_gry: i_hi = i_lo + 1 else: i_hi = i_lo i_lo = i_lo - 1 dlogsig = log_sig[i_hi] - log_sig[i_lo] frac = (log_gry - log_sig[i_lo]) / dlogsig dlogE = log_E[i_hi] - log_E[i_lo] grey_E = 10.0**(log_E[i_lo] + frac * dlogE) self.E.He2 = grey_E * rad_src.U.eV self.E.He2.__doc__ = 'grey He2 energy' class OpticallyThinQuantities: """ This class stores optically thin characterizations of the radiation field. """ def __init__ ( self, rad_src ): self.PC = physical.PhysicalConstants() # dont need to integrate for monochromatic spectra #----------------------------------------------------------------- if rad_src.monochromatic: self.Fu = rad_src.Fu self.u = self.Fu / self.PC.c self.Fn = self.Fu / rad_src.E self.n = self.Fn / self.PC.c self.H1i = self.Fn * rad_src.sigma.H1 self.H1h = self.H1i * (rad_src.E - rad_src.th.E_H1) self.He1i = self.Fn * rad_src.sigma.He1 self.He1h = self.He1i * (rad_src.E - rad_src.th.E_He1) self.He2i = self.Fn * rad_src.sigma.He2 self.He2h = self.He2i * (rad_src.E - rad_src.th.E_He2) # integrate for polychromatic spectra #----------------------------------------------------------------- else: self.Fu = utils.trap( rad_src.E, rad_src._dFu_over_dE ) self.u = utils.trap( rad_src.E, rad_src._du_over_dE ) self.Fn = utils.trap( rad_src.E, rad_src._dFn_over_dE ) self.n = utils.trap( rad_src.E, rad_src._dn_over_dE ) self.H1i = utils.trap( rad_src.E, rad_src._dH1i_over_dE ) self.H1h = utils.trap( rad_src.E, rad_src._dH1h_over_dE ) self.He1i = utils.trap(rad_src.E, rad_src._dHe1i_over_dE) self.He1h = utils.trap(rad_src.E, rad_src._dHe1h_over_dE) self.He2i = utils.trap(rad_src.E, rad_src._dHe2i_over_dE) self.He2h = utils.trap(rad_src.E, rad_src._dHe2h_over_dE) # set units and docstrings #----------------------------------------------------------------- self.Fu.units = 'erg/s/cm^2' self.Fu.__doc__ = 'optically thin energy flux' self.u.units = 'erg/cm^3' self.u.__doc__ = 'optically thin energy density' self.Fn.units = '1/s/cm^2' self.Fn.__doc__ = 'optically thin photon number flux' self.n.units = '1/cm^3' self.n.__doc__ = 'optically thin photon number density' self.H1i.units = '1/s' self.H1i.__doc__ = 'optically thin H1 photoionization rate' self.H1h.units = 'erg/s' self.H1h.__doc__ = 'optically thin H1 photoheating rate' self.He1i.units = '1/s' self.He1i.__doc__ = 'optically thin He1 photoionization rate' self.He1h.units = 'erg/s' self.He1h.__doc__ = 'optically thin He1 photoheating rate' self.He2i.units = '1/s' self.He2i.__doc__ = 'optically thin He2 photoionization rate' self.He2h.units = 'erg/s' self.He2h.__doc__ = 'optically thin He2 photoheating rate'