From 70a25c5c2fe96219f2624ff592f97af7cb09fb24 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 30 Nov 2018 16:56:48 +0100
Subject: [PATCH] Applied python code formatting scripts.
---
examples/CoolingRates/plot_cooling_rates.py | 62 +++++----
.../SmallCosmoVolume_with_cooling/plotRhoT.py | 92 +++++++------
.../plotTempEvolution.py | 128 +++++++++++-------
3 files changed, 158 insertions(+), 124 deletions(-)
diff --git a/examples/CoolingRates/plot_cooling_rates.py b/examples/CoolingRates/plot_cooling_rates.py
index 6420835800..cca12468be 100644
--- a/examples/CoolingRates/plot_cooling_rates.py
+++ b/examples/CoolingRates/plot_cooling_rates.py
@@ -1,6 +1,6 @@
# Plots contribution to cooling rates from each of the different metals
# based on cooling_output.dat and cooling_element_*.dat files produced
-# by testCooling.
+# by testCooling.
import matplotlib.pyplot as plt
import numpy as np
@@ -10,42 +10,44 @@ elements = 11
# Declare arrays of internal energy and cooling rate
u = []
-cooling_rate = [[] for i in range(elements+1)]
-Temperature = [[] for i in range(elements+1)]
+cooling_rate = [[] for i in range(elements + 1)]
+Temperature = [[] for i in range(elements + 1)]
# Read in total cooling rate
-file_in = open('cooling_output.dat', 'r')
+file_in = open("cooling_output.dat", "r")
for line in file_in:
- data = line.split()
- u.append(float(data[0]))
- cooling_rate[0].append(-float(data[1]))
+ data = line.split()
+ u.append(float(data[0]))
+ cooling_rate[0].append(-float(data[1]))
file_in.close()
# Read in contributions to cooling rates from each of the elements
for elem in range(elements):
- file_in = open('cooling_element_'+str(elem)+'.dat','r')
- for line in file_in:
- data = line.split()
- cooling_rate[elem+1].append(-float(data[0]))
- file_in.close()
+ file_in = open("cooling_element_" + str(elem) + ".dat", "r")
+ for line in file_in:
+ data = line.split()
+ cooling_rate[elem + 1].append(-float(data[0]))
+ file_in.close()
# Plot
ax = plt.subplot(111)
-p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(u, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-p2, = plt.loglog(u, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'C')
-p3, = plt.loglog(u, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'N')
-p4, = plt.loglog(u, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'O')
-p5, = plt.loglog(u, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Ne')
-p6, = plt.loglog(u, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Mg')
-p7, = plt.loglog(u, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Si')
-p8, = plt.loglog(u, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'S')
-p9, = plt.loglog(u, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Ca')
-p10, = plt.loglog(u, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Fe')
-ax.set_position([0.15,0.15,0.75,0.75])
-plt.xlim([1e3,1e8])
-plt.ylim([1e-24,1e-21])
-plt.xlabel("Temperature ${\\rm{[K]}}$", fontsize = 14)
-plt.ylabel("${\Lambda/n_H^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
-plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-plt.savefig("cooling_rates",dpi=200)
+p0, = plt.loglog(u, cooling_rate[0], linewidth=0.5, color="k", label="Total")
+p1, = plt.loglog(
+ u, cooling_rate[1], linewidth=0.5, color="k", linestyle="--", label="H + He"
+)
+p2, = plt.loglog(u, cooling_rate[3], linewidth=0.5, color="b", label="C")
+p3, = plt.loglog(u, cooling_rate[4], linewidth=0.5, color="g", label="N")
+p4, = plt.loglog(u, cooling_rate[5], linewidth=0.5, color="r", label="O")
+p5, = plt.loglog(u, cooling_rate[6], linewidth=0.5, color="c", label="Ne")
+p6, = plt.loglog(u, cooling_rate[7], linewidth=0.5, color="m", label="Mg")
+p7, = plt.loglog(u, cooling_rate[8], linewidth=0.5, color="y", label="Si")
+p8, = plt.loglog(u, cooling_rate[9], linewidth=0.5, color="lightgray", label="S")
+p9, = plt.loglog(u, cooling_rate[10], linewidth=0.5, color="olive", label="Ca")
+p10, = plt.loglog(u, cooling_rate[11], linewidth=0.5, color="saddlebrown", label="Fe")
+ax.set_position([0.15, 0.15, 0.75, 0.75])
+plt.xlim([1e3, 1e8])
+plt.ylim([1e-24, 1e-21])
+plt.xlabel("Temperature ${\\rm{[K]}}$", fontsize=14)
+plt.ylabel("${\Lambda/n_H^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize=14)
+plt.legend(handles=[p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10])
+plt.savefig("cooling_rates", dpi=200)
diff --git a/examples/SmallCosmoVolume_with_cooling/plotRhoT.py b/examples/SmallCosmoVolume_with_cooling/plotRhoT.py
index 13f20652b6..4ba8ad66da 100644
--- a/examples/SmallCosmoVolume_with_cooling/plotRhoT.py
+++ b/examples/SmallCosmoVolume_with_cooling/plotRhoT.py
@@ -24,37 +24,39 @@ k_in_J_K = 1.38064852e-23
mH_in_kg = 1.6737236e-27
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.15,
-'figure.subplot.right' : 0.99,
-'figure.subplot.bottom' : 0.13,
-'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.15,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.13,
+ "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"]})
snap = int(sys.argv[1])
# Read the simulation data
-sim = h5py.File("snap_%04d.hdf5"%snap, "r")
+sim = h5py.File("snap_%04d.hdf5" % snap, "r")
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
z = sim["/Cosmology"].attrs["Redshift"][0]
@@ -65,7 +67,9 @@ 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,40 +89,40 @@ 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
+ return ret
rho = sim["/PartType0/Density"][:]
u = sim["/PartType0/InternalEnergy"][:]
# Compute the temperature
-u *= (unit_length_in_si**2 / unit_time_in_si**2)
-u /= a**(3 * (gas_gamma - 1.))
+u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
+u /= a ** (3 * (gas_gamma - 1.0))
Temp = T(u)
# Compute the physical density
-rho *= (unit_mass_in_cgs / unit_length_in_cgs**3)
-rho /= a**3
+rho *= unit_mass_in_cgs / unit_length_in_cgs ** 3
+rho /= a ** 3
rho /= mH_in_kg
# Life is better in log-space
@@ -134,7 +138,7 @@ log_T_max = 8
bins_x = np.linspace(log_rho_min, log_rho_max, 54)
bins_y = np.linspace(log_T_min, log_T_max, 54)
-H,_,_ = histogram2d(log_rho, log_T, bins=[bins_x, bins_y], normed=True)
+H, _, _ = histogram2d(log_rho, log_T, bins=[bins_x, bins_y], normed=True)
# Plot the interesting quantities
@@ -142,16 +146,18 @@ figure()
pcolormesh(bins_x, bins_y, np.log10(H).T)
-text(-5, 8., "$z=%.2f$"%z)
+text(-5, 8.0, "$z=%.2f$" % z)
-xticks([-5, -4, -3, -2, -1, 0, 1, 2, 3],
- ["", "$10^{-4}$", "", "$10^{-2}$", "", "$10^0$", "", "$10^2$", ""])
-yticks([2, 3, 4, 5, 6, 7, 8],
- ["$10^{2}$", "", "$10^{4}$", "", "$10^{6}$", "", "$10^8$"])
+xticks(
+ [-5, -4, -3, -2, -1, 0, 1, 2, 3],
+ ["", "$10^{-4}$", "", "$10^{-2}$", "", "$10^0$", "", "$10^2$", ""],
+)
+yticks(
+ [2, 3, 4, 5, 6, 7, 8], ["$10^{2}$", "", "$10^{4}$", "", "$10^{6}$", "", "$10^8$"]
+)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
xlim(-5.2, 3.2)
ylim(1, 8.5)
-savefig("rhoT_%04d.png"%snap, dpi=200)
-
+savefig("rhoT_%04d.png" % snap, dpi=200)
diff --git a/examples/SmallCosmoVolume_with_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume_with_cooling/plotTempEvolution.py
index aa6c5df5fe..66057449ce 100644
--- a/examples/SmallCosmoVolume_with_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume_with_cooling/plotTempEvolution.py
@@ -27,32 +27,34 @@ mH_in_kg = 1.6737236e-27
n_snapshots = 200
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("snap_0000.hdf5", "r")
@@ -64,7 +66,9 @@ 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"]
@@ -84,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
@@ -118,7 +123,7 @@ T_max = np.zeros(n_snapshots)
# Loop over all the snapshots
for i in range(n_snapshots):
- sim = h5py.File("snap_%04d.hdf5"%i, "r")
+ sim = h5py.File("snap_%04d.hdf5" % i, "r")
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
@@ -126,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
@@ -141,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()
@@ -179,4 +206,3 @@ xlim(9e-3, 1.1)
ylim(20, 2.5e7)
savefig("Temperature_evolution.png", dpi=200)
-
--
GitLab