diff --git a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
index 90de6cb712744359dbdfbf07cc4ed81546ea38bf..dab4b2c90a7b751c8d143ed38c614473c951988a 100644
--- a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
+++ b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
@@ -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)
-