import numpy as np
from yt.utilities.on_demand_imports import _h5py as h5py
[docs]
def h5rd(fname, path, dtype=None):
"""Read Data. Return a dataset located at <path> in file <fname> as
a numpy array.
e.g. rd( fname, '/PartType0/Coordinates' )."""
data = None
fid = h5py.h5f.open(fname.encode("latin-1"), h5py.h5f.ACC_RDONLY)
dg = h5py.h5d.open(fid, path.encode("ascii"))
if dtype is None:
dtype = dg.dtype
data = np.zeros(dg.shape, dtype=dtype)
dg.read(h5py.h5s.ALL, h5py.h5s.ALL, data)
fid.close()
return data
[docs]
class IonTableSpectrum:
"""A class to handle the HM01 spectra in the OWLS ionization tables."""
def __init__(self, ion_file):
where = "/header/spectrum/gammahi"
self.GH1 = h5rd(ion_file, where) # GH1[1/s]
where = "/header/spectrum/logenergy_ryd"
self.logryd = h5rd(ion_file, where) # E[ryd]
where = "/header/spectrum/logflux"
self.logflux = h5rd(ion_file, where) # J[ergs/s/Hz/Sr/cm^2]
where = "/header/spectrum/redshift"
self.z = h5rd(ion_file, where) # z
[docs]
def return_table_GH1_at_z(self, z):
# find redshift indices
# -----------------------------------------------------------------
i_zlo = np.argmin(np.abs(self.z - z))
if self.z[i_zlo] < z:
i_zhi = i_zlo + 1
else:
i_zhi = i_zlo
i_zlo = i_zlo - 1
z_frac = (z - self.z[i_zlo]) / (self.z[i_zhi] - self.z[i_zlo])
# find GH1 from table
# -----------------------------------------------------------------
logGH1_all = np.log10(self.GH1)
dlog_GH1 = logGH1_all[i_zhi] - logGH1_all[i_zlo]
logGH1_table = logGH1_all[i_zlo] + z_frac * dlog_GH1
GH1_table = 10.0**logGH1_table
return GH1_table
[docs]
class IonTableOWLS:
"""A class to handle OWLS ionization tables."""
DELTA_nH = 0.25
DELTA_T = 0.1
def __init__(self, ion_file):
self.ion_file = ion_file
# ionbal is indexed like [nH, T, z]
# nH and T are log quantities
# ---------------------------------------------------------------
self.nH = h5rd(ion_file, "/logd") # log nH [cm^-3]
self.T = h5rd(ion_file, "/logt") # log T [K]
self.z = h5rd(ion_file, "/redshift") # z
# read the ionization fractions
# linear values stored in file so take log here
# ionbal is the ionization balance (i.e. fraction)
# ---------------------------------------------------------------
self.ionbal = h5rd(ion_file, "/ionbal").astype(np.float64)
self.ionbal_orig = self.ionbal.copy()
ipositive = self.ionbal > 0.0
izero = np.logical_not(ipositive)
self.ionbal[izero] = self.ionbal[ipositive].min()
self.ionbal = np.log10(self.ionbal)
# load in background spectrum
# ---------------------------------------------------------------
self.spectrum = IonTableSpectrum(ion_file)
# calculate the spacing along each dimension
# ---------------------------------------------------------------
self.dnH = self.nH[1:] - self.nH[0:-1]
self.dT = self.T[1:] - self.T[0:-1]
self.dz = self.z[1:] - self.z[0:-1]
self.order_str = "[log nH, log T, z]"
# sets iz and fz
# -----------------------------------------------------
[docs]
def set_iz(self, z):
if z <= self.z[0]:
self.iz = 0
self.fz = 0.0
elif z >= self.z[-1]:
self.iz = len(self.z) - 2
self.fz = 1.0
else:
for iz in range(len(self.z) - 1):
if z < self.z[iz + 1]:
self.iz = iz
self.fz = (z - self.z[iz]) / self.dz[iz]
break
# interpolate the table at a fixed redshift for the input
# values of nH and T ( input should be log ). A simple
# tri-linear interpolation is used.
# -----------------------------------------------------
[docs]
def interp(self, nH, T):
nH = np.array(nH)
T = np.array(T)
if nH.size != T.size:
raise ValueError(" owls_ion_tables: array size mismatch !!! ")
# field discovery will have nH.size == 1 and T.size == 1
# in that case we simply return 1.0
if nH.size == 1 and T.size == 1:
ionfrac = 1.0
return ionfrac
# find inH and fnH
# -----------------------------------------------------
x_nH = (nH - self.nH[0]) / self.DELTA_nH
x_nH_clip = np.clip(x_nH, 0.0, self.nH.size - 1.001)
fnH, inH = np.modf(x_nH_clip)
inH = inH.astype(np.int32)
# find iT and fT
# -----------------------------------------------------
x_T = (T - self.T[0]) / self.DELTA_T
x_T_clip = np.clip(x_T, 0.0, self.T.size - 1.001)
fT, iT = np.modf(x_T_clip)
iT = iT.astype(np.int32)
# short names for previously calculated iz and fz
# -----------------------------------------------------
iz = self.iz
fz = self.fz
# calculate interpolated value
# use tri-linear interpolation on the log values
# -----------------------------------------------------
ionfrac = (
self.ionbal[inH, iT, iz] * (1 - fnH) * (1 - fT) * (1 - fz)
+ self.ionbal[inH + 1, iT, iz] * (fnH) * (1 - fT) * (1 - fz)
+ self.ionbal[inH, iT + 1, iz] * (1 - fnH) * (fT) * (1 - fz)
+ self.ionbal[inH, iT, iz + 1] * (1 - fnH) * (1 - fT) * (fz)
+ self.ionbal[inH + 1, iT, iz + 1] * (fnH) * (1 - fT) * (fz)
+ self.ionbal[inH, iT + 1, iz + 1] * (1 - fnH) * (fT) * (fz)
+ self.ionbal[inH + 1, iT + 1, iz] * (fnH) * (fT) * (1 - fz)
+ self.ionbal[inH + 1, iT + 1, iz + 1] * (fnH) * (fT) * (fz)
)
return 10**ionfrac