Commit a862327b authored by Josh Borrow's avatar Josh Borrow
Browse files

Fixed unused variable and formatted

parent c4838d9b
......@@ -28,32 +28,34 @@ n_snapshots = 153
snapname = "santabarbara"
import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
import os.path
# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 9,
'legend.fontsize': 9,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize' : (3.15,3.15),
'figure.subplot.left' : 0.14,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.12,
'figure.subplot.top' : 0.99,
'figure.subplot.wspace' : 0.15,
'figure.subplot.hspace' : 0.12,
'lines.markersize' : 6,
'lines.linewidth' : 2.,
'text.latex.unicode': True
params = {
"axes.labelsize": 10,
"axes.titlesize": 10,
"font.size": 9,
"legend.fontsize": 9,
"xtick.labelsize": 10,
"ytick.labelsize": 10,
"text.usetex": True,
"figure.figsize": (3.15, 3.15),
"figure.subplot.left": 0.14,
"figure.subplot.right": 0.99,
"figure.subplot.bottom": 0.12,
"figure.subplot.top": 0.99,
"figure.subplot.wspace": 0.15,
"figure.subplot.hspace": 0.12,
"lines.markersize": 6,
"lines.linewidth": 2.0,
"text.latex.unicode": True,
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
# Read the simulation data
sim = h5py.File("%s_0000.hdf5" % snapname, "r")
......@@ -63,9 +65,10 @@ scheme = sim["/HydroScheme"].attrs["Scheme"][0]
kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
H_transition_temp = sim["/HydroScheme"].attrs[
"Hydrogen ionization transition temperature"
][0]
T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
git = sim["Code"].attrs["Git Revision"]
......@@ -85,25 +88,26 @@ unit_time_in_si = unit_time_in_cgs
# Primoridal ean molecular weight as a function of temperature
def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
if T > T_trans:
return 4. / (8. - 5. * (1. - H_frac))
return 4.0 / (8.0 - 5.0 * (1.0 - H_frac))
else:
return 4. / (1. + 3. * H_frac)
return 4.0 / (1.0 + 3.0 * H_frac)
# Temperature of some primoridal gas with a given internal energy
def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
T_over_mu = (gas_gamma - 1.0) * u * mH_in_kg / k_in_J_K
ret = np.ones(np.size(u)) * T_trans
# Enough energy to be ionized?
mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1, H_frac, T_trans))
if np.sum(mask_ionized) > 0:
ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10, H_frac, T_trans)
mask_ionized = T_over_mu > (T_trans + 1) / mu(T_trans + 1, H_frac, T_trans)
if np.sum(mask_ionized) > 0:
ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans * 10, H_frac, T_trans)
# Neutral gas?
mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1), H_frac, T_trans))
if np.sum(mask_neutral) > 0:
mask_neutral = T_over_mu < (T_trans - 1) / mu((T_trans - 1), H_frac, T_trans)
if np.sum(mask_neutral) > 0:
ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0, H_frac, T_trans)
return ret
......@@ -119,7 +123,7 @@ T_max = np.zeros(n_snapshots)
# Loop over all the snapshots
for i in range(n_snapshots):
sim = h5py.File("%s_%04d.hdf5"% (snapname, i), "r")
sim = h5py.File("%s_%04d.hdf5" % (snapname, i), "r")
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
......@@ -127,8 +131,8 @@ for i in range(n_snapshots):
u = sim["/PartType0/InternalEnergy"][:]
# Compute the temperature
u *= (unit_length_in_si**2 / unit_time_in_si**2)
u /= a[i]**(3 * (gas_gamma - 1.))
u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
u /= a[i] ** (3 * (gas_gamma - 1.0))
Temp = T(u)
# Gather statistics
......@@ -142,34 +146,56 @@ for i in range(n_snapshots):
# CMB evolution
a_evol = np.logspace(-3, 0, 60)
T_cmb = (1. / a_evol)**2 * 2.72
T_cmb = (1.0 / a_evol) ** 2 * 2.72
# Plot the interesting quantities
figure()
subplot(111, xscale="log", yscale="log")
fill_between(a, T_mean-T_std, T_mean+T_std, color='C0', alpha=0.1)
plot(a, T_max, ls='-.', color='C0', lw=1., label="${\\rm max}~T$")
plot(a, T_min, ls=':', color='C0', lw=1., label="${\\rm min}~T$")
plot(a, T_mean, color='C0', label="${\\rm mean}~T$", lw=1.5)
fill_between(a, 10**(T_log_mean-T_log_std), 10**(T_log_mean+T_log_std), color='C1', alpha=0.1)
plot(a, 10**T_log_mean, color='C1', label="${\\rm mean}~{\\rm log} T$", lw=1.5)
plot(a, T_median, color='C2', label="${\\rm median}~T$", lw=1.5)
fill_between(a, T_mean - T_std, T_mean + T_std, color="C0", alpha=0.1)
plot(a, T_max, ls="-.", color="C0", lw=1.0, label="${\\rm max}~T$")
plot(a, T_min, ls=":", color="C0", lw=1.0, label="${\\rm min}~T$")
plot(a, T_mean, color="C0", label="${\\rm mean}~T$", lw=1.5)
fill_between(
a,
10 ** (T_log_mean - T_log_std),
10 ** (T_log_mean + T_log_std),
color="C1",
alpha=0.1,
)
plot(a, 10 ** T_log_mean, color="C1", label="${\\rm mean}~{\\rm log} T$", lw=1.5)
plot(a, T_median, color="C2", label="${\\rm median}~T$", lw=1.5)
legend(loc="upper left", frameon=False, handlelength=1.5)
# Expected lines
plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], 'k--', lw=0.5, alpha=0.7)
text(2.5e-2, H_transition_temp*1.07, "$T_{\\rm HII\\rightarrow HI}$", va="bottom", alpha=0.7, fontsize=8)
plot([1e-10, 1e10], [T_minimal, T_minimal], 'k--', lw=0.5, alpha=0.7)
text(1e-2, T_minimal*0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
plot(a_evol, T_cmb, 'k--', lw=0.5, alpha=0.7)
text(a_evol[20], T_cmb[20]*0.55, "$(1+z)^2\\times T_{\\rm CMB,0}$", rotation=-34, alpha=0.7, fontsize=8, va="top", bbox=dict(facecolor='w', edgecolor='none', pad=1.0, alpha=0.9))
redshift_ticks = np.array([0., 1., 2., 5., 10., 20., 50., 100.])
plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7)
text(
2.5e-2,
H_transition_temp * 1.07,
"$T_{\\rm HII\\rightarrow HI}$",
va="bottom",
alpha=0.7,
fontsize=8,
)
plot([1e-10, 1e10], [T_minimal, T_minimal], "k--", lw=0.5, alpha=0.7)
text(1e-2, T_minimal * 0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
plot(a_evol, T_cmb, "k--", lw=0.5, alpha=0.7)
text(
a_evol[20],
T_cmb[20] * 0.55,
"$(1+z)^2\\times T_{\\rm CMB,0}$",
rotation=-34,
alpha=0.7,
fontsize=8,
va="top",
bbox=dict(facecolor="w", edgecolor="none", pad=1.0, alpha=0.9),
)
redshift_ticks = np.array([0.0, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0])
redshift_labels = ["$0$", "$1$", "$2$", "$5$", "$10$", "$20$", "$50$", "$100$"]
a_ticks = 1. / (redshift_ticks + 1.)
a_ticks = 1.0 / (redshift_ticks + 1.0)
xticks(a_ticks, redshift_labels)
minorticks_off()
......@@ -180,4 +206,3 @@ xlim(9e-3, 1.1)
ylim(20, 2.5e7)
savefig("Temperature_evolution.png", dpi=200)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment