Author: Britton Smith
Modified by : Andrew Emerick
In [1]:
import h5py
from matplotlib import pyplot
%matplotlib inline
import numpy as np
LW_model = 'Qin2020'
In [2]:
pyplot.rcParams['figure.figsize'] = (10, 6)
pyplot.rcParams['font.size'] = 14
In [3]:
from make_table import k31_JW2012, k31_RFT14, k31_Qin2020
In [4]:
def interpvals(xtab, ytab, x, log=True):
i = np.digitize(x, xtab) - 1
i = np.clip(i, 0, xtab.size-2)
if log:
m = np.log10(ytab[i+1] / ytab[i]) / np.log10(xtab[i+1] / xtab[i])
return np.power(10, m * np.log10(x / xtab[i]) + np.log10(ytab[i]))
else:
m = (ytab[i+1] - ytab[i]) / (xtab[i+1] - xtab[i])
return m * (x - xtab[i]) + ytab[i]
In [5]:
def load_rates(filename, group, rates, zero_val=1e-50):
print ("Loading rates from %s: %s" % (filename, rates))
data = {}
fh = h5py.File(filename, 'r')
g = fh['UVBRates']
data['z'] = g['z'].value
for rate in rates:
data[rate] = g[group][rate].value
data[rate] = np.clip(data[rate], 1e-50, np.inf)
fh.close()
return data
In [6]:
def plot_rates(z, filenames, group, rates):
pyplot.xscale('log')
pyplot.yscale('log')
lss = ['-', '--', ':']
cmap = pyplot.cm.jet
tdata = dict((filename, load_rates(filename, group, rates))
for filename in filenames)
for ir, rate in enumerate(rates):
for ifn, fn in enumerate(filenames):
ztab = tdata[fn]['z']
my_rate = interpvals(ztab+1, tdata[fn][rate], z+1)
if ifn == 0:
label = rate
else:
label = None
pyplot.plot(z+1, my_rate, linestyle=lss[ifn],
color=cmap((ir+1)/len(rates)),
label=label)
pyplot.xlim(z[0]+1, z[-1]+1)
pyplot.xlabel('z+1')
pyplot.ylabel('rates [CGS]')
pyplot.legend(loc='best')
In [7]:
z = np.linspace(0, 50, 100)
plot_rates(z, ['CloudyData_UVB=HM2012.h5',
'CloudyData_HM2012_highz.h5'],
'Chemistry', ['k24', 'k25', 'k26', 'k29', 'k30'])
pyplot.ylim(1e-29, 1e-11)
Out[7]:
In [8]:
z = np.linspace(0, 50, 100)
plot_rates(z, ['CloudyData_UVB=HM2012.h5',
'CloudyData_HM2012_highz.h5'],
'Photoheating', ['piHI', 'piHeI', 'piHeII'])
pyplot.ylim(1e-26, 1e-11)
Out[8]:
In [9]:
z = np.linspace(0, 50, 100)
pyplot.plot(z, k31_RFT14(z), color='red', label='k31 JHW')
pyplot.plot(z, k31_JW2012(z), color='red', ls = ':', label='k31 JW2012')
pyplot.plot(z, k31_Qin2020(z), color='red', ls = '-.', label='k31 Qin2020')
plot_rates(z, ['CloudyData_UVB=HM2012.h5',
'CloudyData_HM2012_highz.h5'],
'Chemistry', ['k27', 'k28', 'k31'])
pyplot.ylim(1e-19, 1e-7)
Out[9]:
In [10]:
z = np.linspace(0, 50, 100)
plot_rates(z, ['CloudyData_UVB=HM2012_shielded.h5',
'CloudyData_HM2012_highz_shielded.h5'],
'Chemistry', ['k24', 'k25', 'k26', 'k29', 'k30'])
pyplot.ylim(1e-29, 1e-11)
Out[10]:
In [11]:
z = np.linspace(0, 50, 100)
plot_rates(z, ['CloudyData_UVB=HM2012_shielded.h5',
'CloudyData_HM2012_highz_shielded.h5'],
'Photoheating', ['piHI', 'piHeI', 'piHeII'])
pyplot.ylim(1e-26, 1e-11)
Out[11]:
In [13]:
z = np.linspace(0, 50, 100)
plot_rates(z, ['CloudyData_UVB=HM2012_shielded.h5',
'CloudyData_HM2012_highz_shielded.h5'],
'Chemistry', ['k27', 'k28', 'k31'])
pyplot.ylim(1e-19, 1e-7)
pyplot.plot(z, k31_RFT14(z), color='red', label='k31 JHW')
pyplot.plot(z, k31_JW2012(z), color='red', ls = ':', label='k31 JW2012')
pyplot.plot(z, k31_Qin2020(z), color='red', ls = '-.', label='k31 Qin2020')
Out[13]:
In [ ]: