Source code for rabacus.solvers.calc_slab

""" Solves a slab with plane parallel radiation incident from one side. """

import numpy as np

import one_zone as ozn
from rabacus.atomic import chemistry
from rabacus.atomic import cooling
from rabacus.constants import units
from rabacus.utils import utils



__all__ = ['Slab']


[docs]class Slab: """ Solves a slab geometry with plane parallel radiation incident from one side. The slab is divided into `Nl` layers. Args: `Edges` (array): Distance to layer edges `T` (array): Temperature in each layer `nH` (array): hydrogen number density in each layer `nHe` (array): helium number density in each layer `rad_src` (:class:`~rabacus.rad_src.plane.PlaneSource`): Plane source .. note:: The arrays `T`, `nH`, and `nHe` must all be the same size. This size determines the number of shells, `Nl`. `Edges` determines the positions of the layer edges and must have `Nl+1` entries. Kwargs: `rec_meth` (string): How to treat recombinations {``fixed``, ``thresh``} `fixed_fcA` (float): If `rec_meth` = ``fixed``, constant caseA fraction `thresh_P_A` (float): If `rec_meth` = ``thresh``, this is the probability of absorption, P = 1 - exp(-tau), at which the transition to caseB rates begins. `thresh_P_B` (float): If `rec_meth` = ``thresh``, this is the probability of absorption, P = 1 - exp(-tau), at which the transition to caseB ends. `thresh_xmfp` (float): Determines the distance to probe when calculating optical depth for the caseA to caseB transition. Specifies a multiple of the mean free path for fully neutral gas. In equation form, L = thresh_xmfp / (n * sigma) `atomic_fit_name` (string): Source for atomic rate fits {``hg97``} `find_Teq` (bool): If ``False``, use fixed input T, if ``True`` solve for equilibrium T `z` (float): Redshift, only need if `find_Teq` = ``True`` `verbose` (bool): Verbose output? `tol` (float): tolerance for all convergence tests `thin` (bool): if ``True`` only solves optically thin Attributes: `U` (:class:`~rabacus.constants.units.Units`) `z_c` (array): distance from surface of slab to center of layer `dz` (array): thickness of layer `NH_c` (array): H column density from surface of slab to center of layer `dNH` (array): H column density through layer `NH1_c` (array): HI column density from surface of slab to center of layer `dNH1` (array): HI column density through layer `H1i` (array): HI photo-ionization rate `H1h` (array): HI photo-heating rate `xH1` (array): H neutral fraction nHI / nH `xH2` (array): H ionized fraction nHII / nH `ne` (array): electron number density `fcA_H2` (array): HII case A fraction `cool` (array): cooling rate [erg / (cm^3 K)] `heat` (array): heating rate [erg / (cm^3 K)] `heatH1` (array): contribution to `heat` from H1 photoheating `dtauH1_th` (array): HI optical depth at H1 ionizing threshold through layer `tauH1_th_lo` (array): HI optical depth below this layer `tauH1_th_hi` (array): HI optical depth above this layer `NH1_thru` (float): HI column density through entire slab `Nl` (int): Number of layers `itr` (int): Number of iterations to converge .. note:: For many of the attributes above, there are analagous versions for helium. We also note that the source of the atomic rates fit is stored in the variable `atomic_fit_name`, but the source of the photoionization cross section fits are stored in the point source object in the variable `px_fit_type`. Examples: The following code will solve a monochromatic source incident onto a slab with fixed recombination rates and temperatues. In this case we use a Haardt and Madau 2012 spectrum to calculate a grey photon energy which is then used to create a monochromatic plane source :: import numpy as np import rabacus as ra # create Haardt and Madau 2012 source #--------------------------------------------------------------- q_min = 1.0; q_max = 4.0e2; z = 3.0 hm12 = ra.rad_src.PlaneSource( q_min, q_max, 'hm12', z=z ) # get grey photon energy and create a monochromatic plane source #--------------------------------------------------------------- q_mono = hm12.grey.E.He2 / hm12.PX.Eth_H1 mono = ra.rad_src.PlaneSource( q_mono, q_mono, 'monochromatic' ) mono.normalize_H1i( hm12.thin.H1i ) # describe slab #--------------------------------------------------------------- Nl = 512; Yp = 0.24 nH = np.ones(Nl) * 1.5e-2 / ra.U.cm**3 nHe = nH * Yp * 0.25 / (1.0-Yp) Tslab = np.ones(Nl) * 1.0e4 * ra.U.K Lslab = 200.0 * ra.U.kpc Edges = np.linspace( 0.0 * ra.U.kpc, Lslab, Nl+1 ) # solve slab #--------------------------------------------------------------- slab = ra.solvers.Slab( Edges, Tslab, nH, nHe, mono, rec_meth='fixed', fixed_fcA=1.0 ) """ def __init__(self, Edges, T, nH, nHe, rad_src, # rec_meth = "fixed", fixed_fcA = 1.0, thresh_P_A = 0.01, thresh_P_B = 0.99, thresh_xmfp = 3.0e1, atomic_fit_name = "hg97", find_Teq = False, z = None, verbose = False, tol = 1.0e-10, thin = False, ): # attach units #----------------------------------------------- self.U = units.Units() # check input #----------------------------------------------- assert( Edges.size == T.size + 1 ) assert( T.shape == nH.shape == nHe.shape ) assert( T.shape == (T.size,) ) assert( rec_meth == 'fixed' or rec_meth == 'thresh' ) if find_Teq and z==None: msg = 'need to supply redshift if finding equilibrium temperature' raise utils.InputError(msg) if rad_src.source_type != 'plane': msg = 'source type needs to be plane' raise utils.InputError(msg) # set units #----------------------------------------------- Edges.units = 'cm' T.units = 'K' nH.units = '1.0/cm^3' nHe.units = '1.0/cm^3' # attach input #----------------------------------------------- self.Edges = Edges.copy() self.T = T.copy() self.nH = nH.copy() self.nHe = nHe.copy() self.rad_src = rad_src self.rec_meth = rec_meth if self.rec_meth == 'fixed': self.fixed_fcA = fixed_fcA elif self.rec_meth == 'thresh': assert thresh_xmfp > 0.0 assert thresh_P_B > thresh_P_A self.thresh_xmfp = thresh_xmfp self.thresh_P_A = thresh_P_A self.thresh_P_B = thresh_P_B self.thresh_dPAB = thresh_P_B - thresh_P_A self.thresh_tau_A = -np.log( 1.0 - thresh_P_A ) self.thresh_tau_B = -np.log( 1.0 - thresh_P_B ) self.atomic_fit_name = atomic_fit_name self.find_Teq = find_Teq self.verbose = verbose self.tol = tol self.thin = thin if find_Teq: self.z = z # initialize slab #----------------------------------------------- self.init_slab() self.set_optically_thin() if self.thin: return # solve slab (sweep until convergence) #----------------------------------------------------------- conv_old = np.sum( self.ne ) not_converged = True self.itr = 0 while not_converged: self.sweep_slab() conv_new = np.sum( self.ne ) if np.abs( conv_new/conv_old - 1.0 ) < self.tol: not_converged = False conv_old = conv_new self.itr += 1 # finalize slab #----------------------------------------------------------- self.finalize_slab()
[docs] def init_slab( self ): """ Initialize slab values. """ if self.verbose: print 'begin initialization' # instantiate atomic rates (optically thin) #----------------------------------------------- if self.rec_meth == 'fixed': fcA_H2 = self.fixed_fcA fcA_He2 = self.fixed_fcA fcA_He3 = self.fixed_fcA elif self.rec_meth == 'thresh': fcA_H2 = 1.0 fcA_He2 = 1.0 fcA_He3 = 1.0 kchem = chemistry.ChemistryRates( self.T[0], fcA_H2, fcA_He2, fcA_He3, H1i = self.rad_src.thin.H1i, He1i = self.rad_src.thin.He1i, He2i = self.rad_src.thin.He2i, fit_name = self.atomic_fit_name ) kcool = cooling.CoolingRates( self.T[0], fcA_H2, fcA_He2, fcA_He3, H1h = self.rad_src.thin.H1h, He1h = self.rad_src.thin.He1h, He2h = self.rad_src.thin.He2h, fit_name = self.atomic_fit_name ) self.kchem = kchem self.kcool = kcool if self.verbose: print ' created kchem and kcool' # setup arrays #----------------------------------------------- self.Nl = self.nH.size Nl = self.Nl self.dz = self.Edges[1:] - self.Edges[0:-1] self.z_c = self.Edges[0:-1] + 0.5 * self.dz self.xH1 = np.zeros(Nl) self.dNH1 = np.zeros(Nl) / self.U.cm**2 self.dtauH1_th = np.zeros(Nl) self.tauH1_th_lo = np.zeros(Nl) # H1 optical depth below layer self.tauH1_th_hi = np.zeros(Nl) # H1 optical depth above layer self.H1i = np.zeros(Nl) / self.U.s self.H1h = np.zeros(Nl) * self.U.eV / self.U.s self.fcA_H2 = np.zeros(Nl) self.xHe1 = np.zeros(Nl) self.dNHe1 = np.zeros(Nl) / self.U.cm**2 self.dtauHe1_th = np.zeros(Nl) self.tauHe1_th_lo = np.zeros(Nl) # He1 optical depth below layer self.tauHe1_th_hi = np.zeros(Nl) # He1 optical depth above layer self.He1i = np.zeros(Nl) / self.U.s self.He1h = np.zeros(Nl) * self.U.eV / self.U.s self.fcA_He2 = np.zeros(Nl) self.xHe2 = np.zeros(Nl) self.dNHe2 = np.zeros(Nl) / self.U.cm**2 self.dtauHe2_th = np.zeros(Nl) self.tauHe2_th_lo = np.zeros(Nl) # He2 optical depth below layer self.tauHe2_th_hi = np.zeros(Nl) # He2 optical depth above layer self.He2i = np.zeros(Nl) / self.U.s self.He2h = np.zeros(Nl) * self.U.eV / self.U.s self.fcA_He3 = np.zeros(Nl) self.cool = np.zeros(Nl) * self.U.eV / (self.U.s * self.U.cm**3) self.heat = np.zeros(Nl) * self.U.eV / (self.U.s * self.U.cm**3) self.heatH1 = np.zeros(Nl) * self.U.eV / (self.U.s * self.U.cm**3) self.heatHe1 = np.zeros(Nl) * self.U.eV / (self.U.s * self.U.cm**3) self.heatHe2 = np.zeros(Nl) * self.U.eV / (self.U.s * self.U.cm**3) self.xH2 = np.zeros(Nl) self.xHe3 = np.zeros(Nl) self.ne = np.zeros(Nl) / self.U.cm**3 if self.verbose: print ' created simple arrays' # calc NH and NHe between S and center of each layer #---------------------------------------- self.dNH = self.dz * self.nH self.dNHe = self.dz * self.nHe self.NH_c = np.zeros(Nl) / self.U.cm**2 self.NHe_c = np.zeros(Nl) / self.U.cm**2 for i in xrange(Nl): self.NH_c[i] = np.sum( self.dNH[0:i] ) + 0.5 * self.dNH[i] self.NHe_c[i] = np.sum( self.dNHe[0:i] ) + 0.5 * self.dNHe[i] if self.verbose: print ' created NH_c and NHe_c' # calc mean free path when fully neutral for each shell # also calculate the bounds for each shell #-------------------------------------------------------- if self.rec_meth == 'thresh': L = self.Edges[-1] self.mfp_H1 = 1.0 / ( self.nH * self.rad_src.th.sigma_H1 ) self.L_H1 = self.mfp_H1 * self.thresh_xmfp if np.any( self.L_H1 > L ): indx = np.where( self.L_H1 > L ) self.L_H1[indx] = L self.mfp_He1 = 1.0 / ( self.nHe * self.rad_src.th.sigma_He1 ) self.L_He1 = self.mfp_He1 * self.thresh_xmfp if np.any( self.L_He1 > L ): indx = np.where( self.L_He1 > L ) self.L_He1[indx] = L self.mfp_He2 = 1.0 / ( self.nHe * self.rad_src.th.sigma_He2 ) self.L_He2 = self.mfp_He2 * self.thresh_xmfp if np.any( self.L_He2 > L ): indx = np.where( self.L_He2 > L ) self.L_He2[indx] = L self.ii_H1 = np.zeros(Nl, dtype=int) self.ff_H1 = np.zeros(Nl, dtype=int) self.ii_He1 = np.zeros(Nl, dtype=int) self.ff_He1 = np.zeros(Nl, dtype=int) self.ii_He2 = np.zeros(Nl, dtype=int) self.ff_He2 = np.zeros(Nl, dtype=int) for i in xrange(Nl): self.ii_H1[i], self.ff_H1[i] = \ self.return_bounds( i, self.L_H1[i] ) self.ii_He1[i], self.ff_He1[i] = \ self.return_bounds( i, self.L_He1[i] ) self.ii_He2[i], self.ff_He2[i] = \ self.return_bounds( i, self.L_He2[i] ) if self.verbose: print ' created index arrays for fcA' if self.verbose: print 'initialization complete'
[docs] def return_bounds( self, i, Ltarget ): """ Given a shell number and a distance, return the indices such that self.r[i] - self.r[ii] > Ltarget and self.r[ff] - self.r[i] > Ltarget. """ L = 0.0 * self.U.cm k = 1 while L < Ltarget: ii = i-k if ii < 0: ii = 0 L = 1.5 * Ltarget else: L = L + self.dz[ii] k += 1 L = 0.0 * self.U.cm k = 1 while L < Ltarget: ff = i+k if ff >= self.Nl: ff = self.Nl-1 L = 1.5 * Ltarget else: L = L + self.dz[ff] k += 1 return ii,ff
[docs] def set_optically_thin( self ): """ Set optically thin ionzation state """ # initialize a chemistry object #---------------------------------------- if self.rec_meth == 'fixed': self.fcA_H2 = np.ones(self.Nl) * self.fixed_fcA self.fcA_He2 = np.ones(self.Nl) * self.fixed_fcA self.fcA_He3 = np.ones(self.Nl) * self.fixed_fcA else: self.fcA_H2 = np.ones(self.Nl) self.fcA_He2 = np.ones(self.Nl) self.fcA_He3 = np.ones(self.Nl) kchem = chemistry.ChemistryRates( self.T, self.fcA_H2, self.fcA_He2, self.fcA_He3, H1i = np.ones(self.Nl) * self.rad_src.thin.H1i, He1i = np.ones(self.Nl) * self.rad_src.thin.He1i, He2i = np.ones(self.Nl) * self.rad_src.thin.He2i, fit_name = self.atomic_fit_name, ) kcool = cooling.CoolingRates( self.T, self.fcA_H2, self.fcA_He2, self.fcA_He3, H1h = np.ones(self.Nl) * self.rad_src.thin.H1h, He1h = np.ones(self.Nl) * self.rad_src.thin.He1h, He2h = np.ones(self.Nl) * self.rad_src.thin.He2h, fit_name = self.atomic_fit_name, ) if self.find_Teq: x = ozn.Solve_PCTE( self.nH, self.nHe, kchem, kcool, self.z, self.tol ) self.T = x.Teq else: x = ozn.Solve_PCE( self.nH, self.nHe, kchem, self.tol ) self.xH1 = x.H1 self.xH2 = x.H2 self.xHe1 = x.He1 self.xHe2 = x.He2 self.xHe3 = x.He3 # summarize #---------------------------------------- self.dNH1 = self.dNH * self.xH1 self.dNHe1 = self.dNHe * self.xHe1 self.dNHe2 = self.dNHe * self.xHe2 self.dtauH1_th = self.dNH1 * self.rad_src.th.sigma_H1 self.dtauHe1_th = self.dNHe1 * self.rad_src.th.sigma_He1 self.dtauHe2_th = self.dNHe2 * self.rad_src.th.sigma_He2 self.ne = self.xH2 * self.nH + \ ( self.xHe2 + 2.0 * self.xHe3 ) * self.nHe
[docs] def set_taus( self, i ): """ Calculate optical depth at the ionization thresholds of each species above and below this layer. Does not include a contribution from the layer itself. Args: `i` (int): layer index """ self.tauH1_th_lo[i] = np.sum( self.dtauH1_th[0:i] ) self.tauH1_th_hi[i] = np.sum( self.dtauH1_th[i+1:self.Nl] ) self.tauHe1_th_lo[i] = np.sum( self.dtauHe1_th[0:i] ) self.tauHe1_th_hi[i] = np.sum( self.dtauHe1_th[i+1:self.Nl] ) self.tauHe2_th_lo[i] = np.sum( self.dtauHe2_th[0:i] ) self.tauHe2_th_hi[i] = np.sum( self.dtauHe2_th[i+1:self.Nl] )
[docs] def return_fcA( self, tau_th ): """ Given the optical depth at nu_th for any given species, returns the case A fraction. """ if tau_th < self.thresh_tau_A: fcA = 1.0 elif tau_th > self.thresh_tau_B: fcA = 0.0 else: P = 1.0 - np.exp( -tau_th ) delta = P - self.thresh_P_A fcB = delta / self.thresh_dPAB fcA = 1.0 - fcB return fcA
[docs] def set_fcA( self, i ): """ Calculate case A fraction for each species in this layer. Args: `i` (int): layer index """ if self.rec_meth == 'fixed': self.fcA_H2[i] = self.fixed_fcA self.fcA_He2[i] = self.fixed_fcA self.fcA_He3[i] = self.fixed_fcA elif self.rec_meth == 'thresh': ii = self.ii_H1[i]; ff = self.ff_H1[i] tau_lo = np.sum( self.dtauH1_th[ii:i] ) fcA_lo = self.return_fcA( tau_lo ) tau_hi = np.sum( self.dtauH1_th[i+1:ff+1] ) fcA_hi = self.return_fcA( tau_hi ) self.fcA_H2[i] = (fcA_lo + fcA_hi) * 0.5 ii = self.ii_He1[i]; ff = self.ff_He1[i] tau_lo = np.sum( self.dtauHe1_th[ii:i] ) fcA_lo = self.return_fcA( tau_lo ) tau_hi = np.sum( self.dtauHe1_th[i+1:ff+1] ) fcA_hi = self.return_fcA( tau_hi ) self.fcA_He2[i] = (fcA_lo + fcA_hi) * 0.5 ii = self.ii_He2[i]; ff = self.ff_He2[i] tau_lo = np.sum( self.dtauHe2_th[ii:i] ) fcA_lo = self.return_fcA( tau_lo ) tau_hi = np.sum( self.dtauHe2_th[i+1:ff+1] ) fcA_hi = self.return_fcA( tau_hi ) self.fcA_He3[i] = (fcA_lo + fcA_hi) * 0.5
[docs] def set_photoion_rates( self, i ): """ set photoionization rates for this layer Args: `i` (int): layer index """ tauH1_th = np.float( self.tauH1_th_lo[i] + 0.5 * self.dtauH1_th[i] ) tauHe1_th = np.float( self.tauHe1_th_lo[i] + 0.5 * self.dtauHe1_th[i] ) tauHe2_th = np.float( self.tauHe2_th_lo[i] + 0.5 * self.dtauHe2_th[i] ) self.H1i[i] = self.rad_src.shld_H1i( tauH1_th, tauHe1_th, tauHe2_th ) self.He1i[i] = self.rad_src.shld_He1i( tauH1_th, tauHe1_th, tauHe2_th ) self.He2i[i] = self.rad_src.shld_He2i( tauH1_th, tauHe1_th, tauHe2_th )
[docs] def set_photoheat_rates( self, i ): """ set photoheating rates for this layer Args: `i` (int): layer index """ tauH1_th = np.float( self.tauH1_th_lo[i] + 0.5 * self.dtauH1_th[i] ) tauHe1_th = np.float( self.tauHe1_th_lo[i] + 0.5 * self.dtauHe1_th[i] ) tauHe2_th = np.float( self.tauHe2_th_lo[i] + 0.5 * self.dtauHe2_th[i] ) self.H1h[i] = self.rad_src.shld_H1h( tauH1_th, tauHe1_th, tauHe2_th ) self.He1h[i] = self.rad_src.shld_He1h( tauH1_th, tauHe1_th, tauHe2_th ) self.He2h[i] = self.rad_src.shld_He2h( tauH1_th, tauHe1_th, tauHe2_th )
[docs] def sweep_slab(self): """ Performs one sweep through slab. """ for i in xrange(self.Nl): if self.verbose: print 'i = ', i # calc tauXX above and below this layer # (constant during iterations) #---------------------------------------- self.set_taus( i ) # calculate fcA (average over directions) # (constant during iterations) #---------------------------------------- self.set_fcA( i ) # iterate until we have convergence in the optical # depth through this layer #-------------------------------------------------- conv_old = self.dtauH1_th[i] + \ self.dtauHe1_th[i] + self.dtauHe2_th[i] not_converged = True itr = 0 while not_converged: # calculate photoion / chemistry rates #---------------------------------------- self.set_photoion_rates( i ) self.kchem.set( self.T[i], self.fcA_H2[i], self.fcA_He2[i], self.fcA_He3[i], H1i = self.H1i[i], He1i = self.He1i[i], He2i = self.He2i[i] ) # calculate photoheat / cooling rates #---------------------------------------- self.set_photoheat_rates( i ) self.kcool.set( self.T[i], self.fcA_H2[i], self.fcA_He2[i], self.fcA_He3[i], H1h = self.H1h[i], He1h = self.He1h[i], He2h = self.He2h[i] ) # if we are finding the equilibrium temperature # we need to call Solve_PCTE #---------------------------------------- if self.find_Teq: x = ozn.Solve_PCTE( np.ones(1) * self.nH[i], np.ones(1) * self.nHe[i], self.kchem, self.kcool, self.z, self.tol ) self.T[i] = x.Teq # otherwise we call Solve_PCE #---------------------------------------- else: x = ozn.Solve_PCE( np.ones(1) * self.nH[i], np.ones(1) * self.nHe[i], self.kchem, self.tol ) # set the ionization fractions in the object #---------------------------------------- self.xH1[i] = x.H1 self.xH2[i] = x.H2 self.xHe1[i] = x.He1 self.xHe2[i] = x.He2 self.xHe3[i] = x.He3 # calculate heating rates in layer #---------------------------------------- self.heatH1[i] = self.H1h[i] * self.nH[i] * self.xH1[i] self.heatHe1[i] = self.He1h[i] * self.nHe[i] * self.xHe1[i] self.heatHe2[i] = self.He2h[i] * self.nHe[i] * self.xHe2[i] self.heat[i] = self.heatH1[i] + \ self.heatHe1[i] + self.heatHe2[i] # calculate electron density in layer #---------------------------------------- self.ne[i] = self.xH2[i] * self.nH[i] + \ ( self.xHe2[i] + 2 * self.xHe3[i] ) * self.nHe[i] # calculate tauXX through layer #---------------------------------------- self.dNH1[i] = self.dNH[i] * self.xH1[i] self.dNHe1[i] = self.dNHe[i] * self.xHe1[i] self.dNHe2[i] = self.dNHe[i] * self.xHe2[i] self.dtauH1_th[i] = self.dNH1[i] * self.rad_src.th.sigma_H1 self.dtauHe1_th[i] = self.dNHe1[i] * self.rad_src.th.sigma_He1 self.dtauHe2_th[i] = self.dNHe2[i] * self.rad_src.th.sigma_He2 # check convergence #---------------------------------------- conv_new = self.dtauH1_th[i] + \ self.dtauHe1_th[i] + self.dtauHe2_th[i] if np.abs( conv_new/conv_old - 1.0 ) < self.tol: not_converged = False conv_old = conv_new itr += 1
[docs] def finalize_slab( self ): """ Calculate some final values using fully solved slab. """ # calculate column densities up to each layer #----------------------------------------------------------- self.NH1_c = np.zeros(self.Nl) / self.U.cm**2 self.NHe1_c = np.zeros(self.Nl) / self.U.cm**2 self.NHe2_c = np.zeros(self.Nl) / self.U.cm**2 for i in xrange(self.Nl): self.NH1_c[i] = np.sum( self.dNH1[0:i] ) + 0.5 * self.dNH1[i] self.NHe1_c[i] = np.sum( self.dNHe1[0:i] ) + 0.5 * self.dNHe1[i] self.NHe2_c[i] = np.sum( self.dNHe2[0:i] ) + 0.5 * self.dNHe2[i] # calculate column densities through whole slab #----------------------------------------------------------- self.NH1_thru = np.sum( self.dNH * self.xH1 ) self.logNH1_thru = np.log10( self.NH1_thru.magnitude ) self.NHe1_thru = np.sum( self.dNHe * self.xHe1 ) self.logNHe1_thru = np.log10( self.NHe1_thru.magnitude ) self.NHe2_thru = np.sum( self.dNHe * self.xHe2 ) self.logNHe2_thru = np.log10( self.NHe2_thru.magnitude ) # set preferred units #----------------------------------------------------------- self.z_c.units = 'kpc' # delete things we don't need #----------------------------------------------------------- del(self.kchem) del(self.kcool) # calculate effective photoionization rates # (i.e. the photoionization rates that result in # the same ionization fractions if caseA rates # are used). Note, we assume that the recombination # radiation is peaked near the ionization thresholds # and therefore does not alter the temperature. #------------------------------------------------ kchemA = chemistry.ChemistryRates( self.T, fcA_H2 = 1.0, fcA_He2 = 1.0, fcA_He3 = 1.0, H1i = self.H1i, He1i = self.He1i, He2i = self.He2i, fit_name = self.atomic_fit_name ) self.H1i_eff = kchemA.reH2 * self.ne * self.xH2 / self.xH1 - \ kchemA.ciH1 * self.ne self.He1i_eff = kchemA.reHe2 * self.ne * self.xHe2 / self.xHe1 - \ kchemA.ciHe1 * self.ne self.He2i_eff = kchemA.reHe3 * self.ne * self.xHe3 / self.xHe2 - \ kchemA.ciHe2 * self.ne