Source code for rabacus.f2py.slab_plane

""" Wrapper to make the fortran calling identical to the python calling. """ 

import rabacus_fc
import numpy as np
from slab_base import SlabBase
from rabacus.constants import units
from rabacus.utils import utils
from rabacus.atomic import chemistry
from rabacus.atomic import cooling



__all__ = ['SlabPln', 'Slab2Pln']



[docs]class SlabPln(SlabBase): """ 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", "ray"} `fixed_fcA` (float): If `rec_meth` = "fixed", constant caseA fraction `atomic_fit_name` (string): Source for atomic rate fits {"hg97", "enzo13"} `find_Teq` (bool): If ``False``, use fixed input T, if ``True`` solve for equilibrium T `z` (float): Redshift, only need if `find_Teq` = ``True`` `Hz` (float): Hubble parameter at `z` if Hubble cooling is desired. `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 C to center of layer `dz` (array): radial thickness of layer `NH_c` (array): H column density from C to center of layer `dNH` (array): H column density through layer `NH1_c` (array): HI column density from C 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 photo-heating `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 from C to R `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`. """ def __init__(self, Edges, T, nH, nHe, rad_src, # rec_meth = "fixed", fixed_fcA = 1.0, atomic_fit_name = "hg97", find_Teq = False, z = None, Hz = None, verbose = False, tol = 1.0e-4, thin = False, ): # 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 self.fixed_fcA = fixed_fcA 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 else: self.z = 0 # call base class init #----------------------------------------------- super(SlabPln, self).__init__() if Hz == None: self.Hz = 0.0 / self.U.s else: self.Hz = Hz self.Hz.units = '1/s' # call the fortran solver routine #----------------------------------------------- self.i_2side = 0 (xH1, xH2, xHe1, xHe2, xHe3, H1i_src, He1i_src, He2i_src, H1i_rec, He1i_rec, He2i_rec, H1h_src, He1h_src, He2h_src, H1h_rec, He1h_rec, He2h_rec) = \ rabacus_fc.slab_plane.slab_plane_solve( self.Edges, self.nH, self.nHe, self.T, self.E_eV, self.shape, self.i_2side, self.i_rec_meth, self.fixed_fcA, self.i_photo_fit, self.i_rate_fit, self.i_find_Teq, self.i_thin, self.z, self.Hz, self.tol, self.Nl, self.Nnu ) # attach output with units. #----------------------------------------------- self.xH1 = xH1 self.xH2 = xH2 self.xHe1 = xHe1 self.xHe2 = xHe2 self.xHe3 = xHe3 self.H1i_src = H1i_src / self.U.s self.He1i_src = He1i_src / self.U.s self.He2i_src = He2i_src / self.U.s self.H1i_rec = H1i_rec / self.U.s self.He1i_rec = He1i_rec / self.U.s self.He2i_rec = He2i_rec / self.U.s self.H1i = self.H1i_src + self.H1i_rec self.He1i = self.He1i_src + self.He1i_rec self.He2i = self.He2i_src + self.He2i_rec self.H1h_src = H1h_src * self.U.erg / self.U.s self.He1h_src = He1h_src * self.U.erg / self.U.s self.He2h_src = He2h_src * self.U.erg / self.U.s self.H1h_rec = H1h_rec * self.U.erg / self.U.s self.He1h_rec = He1h_rec * self.U.erg / self.U.s self.He2h_rec = He2h_rec * self.U.erg / self.U.s self.H1h = self.H1h_src + self.H1h_rec self.He1h = self.He1h_src + self.He1h_rec self.He2h = self.He2h_src + self.He2h_rec # do post-calculation #----------------------------------------------- super(SlabPln,self).__post__()
[docs]class Slab2Pln(SlabBase): """ Solves a slab geometry with plane parallel radiation incident from both sides. 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", "ray"} `fixed_fcA` (float): If `rec_meth` = "fixed", constant caseA fraction `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`` `Hz` (float): Hubble parameter at `z` if Hubble cooling is desired. `verbose` (bool): Verbose output? `tol` (float): tolerance for all convergence tests `thin` (bool): if ``True`` only solves optically thin `i_rec_ems_prf` (int): if ``0`` recombination radiation is emitted isotropically (i.e. more attenuation), if ``1`` recombination radiation is emitted perpendicular to the layers. 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`. """ def __init__(self, Edges, T, nH, nHe, rad_src, # rec_meth = "fixed", fixed_fcA = 1.0, atomic_fit_name = "hg97", find_Teq = False, z = None, Hz = None, verbose = False, tol = 1.0e-4, thin = False, i_rec_ems_prf = 0 ): # initialize #----------------------------------------------- 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 self.fixed_fcA = fixed_fcA self.atomic_fit_name = atomic_fit_name self.find_Teq = find_Teq self.verbose = verbose self.tol = tol self.thin = thin self.i_rec_ems_prf = i_rec_ems_prf if find_Teq: self.z = z else: self.z = 0 # call base class init #----------------------------------------------- super(Slab2Pln,self).__init__() if Hz == None: self.Hz = 0.0 / self.U.s else: self.Hz = Hz self.Hz.units = '1/s' # check for even number of layers #----------------------------------------------- if self.nH.size % 2 != 0: msg = 'the number of layers must be even' raise utils.InputError(msg) # call the fortran solver routine #----------------------------------------------- self.i_2side = 1 (xH1, xH2, xHe1, xHe2, xHe3, H1i_src, He1i_src, He2i_src, H1i_rec, He1i_rec, He2i_rec, H1h_src, He1h_src, He2h_src, H1h_rec, He1h_rec, He2h_rec) = \ rabacus_fc.slab_plane.slab_plane_solve( self.Edges, self.nH, self.nHe, self.T, self.E_eV, self.shape, self.i_2side, self.i_rec_meth, self.fixed_fcA, self.i_photo_fit, self.i_rate_fit, self.i_find_Teq, self.i_thin, self.z, self.Hz, self.tol, self.Nl, self.Nnu ) # attach output with units. #----------------------------------------------- self.xH1 = xH1 self.xH2 = xH2 self.xHe1 = xHe1 self.xHe2 = xHe2 self.xHe3 = xHe3 self.H1i_src = H1i_src / self.U.s self.He1i_src = He1i_src / self.U.s self.He2i_src = He2i_src / self.U.s self.H1i_rec = H1i_rec / self.U.s self.He1i_rec = He1i_rec / self.U.s self.He2i_rec = He2i_rec / self.U.s self.H1i = self.H1i_src + self.H1i_rec self.He1i = self.He1i_src + self.He1i_rec self.He2i = self.He2i_src + self.He2i_rec self.H1h_src = H1h_src * self.U.erg / self.U.s self.He1h_src = He1h_src * self.U.erg / self.U.s self.He2h_src = He2h_src * self.U.erg / self.U.s self.H1h_rec = H1h_rec * self.U.erg / self.U.s self.He1h_rec = He1h_rec * self.U.erg / self.U.s self.He2h_rec = He2h_rec * self.U.erg / self.U.s self.H1h = self.H1h_src + self.H1h_rec self.He1h = self.He1h_src + self.He1h_rec self.He2h = self.He2h_src + self.He2h_rec # do post-calculation #----------------------------------------------- super(Slab2Pln,self).__post__()