diff --git a/.gitignore b/.gitignore
index c50b672a79a52ca588ca30ada2ccfe9f3fc59ad1..5521af273e91ec3a786e1272b8d5bb75a8915d3e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -69,6 +69,10 @@ examples/Cooling/CoolingRates/cooling_output.dat
 examples/SubgridTests/StellarEvolution/StellarEvolutionSolution*
 examples/SubgridTests/CosmologicalStellarEvolution/StellarEvolutionSolution*
 examples/SmallCosmoVolume/SmallCosmoVolume_DM/power_spectra
+examples/SmallCosmoVolume/SmallCosmoVolume_cooling/snapshots/
+examples/SmallCosmoVolume/SmallCosmoVolume_hydro/snapshots/
+examples/SmallCosmoVolume/SmallCosmoVolume_cooling/CloudyData_UVB=HM2012.h5
+examples/SmallCosmoVolume/SmallCosmoVolume_cooling/CloudyData_UVB=HM2012_shielded.h5
 
 tests/testActivePair
 tests/testActivePair.sh
diff --git a/examples/Cooling/CoolingBox/plotEnergy.py b/examples/Cooling/CoolingBox/plotEnergy.py
index 29b9a2594798da0a589faf18bbc8dd0295917d0d..64d34c6127244e6ed920cc7f4b44676f034a8904 100644
--- a/examples/Cooling/CoolingBox/plotEnergy.py
+++ b/examples/Cooling/CoolingBox/plotEnergy.py
@@ -1,38 +1,35 @@
-from h5py import File
-import numpy as np
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 import matplotlib
-from glob import glob
 
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
+import numpy as np
+from glob import glob
+import h5py
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (5, 5),
-    "figure.subplot.left": 0.145,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.11,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-plt.rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Some constants in cgs units
 k_b_cgs = 1.38e-16  # boltzmann
 m_h_cgs = 1.67e-24  # proton mass
 
-
 # File containing the total energy
 stats_filename = "./statistics.txt"
 
@@ -40,7 +37,7 @@ stats_filename = "./statistics.txt"
 snap_filename = "coolingBox_0000.hdf5"
 
 # Read the initial state of the gas
-f = File(snap_filename, "r")
+f = h5py.File(snap_filename, "r")
 
 # Read the units parameters from the snapshot
 units = f["InternalCodeUnits"]
@@ -88,7 +85,7 @@ N = len(files)
 temp_snap = np.zeros(N)
 time_snap_cgs = np.zeros(N)
 for i in range(N):
-    snap = File(files[i], "r")
+    snap = h5py.File(files[i], "r")
     u = snap["/PartType0/InternalEnergies"][:] * snap["/PartType0/Masses"][:]
     u = sum(u) / total_mass[0]
     temp_snap[i] = energyUnits(u)
@@ -97,7 +94,7 @@ for i in range(N):
 
 plt.figure()
 
-Myr_in_s = 3.15e13
+Myr_in_s = 1e6 * 365.25 * 24.0 * 60.0 * 60.0
 plt.plot(time / Myr_in_s, total_energy_cgs, "r-", lw=1.6, label="Gas total energy")
 # statistics and snapshots may not be at same timestep and frequency
 plt.plot(time_snap_cgs / Myr_in_s, temp_snap, "rD", ms=3)
@@ -114,4 +111,6 @@ plt.legend(loc="right", fontsize=8, frameon=False, handlelength=3, ncol=1)
 plt.xlabel("${\\rm{Time~[Myr]}}$", labelpad=0)
 plt.ylabel("${\\rm{Internal ~Energy ~(u ~m_H / k_B) ~[K]}}$")
 
+plt.tight_layout()
+
 plt.savefig("energy.png", dpi=200)
diff --git a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
index 7f153b9ad39b6e877022f39b31481ab37df80993..cc25919e0ff30e32f2996f3f9205ca6d47d90833 100644
--- a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
+++ b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
@@ -25,6 +25,7 @@ T_i = 100.0  # Initial temperature of the gas (in K)
 z_c = 1.0  # Redshift of caustic formation (non-linear collapse)
 z_i = 100.0  # Initial redshift
 gas_gamma = 5.0 / 3.0  # Gas adiabatic index
+N_output = 119
 
 # Physical constants needed for internal energy to temperature conversion
 kB_in_SI = 1.38064852e-23
@@ -33,30 +34,11 @@ mH_in_kg = 1.6737236e-27
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
-import os.path
-
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.06,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.06,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.13,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Read the simulation data
 sim = h5py.File("box_0000.hdf5", "r")
@@ -85,25 +67,25 @@ u_0 *= a ** (-3 * (1 - gas_gamma))
 S_0 = (gas_gamma - 1.0) * u_0 * rho_0 ** (-(gas_gamma - 1.0))
 
 # Mean quantities over time
-z = np.zeros(119)
-a = np.zeros(119)
-S_mean = np.zeros(119)
-S_std = np.zeros(119)
-u_mean = np.zeros(119)
-u_std = np.zeros(119)
-P_mean = np.zeros(119)
-P_std = np.zeros(119)
-rho_mean = np.zeros(119)
-rho_std = np.zeros(119)
-
-vx_mean = np.zeros(119)
-vy_mean = np.zeros(119)
-vz_mean = np.zeros(119)
-vx_std = np.zeros(119)
-vy_std = np.zeros(119)
-vz_std = np.zeros(119)
-
-for i in range(119):
+z = np.zeros(N_output)
+a = np.zeros(N_output)
+S_mean = np.zeros(N_output)
+S_std = np.zeros(N_output)
+u_mean = np.zeros(N_output)
+u_std = np.zeros(N_output)
+P_mean = np.zeros(N_output)
+P_std = np.zeros(N_output)
+rho_mean = np.zeros(N_output)
+rho_std = np.zeros(N_output)
+
+vx_mean = np.zeros(N_output)
+vy_mean = np.zeros(N_output)
+vz_mean = np.zeros(N_output)
+vx_std = np.zeros(N_output)
+vy_std = np.zeros(N_output)
+vz_std = np.zeros(N_output)
+
+for i in range(N_output):
     sim = h5py.File("box_%04d.hdf5" % i, "r")
 
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
@@ -138,117 +120,177 @@ rho_mean_phys = rho_mean / a ** 3
 u_mean_phys = u_mean / a ** (3 * (gas_gamma - 1.0))
 S_mean_phys = S_mean
 
-# Solution in physical coordinates
-# T_solution = np.ones(T) / a
-
-figure()
+vx_mean_phys = vx_mean / a
+vy_mean_phys = vy_mean / a
+vz_mean_phys = vz_mean / a
+vx_std_phys = vx_std / a
+vy_std_phys = vy_std / a
+vz_std_phys = vz_std / a
+
+# Plot the interesting quantities
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=4,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
 
 # Density evolution --------------------------------
-subplot(231)  # , yscale="log")
-semilogx(a, rho_mean / rho_0, "-", color="r", lw=1)
-semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
-text(1e-2, 1.0105, "+1\\%", color="0.6", fontsize=9)
-text(1e-2, 0.9895, "-1\\%", color="0.6", fontsize=9, va="top")
-text(1e-2, 1.015, "$\\rho_0=%.3f$" % rho_0)
-ylim(0.98, 1.02)
-xlim(8e-3, 1.1)
-xlabel("${\\rm Scale-factor}$", labelpad=0.0)
-ylabel("${\\rm Comoving~density}~\\rho / \\rho_0$", labelpad=0.0)
+plt.subplot(231)  # , yscale="log")
+plt.semilogx(a, rho_mean / rho_0, **scatter_props)
+plt.semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
+plt.text(1e-2, 1.0105, "+1%", color="0.6", fontsize=9)
+plt.text(1e-2, 0.9895, "-1%", color="0.6", fontsize=9, va="top")
+plt.text(1e-2, 1.015, "$\\rho_0=%.3f$" % rho_0)
+plt.ylim(0.98, 1.02)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Comoving~density}~\\rho / \\rho_0$", labelpad=0.0)
 
 # Thermal energy evolution --------------------------------
-subplot(232)  # , yscale="log")
-semilogx(a, u_mean / u_0, "-", color="r", lw=1)
-semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
-text(1e-2, 1.0105, "+1\\%", color="0.6", fontsize=9)
-text(1e-2, 0.9895, "-1\\%", color="0.6", fontsize=9, va="top")
-text(1e-2, 1.015, "$u_0=%.3e$" % (u_0))
-ylim(0.98, 1.02)
-xlim(8e-3, 1.1)
-xlabel("${\\rm Scale-factor}$", labelpad=0.0)
-ylabel("${\\rm Comoving~internal~energy}~u / u_0$", labelpad=0.0)
+plt.subplot(232)  # , yscale="log")
+plt.semilogx(a, u_mean / u_0, **scatter_props)
+plt.semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
+plt.text(1e-2, 1.0105, "+1%", color="0.6", fontsize=9)
+plt.text(1e-2, 0.9895, "-1%", color="0.6", fontsize=9, va="top")
+plt.text(1e-2, 1.015, "$u_0=%.3e$" % (u_0))
+plt.ylim(0.98, 1.02)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Comoving~internal~energy}~u / u_0$", labelpad=0.0)
 
 # Entropy evolution --------------------------------
-subplot(233)  # , yscale="log")
-semilogx(a, S_mean / S_0, "-", color="r", lw=1)
-semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
-semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
-text(1e-2, 1.0105, "+1\\%", color="0.6", fontsize=9)
-text(1e-2, 0.9895, "-1\\%", color="0.6", fontsize=9, va="top")
-text(1e-2, 1.015, "$A_0=%.3e$" % (S_0))
-ylim(0.98, 1.02)
-xlim(8e-3, 1.1)
-xlabel("${\\rm Scale-factor}$", labelpad=0.0)
-ylabel("${\\rm Comoving~entropy}~A / A_0$", labelpad=0.0)
+plt.subplot(233)  # , yscale="log")
+plt.semilogx(a, S_mean / S_0, **scatter_props)
+plt.semilogx([1e-10, 1e10], np.ones(2), "-", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 0.99, "--", color="0.6", lw=1.0)
+plt.semilogx([1e-10, 1e10], np.ones(2) * 1.01, "--", color="0.6", lw=1.0)
+plt.text(1e-2, 1.0105, "+1%", color="0.6", fontsize=9)
+plt.text(1e-2, 0.9895, "-1%", color="0.6", fontsize=9, va="top")
+plt.text(1e-2, 1.015, "$A_0=%.3e$" % (S_0))
+plt.ylim(0.98, 1.02)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Comoving~entropy}~A / A_0$", labelpad=0.0)
 
 # Peculiar velocity evolution ---------------------
-subplot(234)
-semilogx(a, vx_mean, "-", color="r", lw=1)
-semilogx(a, vy_mean, "-", color="g", lw=1)
-semilogx(a, vz_mean, "-", color="b", lw=1)
-xlabel("${\\rm Scale-factor}$", labelpad=0.0)
-ylabel("${\\rm Peculiar~velocity~mean}$", labelpad=-5.0)
+plt.subplot(234)
+plt.semilogx(a, vx_mean, **scatter_props)
+plt.semilogx(a, vy_mean, **scatter_props)
+plt.semilogx(a, vz_mean, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Peculiar~velocity~mean}$", labelpad=-5.0)
 
 # Peculiar velocity evolution ---------------------
-subplot(235)
-semilogx(a, vx_std, "--", color="r", lw=1)
-semilogx(a, vy_std, "--", color="g", lw=1)
-semilogx(a, vz_std, "--", color="b", lw=1)
-xlabel("${\\rm Scale-factor}$", labelpad=0.0)
-ylabel("${\\rm Peculiar~velocity~std-dev}$", labelpad=0.0)
+plt.subplot(235)
+plt.semilogx(a, vx_std, **scatter_props)
+plt.semilogx(a, vy_std, **scatter_props)
+plt.semilogx(a, vz_std, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Peculiar~velocity~std-dev}$", labelpad=0.0)
 
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
 
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
+plt.tight_layout()
 
-savefig("ConstantBox_comoving.png", dpi=200)
+plt.savefig("ConstantBox_comoving.png", dpi=200)
 
 
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
 
 # Density evolution --------------------------------
-subplot(231)  # , yscale="log")
-loglog(a, rho_mean_phys, "-", color="r", lw=1)
-xlabel("${\\rm Scale-factor}$")
-ylabel("${\\rm Physical~density}$")
+plt.subplot(231)  # , yscale="log")
+plt.loglog(a, rho_mean_phys, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$")
+plt.ylabel("${\\rm Physical~density}$")
 
 # Thermal energy evolution --------------------------------
-subplot(232)  # , yscale="log")
-loglog(a, u_mean_phys, "-", color="r", lw=1)
-xlabel("${\\rm Scale-factor}$")
-ylabel("${\\rm Physical~internal~energy}$")
+plt.subplot(232)  # , yscale="log")
+plt.loglog(a, u_mean_phys, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$")
+plt.ylabel("${\\rm Physical~internal~energy}$")
 
 # Entropy evolution --------------------------------
-subplot(233)  # , yscale="log")
-semilogx(a, S_mean_phys, "-", color="r", lw=1)
-xlabel("${\\rm Scale-factor}$")
-ylabel("${\\rm Physical~entropy}$")
+plt.subplot(233)  # , yscale="log")
+plt.semilogx(a, S_mean_phys, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$")
+plt.ylabel("${\\rm Physical~entropy}$")
+
+# Peculiar velocity evolution ---------------------
+plt.subplot(234)
+plt.semilogx(a, vx_mean_phys, **scatter_props)
+plt.semilogx(a, vy_mean_phys, **scatter_props)
+plt.semilogx(a, vz_mean_phys, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Peculiar~velocity~mean}$", labelpad=-5.0)
+
+# Peculiar velocity evolution ---------------------
+plt.subplot(235)
+plt.semilogx(a, vx_std_phys, **scatter_props)
+plt.semilogx(a, vy_std_phys, **scatter_props)
+plt.semilogx(a, vz_std_phys, **scatter_props)
+plt.xlim(8e-3, 1.1)
+plt.xlabel("${\\rm Scale-factor}$", labelpad=0.0)
+plt.ylabel("${\\rm Peculiar~velocity~std-dev}$", labelpad=0.0)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-savefig("ConstantBox_physical.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("ConstantBox_physical.png", dpi=200)
diff --git a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
index 2fc0080cab691cf15460b370f8548e57516bf1de..293f93792a2cf4d7d0af86557a15603700df07e5 100644
--- a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
@@ -32,30 +32,13 @@ mH_in_kg = 1.6737236e-27
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 import os.path
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -84,7 +67,7 @@ P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 m = sim["/PartType0/Masses"][:]
 try:
-    phi = sim["/PartType0/Potential"][:]
+    phi = sim["/PartType0/Potentials"][:]
 except KeyError:
     # We didn't write the potential, try to go on without
     print("Couldn't find potential in your output file")
@@ -118,12 +101,12 @@ lambda_i = boxSize  # wavelength of the perturbation
 
 
 # Solution taken from Springel 2010. Eqs. 127 - 130
-q = linspace(-0.5 * lambda_i, 0.5 * lambda_i, 256)
-k_i = 2.0 * pi / lambda_i
+q = np.linspace(-0.5 * lambda_i, 0.5 * lambda_i, 256)
+k_i = 2.0 * np.pi / lambda_i
 zfac = (1.0 + z_c) / (1.0 + redshift)
-x_s = q - zfac * sin(k_i * q) / k_i
-rho_s = rho_0 / (1.0 - zfac * cos(k_i * q))
-v_s = -H_0 * (1.0 + z_c) / sqrt(1.0 + redshift) * sin(k_i * q) / k_i
+x_s = q - zfac * np.sin(k_i * q) / k_i
+rho_s = rho_0 / (1.0 - zfac * np.cos(k_i * q))
+v_s = -H_0 * (1.0 + z_c) / np.sqrt(1.0 + redshift) * np.sin(k_i * q) / k_i
 T_s = T_i * (((1.0 + redshift) / (1.0 + z_i)) ** 3.0 * rho_s / rho_0) ** (2.0 / 3.0)
 
 
@@ -136,37 +119,51 @@ unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=4,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
 
 # Velocity profile --------------------------------
-subplot(231)
+plt.subplot(231)
 if np.size(x_g) > 1:
-    plot(x_g, v_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
-plot(x, v, ".", color="r", ms=4.0)
-plot(x_s, v_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
-ylabel("${\\rm{Peculiar~velocity}}~v_x$", labelpad=0)
+    plt.plot(x_g, v_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
+plt.plot(x, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Peculiar~velocity}}~v_x$", labelpad=0)
 
 
 # Density profile --------------------------------
-subplot(232)  # , yscale="log")
+plt.subplot(232)  # , yscale="log")
 if np.size(x_g) > 1:
-    plot(x_g, rho_g / rho_0, "s", color="g", alpha=0.8, lw=1.2, ms=4)
-plot(x, rho / rho_0, ".", color="r", ms=4.0)
-plot(x_s, rho_s / rho_0, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho / \\rho_0$", labelpad=0)
+    plt.plot(x_g, rho_g / rho_0, "s", color="g", alpha=0.8, lw=1.2, ms=4)
+plt.plot(x, rho / rho_0, **scatter_props)
+plt.plot(x_s, rho_s / rho_0, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho / \\rho_0$", labelpad=0)
 
 # Potential profile --------------------------------
-subplot(233)
+plt.subplot(233)
 if np.size(x_g) > 1:
-    plot(x_g, phi_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
-plot(x, phi, ".", color="r", ms=4.0)
-xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
-ylabel("${\\rm{Potential}}~\\phi$", labelpad=0)
+    plt.plot(x_g, phi_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
+plt.plot(x, phi, **scatter_props)
+plt.xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Potential}}~\\phi$", labelpad=0)
 
 # Temperature profile -------------------------
-subplot(234)  # , yscale="log")
+plt.subplot(234)  # , yscale="log")
 u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
 u_g *= unit_length_in_si ** 2 / unit_time_in_si ** 2
 u /= a ** (3 * (gas_gamma - 1.0))
@@ -175,30 +172,41 @@ T = (gas_gamma - 1.0) * u * mH_in_kg / k_in_J_K
 T_g = (gas_gamma - 1.0) * u_g * mH_in_kg / k_in_J_K
 print("z = {0:.2f}, T_avg = {1:.2f}".format(redshift, T.mean()))
 if np.size(x_g) > 1:
-    plot(x_g, T_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
-plot(x, T, ".", color="r", ms=4.0)
-plot(x_s, T_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
-ylabel("${\\rm{Temperature}}~T$", labelpad=0)
+    plt.plot(x_g, T_g, "s", color="g", alpha=0.8, lw=1.2, ms=4)
+plt.plot(x, T, **scatter_props)
+plt.plot(x_s, T_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Comoving~position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Temperature}}~T$", labelpad=0)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
 
-text(-0.49, 0.9, "Zeldovich pancake at z=%.2f " % (redshift), fontsize=10)
-text(
-    -0.49,
+plt.text(
+    -0.45, 0.9, "Zeldovich pancake at z=%.2f " % (redshift), fontsize=text_fontsize
+)
+plt.text(
+    -0.45,
     0.8,
     "adiabatic index $\\gamma=%.2f$, viscosity $\\alpha=%.2f$" % (gas_gamma, alpha),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-savefig("ZeldovichPancake_%.4d.png" % snap, dpi=200)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("ZeldovichPancake_%.4d.png" % snap, dpi=200)
diff --git a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
index 6b89eb4ff04011454ece10c9190a06880ae6b21b..43cb309a6fcf929a4768ed99c0b3593a33723702 100644
--- a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
+++ b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
@@ -24,9 +24,11 @@
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 from scipy import stats
 import h5py
+import sys
 
 # Parameters
 gas_gamma = 5.0 / 3.0  # Polytropic index
@@ -37,25 +39,7 @@ v_R = 0.0  # Velocity right state
 P_L = 1.0  # Pressure left state
 P_R = 0.1  # Pressure right state
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -70,20 +54,20 @@ eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
 coords = sim["/PartType0/Coordinates"]
-x = sqrt(
+x = np.sqrt(
     (coords[:, 0] - 0.5 * boxSize) ** 2
     + (coords[:, 1] - 0.5 * boxSize) ** 2
     + (coords[:, 2] - 0.5 * boxSize) ** 2
 )
 vels = sim["/PartType0/Velocities"]
-v = sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2 + vels[:, 2] ** 2)
+v = np.sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2 + vels[:, 2] ** 2)
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
 # Bin the data
-x_bin_edge = logspace(-3.0, log10(2.0), 100)
+x_bin_edge = np.logspace(-3.0, np.log10(2.0), 100)
 x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
 rho_bin, _, _ = stats.binned_statistic(x, rho, statistic="mean", bins=x_bin_edge)
 v_bin, _, _ = stats.binned_statistic(x, v, statistic="mean", bins=x_bin_edge)
@@ -101,79 +85,117 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin ** 2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin ** 2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin ** 2)
 
-ref = loadtxt("evrardCollapse3D_exact.txt")
+ref = np.loadtxt("evrardCollapse3D_exact.txt")
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-semilogx(x, -v, ".", color="r", ms=0.2)
-semilogx(ref[:, 0], ref[:, 2], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, -v_bin, yerr=v_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
-xlim(1.0e-3, 2.0)
-ylim(-1.7, 0.1)
+plt.subplot(231)
+plt.semilogx(x, -v, **scatter_props)
+plt.semilogx(ref[:, 0], ref[:, 2], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, -v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(-1.7, 0.1)
 
 # Density profile --------------------------------
-subplot(232)
-loglog(x, rho, ".", color="r", ms=0.2)
-loglog(ref[:, 0], ref[:, 1], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(1.0e-3, 2.0)
-ylim(1.0e-2, 1.0e4)
+plt.subplot(232)
+plt.loglog(x, rho, **scatter_props)
+plt.loglog(ref[:, 0], ref[:, 1], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-2, 1.0e4)
 
 # Pressure profile --------------------------------
-subplot(233)
-loglog(x, P, ".", color="r", ms=0.2)
-loglog(ref[:, 0], ref[:, 3], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(1.0e-3, 2.0)
-ylim(1.0e-4, 1.0e3)
+plt.subplot(233)
+plt.loglog(x, P, **scatter_props)
+plt.loglog(ref[:, 0], ref[:, 3], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-4, 1.0e3)
 
 # Internal energy profile -------------------------
-subplot(234)
-loglog(x, u, ".", color="r", ms=0.2)
-loglog(ref[:, 0], ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0), "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(1.0e-3, 2.0)
-ylim(1.0e-2, 2.0)
+plt.subplot(234)
+plt.loglog(x, u, **scatter_props)
+plt.loglog(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0),
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-2, 2.0)
 
 # Entropy profile ---------------------------------
-subplot(235)
-semilogx(x, S, ".", color="r", ms=0.2)
-semilogx(ref[:, 0], ref[:, 3] / ref[:, 1] ** gas_gamma, "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-xlim(1.0e-3, 2.0)
-ylim(0.0, 0.25)
+plt.subplot(235)
+plt.semilogx(x, S, **scatter_props)
+plt.semilogx(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] ** gas_gamma,
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(0.0, 0.25)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(
-    -0.49,
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$Swift$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-tight_layout()
-savefig("EvrardCollapse.png", dpi=200)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("EvrardCollapse.png", dpi=200)
diff --git a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
index 5f1cea4f115358a7348e404a879210a2128518b9..0e143567930b856cb5af5e79ca6682df9974c435 100644
--- a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
+++ b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
@@ -17,126 +17,139 @@
 #
 ##############################################################################
 
-import numpy as np
 import matplotlib
 
 matplotlib.use("Agg")
-import pylab as pl
-import h5py
+import matplotlib.pyplot as plt
+import numpy as np
 import sys
+import h5py
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
-# Parameters
-gamma = 1.4  # Polytropic index
-
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-pl.rcParams.update(params)
-
-# Read the snapshot index from the command line argument
 snap = int(sys.argv[1])
 
 # Open the file and read the relevant data
 file = h5py.File("interactingBlastWaves_{0:04d}.hdf5".format(snap), "r")
 x = file["/PartType0/Coordinates"][:, 0]
-rho = file["/PartType0/Densities"]
+rho = file["/PartType0/Densities"][:]
 v = file["/PartType0/Velocities"][:, 0]
-u = file["/PartType0/InternalEnergies"]
-S = file["/PartType0/Entropies"]
-P = file["/PartType0/Pressures"]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
 kernel = file["/HydroScheme"].attrs["Kernel function"]
 neighbours = file["/HydroScheme"].attrs["Kernel target N_ngb"][0]
 eta = file["/HydroScheme"].attrs["Kernel eta"][0]
+gamma = file["/HydroScheme"].attrs["Adiabatic index"]
 git = file["Code"].attrs["Git Revision"]
 
+if gamma != 1.4:
+    print(
+        "Error: SWIFT was run with the wrong adiabatic index. Should have been 1.4",
+        gamma,
+    )
+    exit()
+
 ref = np.loadtxt("interactingBlastWaves1D_exact.txt")
 
 # Plot the interesting quantities
-fig, ax = pl.subplots(2, 3)
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=4,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
 
 # Velocity profile
-ax[0][0].plot(x, v, "r.", markersize=4.0)
-ax[0][0].plot(ref[:, 0], ref[:, 2], "k--", alpha=0.8, linewidth=1.2)
-ax[0][0].set_xlabel("${\\rm{Position}}~x$", labelpad=0)
-ax[0][0].set_ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
-ax[0][0].set_xlim(0.0, 1.0)
-ax[0][0].set_ylim(-1.0, 15.0)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+plt.xlim(0.0, 1.0)
+plt.ylim(-1.0, 15.0)
 
 # Density profile
-ax[0][1].plot(x, rho, "r.", markersize=4.0)
-ax[0][1].plot(ref[:, 0], ref[:, 1], "k--", alpha=0.8, linewidth=1.2)
-ax[0][1].set_xlabel("${\\rm{Position}}~x$", labelpad=0)
-ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-ax[0][1].set_xlim(0.0, 1.0)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Pressure profile
-ax[0][2].plot(x, P, "r.", markersize=4.0)
-ax[0][2].plot(ref[:, 0], ref[:, 3], "k--", alpha=0.8, linewidth=1.2)
-ax[0][2].set_xlabel("${\\rm{Position}}~x$", labelpad=0)
-ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-ax[0][2].set_xlim(0.0, 1.0)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Internal energy profile
-ax[1][0].plot(x, u, "r.", markersize=4.0)
-ax[1][0].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] / (gamma - 1.0), "k--", alpha=0.8, linewidth=1.2
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gamma - 1.0),
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
 )
-ax[1][0].set_xlabel("${\\rm{Position}}~x$", labelpad=0)
-ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-ax[1][0].set_xlim(0.0, 1.0)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Entropy profile
-ax[1][1].plot(x, S, "r.", markersize=4.0)
-ax[1][1].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, "k--", alpha=0.8, linewidth=1.2
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(
+    ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, "--", color=line_color, alpha=0.8, lw=1.2
 )
-ax[1][1].set_xlabel("${\\rm{Position}}~x$", labelpad=0)
-ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-ax[1][1].set_xlim(0.0, 1.0)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Run information
-ax[1][2].set_frame_on(False)
-ax[1][2].text(
-    -0.49,
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Interacting blast waves test\nwith $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
-        gamma, time
+        gamma[0], time
     ),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize=10)
-ax[1][2].text(-0.49, 0.4, scheme, fontsize=10)
-ax[1][2].text(-0.49, 0.3, kernel, fontsize=10)
-ax[1][2].text(
-    -0.49,
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
     0.2,
-    "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
-    fontsize=10,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-ax[1][2].set_xlim(-0.5, 0.5)
-ax[1][2].set_ylim(0.0, 1.0)
-ax[1][2].set_xticks([])
-ax[1][2].set_yticks([])
+plt.xlim(-0.5, 0.5)
+plt.ylim(0.0, 1.0)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
 
-pl.tight_layout()
-pl.savefig("InteractingBlastWaves.png", dpi=200)
+plt.savefig("InteractingBlastWaves.png", dpi=200)
diff --git a/examples/HydroTests/InteractingBlastWaves_1D/run.sh b/examples/HydroTests/InteractingBlastWaves_1D/run.sh
index 9139835bdbea88a80026f4a264a44bb60638882a..a3f927491ec07680756072beb378cc7903e374da 100755
--- a/examples/HydroTests/InteractingBlastWaves_1D/run.sh
+++ b/examples/HydroTests/InteractingBlastWaves_1D/run.sh
@@ -8,7 +8,7 @@ then
 fi
 
 # Run SWIFT
-../../../swift --hydro --threads=1 interactingBlastWaves.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=1  --limiter interactingBlastWaves.yml 2>&1 | tee output.log
 
 # Get the high resolution reference solution if not present.
 if [ ! -e interactingBlastWaves1D_exact.txt ]
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
index 41997c716e29bd6e005c0b2a9b88ba949f058fe7..ce5a394585c0b6a07e67d0365141346f123e2f34 100644
--- a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
+++ b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
@@ -35,30 +35,12 @@ rho2 = 1  # Outskirts density
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -83,11 +65,11 @@ S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 
 # Plot the interesting quantities
-figure()
+plt.figure()
 
 
 # Azimuthal velocity profile -----------------------------
-subplot(231)
+plt.subplot(231)
 scatter(
     pos[:, 0],
     pos[:, 1],
@@ -98,89 +80,88 @@ scatter(
     vmin=-1.0,
     vmax=1.0,
 )
-text(
+plt.text(
     0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w"
 )
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial density profile --------------------------------
-subplot(232)
-scatter(
+plt.subplot(232)
+plt.scatter(
     pos[:, 0], pos[:, 1], c=rho, cmap="PuBu", edgecolors="face", s=4, vmin=0.8, vmax=2.2
 )
-text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial pressure profile --------------------------------
-subplot(233)
-scatter(pos[:, 0], pos[:, 1], c=P, cmap="PuBu", edgecolors="face", s=4, vmin=1, vmax=4)
-text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.subplot(233)
+plt.scatter(
+    pos[:, 0], pos[:, 1], c=P, cmap="PuBu", edgecolors="face", s=4, vmin=1, vmax=4
+)
+plt.text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Internal energy profile --------------------------------
-subplot(234)
-scatter(
+plt.subplot(234)
+plt.scatter(
     pos[:, 0], pos[:, 1], c=u, cmap="PuBu", edgecolors="face", s=4, vmin=1.5, vmax=5.0
 )
-text(0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.text(
+    0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w"
+)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial entropy profile --------------------------------
-subplot(235)
-scatter(
+plt.subplot(235)
+plt.scatter(
     pos[:, 0], pos[:, 1], c=S, cmap="PuBu", edgecolors="face", s=4, vmin=0.5, vmax=3.0
 )
-text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
-
-# Image --------------------------------------------------
-# subplot(234)
-# scatter(pos[:,0], pos[:,1], c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=0, vmax=1)
-# text(0.95, 0.95, "$|v|$", ha="right", va="top")
-# xlim(0,1)
-# ylim(0,1)
-# xlabel("$x$", labelpad=0)
-# ylabel("$y$", labelpad=0)
+plt.text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Kelvin-Helmholtz instability at $t=%.2f$" % (time), fontsize=10)
-text(
-    -0.49,
+plt.text(-0.45, 0.9, "Kelvin-Helmholtz instability at $t=%.2f$" % (time), fontsize=10)
+plt.text(
+    -0.45,
     0.8,
-    "Centre:~~~ $(P, \\rho, v) = (%.3f, %.3f, %.3f)$" % (P1, rho1, v1),
+    "Centre: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$" % (P1, rho1, v1),
     fontsize=10,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.7,
     "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$" % (P2, rho2, v2),
     fontsize=10,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-savefig("KelvinHelmholtz.png", dpi=200)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=10)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=10)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=10)
+plt.text(
+    -0.45, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("KelvinHelmholtz.png", dpi=200)
diff --git a/examples/HydroTests/KeplerianRing/plotSolution.py b/examples/HydroTests/KeplerianRing/plotSolution.py
index 239bb639c9488f8f17b76c6ae8679338fd6bf7f9..25c3a846d972e900affda4b14aa2498fabad369c 100644
--- a/examples/HydroTests/KeplerianRing/plotSolution.py
+++ b/examples/HydroTests/KeplerianRing/plotSolution.py
@@ -306,7 +306,7 @@ def plot_extra_info(ax, filename):
     ax.text(-0.49, 0.8, f"Compiler: {compiler_name} {compiler_version}", fontsize=10)
     ax.text(-0.49, 0.7, "Rotations are quoted at $r=1$", fontsize=10)
     ax.plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-    ax.text(-0.49, 0.5, f"$\\textsc{{Swift}}$ {git}", fontsize=10)
+    ax.text(-0.49, 0.5, f"$SWIFT$ {git}", fontsize=10)
     ax.text(-0.49, 0.4, scheme, fontsize=10)
     ax.text(-0.49, 0.3, kernel, fontsize=10)
     ax.text(
diff --git a/examples/HydroTests/Noh_1D/plotSolution.py b/examples/HydroTests/Noh_1D/plotSolution.py
index 3d435c8f7ec3364137f53e0da5fbd44381c358a8..d4be4757a295fa086f00c39ad6e443c74dc63bc6 100644
--- a/examples/HydroTests/Noh_1D/plotSolution.py
+++ b/examples/HydroTests/Noh_1D/plotSolution.py
@@ -29,34 +29,15 @@ v0 = 1
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
 sim = h5py.File("noh_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
@@ -109,77 +90,101 @@ u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
 s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(x, v, ".", color="r", ms=4.0)
-plot(x_s, v_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_x$", labelpad=-4)
-xlim(-0.5, 0.5)
-ylim(-1.2, 1.2)
+plt.subplot(231)
+plt.plot(r, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_x$", labelpad=-4)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-1.2, 1.2)
 
 # Density profile --------------------------------
-subplot(232)
-plot(x, rho, ".", color="r", ms=4.0)
-plot(x_s, rho_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(0.95, 4.4)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0.95, 4.4)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(x, P, ".", color="r", ms=4.0)
-plot(x_s, P_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(-0.05, 1.8)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 1.8)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(x, u, ".", color="r", ms=4.0)
-plot(x_s, u_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(x, S, ".", color="r", ms=4.0)
-plot(x_s, s_s, "--", color="k", alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(-0.5, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(r, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(
-    -0.49,
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Noh problem with  $\\gamma=%.3f$ in 1D at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
-    "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
-    fontsize=10,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)
diff --git a/examples/HydroTests/Noh_2D/plotSolution.py b/examples/HydroTests/Noh_2D/plotSolution.py
index 1d6c62b26c83e73e56d2a5ff524de2cd51654d77..62e5437c1b82c4862137a0cb738a7eeb22726845 100644
--- a/examples/HydroTests/Noh_2D/plotSolution.py
+++ b/examples/HydroTests/Noh_2D/plotSolution.py
@@ -29,31 +29,13 @@ v0 = 1
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
 from scipy import stats
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -127,82 +109,106 @@ u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
 s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(r, v, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, v_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
-xlim(0, 0.5)
-ylim(-1.2, 0.4)
+plt.subplot(231)
+plt.plot(r, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
+plt.xlim(0, 0.5)
+plt.ylim(-1.2, 0.4)
 
 # Density profile --------------------------------
-subplot(232)
-plot(r, rho, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, rho_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(0, 0.5)
-ylim(0.95, 19)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(0.95, 19)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(r, P, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, P_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.5, 11)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.5, 11)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(r, u, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, u_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(r, S, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, s_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(0, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(r, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
 
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.9,
     "Noh problem with  $\\gamma=%.3f$ in 2D at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
-    "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
-    fontsize=10,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)
diff --git a/examples/HydroTests/Noh_3D/plotSolution.py b/examples/HydroTests/Noh_3D/plotSolution.py
index 31f6b520863b857b409bce2f2d7ebbb57aa05155..6b46ed07802b45abb4f4f0debea813f9b96ee4be 100644
--- a/examples/HydroTests/Noh_3D/plotSolution.py
+++ b/examples/HydroTests/Noh_3D/plotSolution.py
@@ -29,35 +29,16 @@ v0 = 1
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
 from scipy import stats
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
 sim = h5py.File("noh_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
@@ -133,82 +114,106 @@ u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
 s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(r, v, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, v_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
-xlim(0, 0.5)
-ylim(-1.2, 0.4)
+plt.subplot(231)
+plt.plot(r, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
+plt.xlim(0, 0.5)
+plt.ylim(-1.2, 0.4)
 
 # Density profile --------------------------------
-subplot(232)
-plot(r, rho, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, rho_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(0, 0.5)
-ylim(0.95, 71)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(0.95, 71)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(r, P, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, P_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.5, 25)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.5, 25)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(r, u, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, u_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(r, S, ".", color="r", ms=0.5, alpha=0.2)
-plot(x_s, s_s, "--", color="k", alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(0, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(r, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(
-    -0.49,
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Noh problem with  $\\gamma=%.3f$ in 3D at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
-    "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
-    fontsize=10,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)
diff --git a/examples/HydroTests/SedovBlast_1D/plotSolution.py b/examples/HydroTests/SedovBlast_1D/plotSolution.py
index 84976689d5e93effaa0d81e1d356848661e58580..93ecd501946b31d17e84376277813f64203abd0c 100644
--- a/examples/HydroTests/SedovBlast_1D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_1D/plotSolution.py
@@ -35,11 +35,13 @@ gas_gamma = 5.0 / 3.0  # Gas polytropic index
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 
 # Plot parameters
-style.use("../../../tools/stylesheets/mnras.mplstyle")
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -228,13 +230,13 @@ s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 
 # Plot the interesting quantities
-figure(figsize=(7, 7 / 1.6))
+plt.figure(figsize=(7, 7 / 1.6))
 
 line_color = "C4"
 
 scatter_props = dict(
     marker=".",
-    ms=1,
+    ms=4,
     markeredgecolor="none",
     alpha=1.0,
     zorder=-1,
@@ -243,88 +245,88 @@ scatter_props = dict(
 )
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(r, v_r, **scatter_props)
-plot(r_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
-xlabel("Radius $r$")
-ylabel("Radialvelocity $v_r$")
-xlim(0, 1.3 * r_shock)
-ylim(-0.2, 3.8)
+plt.subplot(231)
+plt.plot(r, v_r, **scatter_props)
+plt.plot(r_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("Radius $r$")
+plt.ylabel("Radialvelocity $v_r$")
+plt.xlim(0, 1.3 * r_shock)
+plt.ylim(-0.2, 3.8)
 
 # Density profile --------------------------------
-subplot(232)
-plot(r, rho, **scatter_props)
-plot(r_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
-xlabel("Radius $r$")
-ylabel("Density $\\rho$")
-xlim(0, 1.3 * r_shock)
-ylim(-0.2, 5.2)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(r_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("Radius $r$")
+plt.ylabel("Density $\\rho$")
+plt.xlim(0, 1.3 * r_shock)
+plt.ylim(-0.2, 5.2)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(r, P, **scatter_props)
-plot(r_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
-xlabel("Radius $r$")
-ylabel("Pressure $P$")
-xlim(0, 1.3 * r_shock)
-ylim(-1, 12.5)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(r_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("Radius $r$")
+plt.ylabel("Pressure $P$")
+plt.xlim(0, 1.3 * r_shock)
+plt.ylim(-1, 12.5)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(r, u, **scatter_props)
-plot(r_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
-xlabel("Radius $r$")
-ylabel("Internal Energy $u$")
-xlim(0, 1.3 * r_shock)
-ylim(-2, 22)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(r_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("Radius $r$")
+plt.ylabel("Internal Energy $u$")
+plt.xlim(0, 1.3 * r_shock)
+plt.ylim(-2, 22)
 
 # Entropy profile ---------------------------------
-subplot(235)
-xlabel("Radius $r$")
+plt.subplot(235)
+plt.xlabel("Radius $r$")
 if plot_diffusion or plot_viscosity:
     if plot_diffusion:
-        plot(r, diffusion, **scatter_props)
+        plt.plot(r, diffusion, **scatter_props)
 
     if plot_viscosity:
-        plot(r, viscosity, **scatter_props)
+        plt.plot(r, viscosity, **scatter_props)
 
-    ylabel(r"Rate Coefficient $\alpha$", labelpad=0)
-    legend()
+    plt.ylabel(r"Rate Coefficient $\alpha$", labelpad=0)
+    plt.legend()
 else:
-    plot(r, S, **scatter_props)
-    plot(r_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
-    ylabel("Entropy $S$", labelpad=0)
-    ylim(-5, 50)
+    plt.plot(r, S, **scatter_props)
+    plt.plot(r_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+    plt.ylabel("Entropy $S$", labelpad=0)
+    plt.ylim(-5, 50)
 
-xlim(0, 1.3 * r_shock)
+plt.xlim(0, 1.3 * r_shock)
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
 text_fontsize = 5
 
-text(
+plt.text(
     -0.45,
     0.9,
     "Sedov blast with  $\\gamma=%.3f$ in 3D at $t=%.2f$" % (gas_gamma, time),
     fontsize=text_fontsize,
 )
-text(-0.45, 0.8, "Background $\\rho_0=%.2f$" % (rho_0), fontsize=text_fontsize)
-text(-0.45, 0.7, "Energy injected $E_0=%.2f$" % (E_0), fontsize=text_fontsize)
-plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.45, 0.5, "SWIFT %s" % git.decode("utf-8"), fontsize=text_fontsize)
-text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
-text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
-text(
+plt.text(-0.45, 0.8, "Background $\\rho_0=%.2f$" % (rho_0), fontsize=text_fontsize)
+plt.text(-0.45, 0.7, "Energy injected $E_0=%.2f$" % (E_0), fontsize=text_fontsize)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "SWIFT %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
     -0.45,
     0.2,
     "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
     fontsize=text_fontsize,
 )
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
 
-tight_layout()
+plt.tight_layout()
 
-savefig("Sedov.png")
+plt.savefig("Sedov.png")
diff --git a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
index c1a490d0b8a453592abea9815ff3e2bb3aa5f4c9..c9b4f19cd7d5c80ab16ea1f160713f9831b88a6d 100644
--- a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
@@ -24,9 +24,11 @@
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
-from scipy import stats
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.stats as stats
 import h5py
+import sys
 
 # Parameters
 gas_gamma = 5.0 / 3.0  # Polytropic index
@@ -37,26 +39,7 @@ v_R = 0.0  # Velocity right state
 P_L = 1.0  # Pressure left state
 P_R = 0.1  # Pressure right state
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -71,9 +54,9 @@ eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
 coords = sim["/PartType0/Coordinates"]
-x = sqrt((coords[:, 0] - 0.5) ** 2 + (coords[:, 1] - 0.5) ** 2)
+x = np.sqrt((coords[:, 0] - 0.5) ** 2 + (coords[:, 1] - 0.5) ** 2)
 vels = sim["/PartType0/Velocities"]
-v = sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2)
+v = np.sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2)
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
@@ -97,81 +80,113 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin ** 2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin ** 2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin ** 2)
 
-ref = loadtxt("sodShockSpherical2D_exact.txt")
+ref = np.loadtxt("sodShockSpherical2D_exact.txt")
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.5,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(x, v, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 2], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, v_bin, yerr=v_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
 
 # Density profile --------------------------------
-subplot(232)
-plot(x, rho, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 1], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(x, P, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(x, u, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0), "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0),
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(x, S, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3] / ref[:, 1] ** gas_gamma, "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(
+    ref[:, 0], ref[:, 3] / ref[:, 1] ** gas_gamma, color=line_color, alpha=0.8, lw=1.2
+)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(
-    -0.49,
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Sod shock with  $\\gamma=%.3f$ in 2D at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
     "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$" % (P_L, rho_L, v_L),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.7,
     "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$" % (P_R, rho_R, v_R),
-    fontsize=10,
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("SodShock.png", dpi=200)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("SodShock.png", dpi=200)
diff --git a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
index b6eb53957d0f260897718608f619a84a797a8358..a548f8ed7cbacd5c4bae05164fadc18a6c2cf9d1 100644
--- a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
@@ -24,9 +24,11 @@
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
-from scipy import stats
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.stats as stats
 import h5py
+import sys
 
 # Parameters
 gas_gamma = 5.0 / 3.0  # Polytropic index
@@ -37,26 +39,7 @@ v_R = 0.0  # Velocity right state
 P_L = 1.0  # Pressure left state
 P_R = 0.1  # Pressure right state
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
@@ -71,11 +54,11 @@ eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
 coords = sim["/PartType0/Coordinates"]
-x = sqrt(
+x = np.sqrt(
     (coords[:, 0] - 0.5) ** 2 + (coords[:, 1] - 0.5) ** 2 + (coords[:, 2] - 0.5) ** 2
 )
 vels = sim["/PartType0/Velocities"]
-v = sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2 + vels[:, 2] ** 2)
+v = np.sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2 + vels[:, 2] ** 2)
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
@@ -99,81 +82,113 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin ** 2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin ** 2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin ** 2)
 
-ref = loadtxt("sodShockSpherical3D_exact.txt")
+ref = np.loadtxt("sodShockSpherical3D_exact.txt")
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.5,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(x, v, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 2], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, v_bin, yerr=v_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
 
 # Density profile --------------------------------
-subplot(232)
-plot(x, rho, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 1], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(x, P, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(x, u, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0), "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0),
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(x, S, ".", color="r", ms=0.2)
-plot(ref[:, 0], ref[:, 3] / ref[:, 1] ** gas_gamma, "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt=".", ms=8.0, color="b", lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(
+    ref[:, 0], ref[:, 3] / ref[:, 1] ** gas_gamma, color=line_color, alpha=0.8, lw=1.2
+)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
+plt.subplot(236, frameon=False)
 
-text(
-    -0.49,
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
     "Sod shock with  $\\gamma=%.3f$ in 3D at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
     "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$" % (P_L, rho_L, v_L),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-text(
-    -0.49,
+plt.text(
+    -0.45,
     0.7,
     "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$" % (P_R, rho_R, v_R),
-    fontsize=10,
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
 )
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("SodShock.png", dpi=200)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("SodShock.png", dpi=200)
diff --git a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
deleted file mode 100644
index 897aaebf0891cb6a1e2e821a1dafd79a4c22d9c8..0000000000000000000000000000000000000000
--- a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
+++ /dev/null
@@ -1,191 +0,0 @@
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016  Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-##############################################################################
-
-# Computes the analytical solution of the square test
-
-# Parameters
-gas_gamma = 5.0 / 3.0  # Gas adiabatic index
-gamma = 5.0 / 3.0  # Gas adiabatic index
-rho0 = 4  # Gas central density
-rho1 = 1  # Gas outskirt density
-P0 = 2.5  # Gas central pressure
-P1 = 2.5  # Gas central pressure
-vx = 0.0  # Random velocity for all particles
-vy = 0.0
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
-import matplotlib
-
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
-snap = int(sys.argv[1])
-
-# Read the simulation data
-sim = h5py.File("square_%04d.hdf5" % snap, "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-time = sim["/Header"].attrs["Time"][0]
-scheme = sim["/HydroScheme"].attrs["Scheme"]
-kernel = sim["/HydroScheme"].attrs["Kernel function"]
-neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
-eta = sim["/HydroScheme"].attrs["Kernel eta"]
-git = sim["Code"].attrs["Git Revision"]
-
-# Analytical soltion
-centre_x = 0.5 + time * vx
-centre_y = 0.5 + time * vy
-
-while centre_x > 1.0:
-    centre_x -= 1.0
-while centre_x < 0.0:
-    centre_x += 1.0
-while centre_y > 1.0:
-    centre_y -= 1.0
-while centre_y < 0.0:
-    centre_y += 1.0
-
-pos = sim["/PartType0/Coordinates"][:, :]
-vel = sim["/PartType0/Velocities"][:, :]
-v_norm = sqrt(vel[:, 0] ** 2 + vel[:, 1] ** 2)
-rho = sim["/PartType0/Densities"][:]
-u = sim["/PartType0/InternalEnergies"][:]
-S = sim["/PartType0/Entropies"][:]
-P = sim["/PartType0/Pressures"][:]
-x = pos[:, 0] - centre_x
-y = pos[:, 1] - centre_y
-
-# Box wrapping
-x[x > 0.5] -= 1.0
-x[x < -0.5] += 1.0
-y[y > 0.5] -= 1.0
-y[y < -0.5] += 1.0
-
-# Azimuthal velocity profile -----------------------------
-subplot(231)
-scatter(x, y, c=v_norm, cmap="PuBu", edgecolors="face", s=4, vmin=-1.0, vmax=1.0)
-text(0.47, 0.47, "${\\rm{Velocity~norm}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], "--", color="k", alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial density profile --------------------------------
-subplot(232)
-scatter(x, y, c=rho, cmap="PuBu", edgecolors="face", s=4, vmin=0.0, vmax=4.0)
-text(0.47, 0.47, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], "--", color="k", alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial pressure profile --------------------------------
-subplot(233)
-scatter(x, y, c=P, cmap="PuBu", edgecolors="face", s=4, vmin=2, vmax=4)
-text(0.47, 0.47, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], "--", color="k", alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Internal energy profile --------------------------------
-subplot(234)
-scatter(x, y, c=u, cmap="PuBu", edgecolors="face", s=4, vmin=0.5, vmax=4.0)
-text(0.47, 0.47, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], "--", color="k", alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial entropy profile --------------------------------
-subplot(235)
-scatter(x, y, c=S, cmap="PuBu", edgecolors="face", s=4, vmin=0.0, vmax=3.0)
-text(0.47, 0.47, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], "--", color="k", alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], "--", color="k", alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-
-# Information -------------------------------------
-subplot(236, frameon=False)
-
-text(
-    -0.49,
-    0.9,
-    "Square test with $\\gamma=%.3f$ at $t=%.2f$" % (gas_gamma, time),
-    fontsize=10,
-)
-text(-0.49, 0.8, "Centre:~~~ $(P, \\rho) = (%.3f, %.3f)$" % (P0, rho0), fontsize=10)
-text(-0.49, 0.7, "Outskirts: $(P, \\rho) = (%.3f, %.3f)$" % (P1, rho1), fontsize=10)
-plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s" % git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-savefig("SquareTest.png", dpi=200)
diff --git a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
index bf53bc8373338937db2fa4092dcece9e07f528b2..ef6f398ec7105bc946cb375167c93b8cc067c1bc 100644
--- a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
@@ -17,14 +17,14 @@
 #
 ##############################################################################
 
-import numpy as np
 import matplotlib
 
 matplotlib.use("Agg")
-import pylab as pl
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import stats
 import h5py
 import sys
-import scipy.stats as stats
 
 # Parameters
 gamma = 5.0 / 3.0  # Polytropic index
@@ -35,28 +35,8 @@ rhoR = 0.0  # Initial vacuum density
 vR = 0.0  # Initial vacuum velocity
 PR = 0.0  # Initial vacuum pressure
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-pl.rcParams.update(params)
-
-# Read the snapshot index from the command line argument
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+
 snap = int(sys.argv[1])
 
 # Open the file and read the relevant data
@@ -104,104 +84,115 @@ x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
 ref = np.loadtxt("vacuumSpherical2D_exact.txt")
 
 # Plot the interesting quantities
-fig, ax = pl.subplots(2, 3)
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.5,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile
-ax[0][0].plot(x, v, "r.", markersize=0.2)
-ax[0][0].plot(ref[:, 0], ref[:, 2], "k--", alpha=0.8, linewidth=1.2)
-ax[0][0].errorbar(
-    x_bin, v_bin, yerr=v_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[0][0].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][0].set_ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
-ax[0][0].set_xlim(0.0, 0.4)
-ax[0][0].set_ylim(-0.1, 3.2)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.xlim(0.0, 0.4)
+plt.ylim(-0.1, 3.2)
 
 # Density profile
-ax[0][1].plot(x, rho, "r.", markersize=0.2)
-ax[0][1].plot(ref[:, 0], ref[:, 1], "k--", alpha=0.8, linewidth=1.2)
-ax[0][1].errorbar(
-    x_bin,
-    rho_bin,
-    yerr=rho_sigma_bin,
-    fmt=".",
-    markersize=8.0,
-    color="b",
-    linewidth=1.2,
-)
-ax[0][1].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-ax[0][1].set_xlim(0.0, 0.4)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Pressure profile
-ax[0][2].plot(x, P, "r.", markersize=0.2)
-ax[0][2].plot(ref[:, 0], ref[:, 3], "k--", alpha=0.8, linewidth=1.2)
-ax[0][2].errorbar(
-    x_bin, P_bin, yerr=P_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[0][2].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-ax[0][2].set_xlim(0.0, 0.4)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Internal energy profile
-ax[1][0].plot(x, u, "r.", markersize=0.2)
-ax[1][0].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] / (gamma - 1.0), "k--", alpha=0.8, linewidth=1.2
-)
-ax[1][0].errorbar(
-    x_bin, u_bin, yerr=u_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gamma - 1.0),
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
 )
-ax[1][0].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-ax[1][0].set_xlim(0.0, 0.4)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Entropy profile
-ax[1][1].plot(x, S, "r.", markersize=0.2)
-ax[1][1].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, "k--", alpha=0.8, linewidth=1.2
-)
-ax[1][1].errorbar(
-    x_bin, S_bin, yerr=S_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[1][1].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-ax[1][1].set_xlim(0.0, 0.4)
-ax[1][1].set_ylim(0.0, 4.0)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(0.0, 0.4)
+plt.ylim(0.0, 4.0)
 
 # Run information
-ax[1][2].set_frame_on(False)
-ax[1][2].text(
-    -0.49,
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
-    "Vacuum test with $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(gamma, time),
-    fontsize=10,
+    "Vacuum test with $\\gamma={0:.3f}$ in 2D at $t = {1:.2f}$".format(gamma, time),
+    fontsize=text_fontsize,
 )
-ax[1][2].text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
-    "Left:~~ $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PL, rhoL, vL),
-    fontsize=10,
+    "Left: $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PL, rhoL, vL),
+    fontsize=text_fontsize,
 )
-ax[1][2].text(
-    -0.49,
+plt.text(
+    -0.45,
     0.7,
     "Right: $(P_R, \\rho_R, v_R) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PR, rhoR, vR),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize=10)
-ax[1][2].text(-0.49, 0.4, scheme, fontsize=10)
-ax[1][2].text(-0.49, 0.3, kernel, fontsize=10)
-ax[1][2].text(
-    -0.49,
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ {0}".format(git.decode("utf-8")), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
     0.2,
     "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-ax[1][2].set_xlim(-0.5, 0.5)
-ax[1][2].set_ylim(0.0, 1.0)
-ax[1][2].set_xticks([])
-ax[1][2].set_yticks([])
+plt.xlim(-0.5, 0.5)
+plt.ylim(0.0, 1.0)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
 
-pl.tight_layout()
-pl.savefig("Vacuum.png", dpi=200)
+plt.savefig("Vacuum.png", dpi=200)
diff --git a/examples/HydroTests/VacuumSpherical_3D/getGlass.sh b/examples/HydroTests/VacuumSpherical_3D/getGlass.sh
old mode 100644
new mode 100755
diff --git a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
index dc0b491cbc561163796477ac62a460e6805f2442..cb5ab373718cfcd4c3c99717a59d48f8a23501f0 100644
--- a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
@@ -17,14 +17,14 @@
 #
 ##############################################################################
 
-import numpy as np
 import matplotlib
 
 matplotlib.use("Agg")
-import pylab as pl
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import stats
 import h5py
 import sys
-import scipy.stats as stats
 
 # Parameters
 gamma = 5.0 / 3.0  # Polytropic index
@@ -35,28 +35,8 @@ rhoR = 0.0  # Initial vacuum density
 vR = 0.0  # Initial vacuum velocity
 PR = 0.0  # Initial vacuum pressure
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.045,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.05,
-    "figure.subplot.top": 0.99,
-    "figure.subplot.wspace": 0.15,
-    "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-pl.rcParams.update(params)
-
-# Read the snapshot index from the command line argument
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+
 snap = int(sys.argv[1])
 
 # Open the file and read the relevant data
@@ -105,105 +85,115 @@ x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
 
 ref = np.loadtxt("vacuumSpherical3D_exact.txt")
 
-# Plot the interesting quantities
-fig, ax = pl.subplots(2, 3)
+plt.figure(figsize=(7, 7 / 1.6))
 
-# Velocity profile
-ax[0][0].plot(x, v, "r.", markersize=0.2)
-ax[0][0].plot(ref[:, 0], ref[:, 2], "k--", alpha=0.8, linewidth=1.2)
-ax[0][0].errorbar(
-    x_bin, v_bin, yerr=v_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.5,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
 )
-ax[0][0].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][0].set_ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
-ax[0][0].set_xlim(0.0, 0.4)
-ax[0][0].set_ylim(-0.1, 3.2)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
+
+# Velocity profile
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.xlim(0.0, 0.4)
+plt.ylim(-0.1, 3.2)
 
 # Density profile
-ax[0][1].plot(x, rho, "r.", markersize=0.2)
-ax[0][1].plot(ref[:, 0], ref[:, 1], "k--", alpha=0.8, linewidth=1.2)
-ax[0][1].errorbar(
-    x_bin,
-    rho_bin,
-    yerr=rho_sigma_bin,
-    fmt=".",
-    markersize=8.0,
-    color="b",
-    linewidth=1.2,
-)
-ax[0][1].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-ax[0][1].set_xlim(0.0, 0.4)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Pressure profile
-ax[0][2].plot(x, P, "r.", markersize=0.2)
-ax[0][2].plot(ref[:, 0], ref[:, 3], "k--", alpha=0.8, linewidth=1.2)
-ax[0][2].errorbar(
-    x_bin, P_bin, yerr=P_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[0][2].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-ax[0][2].set_xlim(0.0, 0.4)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Internal energy profile
-ax[1][0].plot(x, u, "r.", markersize=0.2)
-ax[1][0].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] / (gamma - 1.0), "k--", alpha=0.8, linewidth=1.2
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gamma - 1.0),
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
 )
-ax[1][0].errorbar(
-    x_bin, u_bin, yerr=u_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[1][0].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-ax[1][0].set_xlim(0.0, 0.4)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0.0, 0.4)
 
 # Entropy profile
-ax[1][1].plot(x, S, "r.", markersize=0.2)
-ax[1][1].plot(
-    ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, "k--", alpha=0.8, linewidth=1.2
-)
-ax[1][1].errorbar(
-    x_bin, S_bin, yerr=S_sigma_bin, fmt=".", markersize=8.0, color="b", linewidth=1.2
-)
-ax[1][1].set_xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-ax[1][1].set_xlim(0.0, 0.4)
-ax[1][1].set_ylim(0.0, 4.0)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(0.0, 0.4)
+plt.ylim(0.0, 4.0)
 
 # Run information
-ax[1][2].set_frame_on(False)
-ax[1][2].text(
-    -0.49,
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
     0.9,
-    "Vacuum test with $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(gamma, time),
-    fontsize=10,
+    "Vacuum test with $\\gamma={0:.3f}$ in 3D at $t = {1:.2f}$".format(gamma, time),
+    fontsize=text_fontsize,
 )
-ax[1][2].text(
-    -0.49,
+plt.text(
+    -0.45,
     0.8,
-    "Left:~~ $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PL, rhoL, vL),
-    fontsize=10,
+    "Left: $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PL, rhoL, vL),
+    fontsize=text_fontsize,
 )
-ax[1][2].text(
-    -0.49,
+plt.text(
+    -0.45,
     0.7,
     "Right: $(P_R, \\rho_R, v_R) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(PR, rhoR, vR),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
-ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize=10)
-ax[1][2].text(-0.49, 0.4, scheme, fontsize=10)
-ax[1][2].text(-0.49, 0.3, kernel, fontsize=10)
-ax[1][2].text(
-    -0.49,
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ {0}".format(git.decode("utf-8")), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
     0.2,
     "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
-    fontsize=10,
+    fontsize=text_fontsize,
 )
-ax[1][2].set_xlim(-0.5, 0.5)
-ax[1][2].set_ylim(0.0, 1.0)
-ax[1][2].set_xticks([])
-ax[1][2].set_yticks([])
+plt.xlim(-0.5, 0.5)
+plt.ylim(0.0, 1.0)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
 
-pl.tight_layout()
-pl.savefig("Vacuum.png", dpi=200)
+plt.savefig("Vacuum.png", dpi=200)
diff --git a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
index 598ebb41adddf23e474df252a739e9801cf62c04..68b07ca6ec7390b3eab4eb90b5cd991949d03326 100644
--- a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
+++ b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
@@ -30,30 +30,11 @@ snapname = "santabarbara"
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 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.0,
-}
-rcParams.update(params)
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Read the simulation data
 sim = h5py.File("%s_0000.hdf5" % snapname, "r")
@@ -63,6 +44,7 @@ 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"
@@ -70,6 +52,10 @@ H_transition_temp = sim["/HydroScheme"].attrs[
 T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
 T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
 git = sim["Code"].attrs["Git Revision"]
+cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"].decode("utf-8")
+
+if cooling_model == "Constant Lambda":
+    Lambda = sim["/SubgridScheme"].attrs["Lambda/n_H^2 [cgs]"][0]
 
 # Cosmological parameters
 H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
@@ -83,7 +69,7 @@ unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
 
-# Primoridal ean molecular weight as a function of temperature
+# Primoridal mean 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.0 / (8.0 - 5.0 * (1.0 - H_frac))
@@ -147,28 +133,46 @@ a_evol = np.logspace(-3, 0, 60)
 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.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(
+plt.figure()
+plt.subplot(111, xscale="log", yscale="log")
+
+plt.fill_between(a, T_mean - T_std, T_mean + T_std, color="C0", alpha=0.1)
+plt.plot(a, T_max, ls="-.", color="C0", lw=1.0, label="${\\rm max}~T$")
+plt.plot(a, T_min, ls=":", color="C0", lw=1.0, label="${\\rm min}~T$")
+plt.plot(a, T_mean, color="C0", label="${\\rm mean}~T$", lw=1.5)
+plt.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)
+plt.plot(a, 10 ** T_log_mean, color="C1", label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plt.plot(a, T_median, color="C2", label="${\\rm median}~T$", lw=1.5)
+
+plt.legend(loc="upper left", frameon=False, handlelength=1.5)
+
+# Cooling model
+if cooling_model == "Constant Lambda":
+    plt.text(
+        1e-2,
+        6e4,
+        "$\Lambda_{\\rm const}/n_{\\rm H}^2 = %.1f\\times10^{%d}~[\\rm{cgs}]$"
+        % (Lambda / 10.0 ** (int(log10(Lambda))), log10(Lambda)),
+        fontsize=7,
+    )
+elif cooling_model == "EAGLE":
+    plt.text(1e-2, 6e4, "EAGLE (Wiersma et al. 2009)")
+elif cooling_model == b"Grackle":
+    plt.text(1e-2, 6e4, "Grackle (Smith et al. 2016)")
+else:
+    plt.text(1e-2, 6e4, "No cooling")
 
 # Expected lines
-plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7)
-text(
+plt.plot(
+    [1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7
+)
+plt.text(
     2.5e-2,
     H_transition_temp * 1.07,
     "$T_{\\rm HII\\rightarrow HI}$",
@@ -176,10 +180,10 @@ text(
     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(
+plt.plot([1e-10, 1e10], [T_minimal, T_minimal], "k--", lw=0.5, alpha=0.7)
+plt.text(1e-2, T_minimal * 0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plt.plot(a_evol, T_cmb, "k--", lw=0.5, alpha=0.7)
+plt.text(
     a_evol[20],
     T_cmb[20] * 0.55,
     "$(1+z)^2\\times T_{\\rm CMB,0}$",
@@ -195,12 +199,14 @@ 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.0 / (redshift_ticks + 1.0)
 
-xticks(a_ticks, redshift_labels)
-minorticks_off()
+plt.xticks(a_ticks, redshift_labels)
+plt.minorticks_off()
+
+plt.xlabel("${\\rm Redshift}~z$", labelpad=0)
+plt.ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+plt.xlim(9e-3, 1.1)
+plt.ylim(20, 2.5e7)
 
-xlabel("${\\rm Redshift}~z$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
-xlim(9e-3, 1.1)
-ylim(20, 2.5e7)
+plt.tight_layout()
 
-savefig("Temperature_evolution.png", dpi=200)
+plt.savefig("Temperature_evolution.png", dpi=200)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py
deleted file mode 100644
index 3eb80e98e8dfef70cf5dedaad876f37a54a587f5..0000000000000000000000000000000000000000
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py
+++ /dev/null
@@ -1,206 +0,0 @@
-################################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-################################################################################
-
-# Computes the temperature evolution of the gas in a cosmological box
-
-# Physical constants needed for internal energy to temperature conversion
-k_in_J_K = 1.38064852e-23
-mH_in_kg = 1.6737236e-27
-
-# Number of snapshots generated
-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.0,
-}
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("snap_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-time = sim["/Header"].attrs["Time"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-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.0 / (8.0 - 5.0 * (1.0 - H_frac))
-    else:
-        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.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)
-
-    # Neutral gas?
-    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
-
-
-z = np.zeros(n_snapshots)
-a = np.zeros(n_snapshots)
-T_mean = np.zeros(n_snapshots)
-T_std = np.zeros(n_snapshots)
-T_log_mean = np.zeros(n_snapshots)
-T_log_std = np.zeros(n_snapshots)
-T_median = np.zeros(n_snapshots)
-T_min = np.zeros(n_snapshots)
-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")
-
-    z[i] = sim["/Cosmology"].attrs["Redshift"][0]
-    a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
-
-    u = sim["/PartType0/InternalEnergies"][:]
-
-    # Compute the temperature
-    u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
-    u /= a[i] ** (3 * (gas_gamma - 1.0))
-    Temp = T(u)
-
-    # Gather statistics
-    T_median[i] = np.median(Temp)
-    T_mean[i] = Temp.mean()
-    T_std[i] = Temp.std()
-    T_log_mean[i] = np.log10(Temp).mean()
-    T_log_std[i] = np.log10(Temp).std()
-    T_min[i] = Temp.min()
-    T_max[i] = Temp.max()
-
-# CMB evolution
-a_evol = np.logspace(-3, 0, 60)
-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.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.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.0 / (redshift_ticks + 1.0)
-
-xticks(a_ticks, redshift_labels)
-minorticks_off()
-
-xlabel("${\\rm Redshift}~z$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
-xlim(9e-3, 1.1)
-ylim(20, 2.5e7)
-
-savefig("Temperature_evolution.png", dpi=200)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
index e5db3348cd0f8c389643cbc12e79e58ccbaf8cbe..d758fd8a9f9ee81050e4336cb5a6421563532f92 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
@@ -26,35 +26,17 @@ mH_in_kg = 1.6737236e-27
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 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.0,
-}
-rcParams.update(params)
+import sys
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
 # Read the simulation data
-sim = h5py.File("snap_%04d.hdf5" % snap, "r")
+sim = h5py.File("snapshots/snap_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 z = sim["/Cosmology"].attrs["Redshift"][0]
@@ -84,7 +66,7 @@ unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
 
-# Primoridal ean molecular weight as a function of temperature
+# Primoridal mean 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.0 / (8.0 - 5.0 * (1.0 - H_frac))
@@ -136,26 +118,28 @@ 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, _, _ = np.histogram2d(log_rho, log_T, bins=[bins_x, bins_y], normed=True)
 
 
 # Plot the interesting quantities
-figure()
+plt.figure()
 
-pcolormesh(bins_x, bins_y, np.log10(H).T)
+plt.pcolormesh(bins_x, bins_y, np.log10(H).T)
 
-text(-5, 8.0, "$z=%.2f$" % z)
+plt.text(-5, 8.0, "$z=%.2f$" % z)
 
-xticks(
+plt.xticks(
     [-5, -4, -3, -2, -1, 0, 1, 2, 3],
     ["", "$10^{-4}$", "", "$10^{-2}$", "", "$10^0$", "", "$10^2$", ""],
 )
-yticks(
+plt.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)
+plt.xlabel("${\\rm Physical~Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+plt.ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+plt.xlim(-5.2, 3.2)
+plt.ylim(1, 8.5)
+
+plt.tight_layout()
 
-savefig("rhoT_%04d.png" % snap, dpi=200)
+plt.savefig("rhoT_%04d.png" % snap, dpi=200)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
index 006f7742f23a8a0400b40ae9618e1d9d70868301..f42dd927d45384e65f7ee2fb8b41d93a3df7b1cc 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -29,33 +29,14 @@ n_snapshots = 200
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 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.0,
-}
-rcParams.update(params)
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Read the simulation data
-sim = h5py.File("snap_0000.hdf5", "r")
+sim = h5py.File("snapshots/snap_0000.hdf5", "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"][0]
@@ -70,7 +51,7 @@ H_transition_temp = sim["/HydroScheme"].attrs[
 T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
 T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
 git = sim["Code"].attrs["Git Revision"]
-cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"]
+cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"].decode("utf-8")
 
 if cooling_model == "Constant Lambda":
     Lambda = sim["/SubgridScheme"].attrs["Lambda/n_H^2 [cgs]"][0]
@@ -125,7 +106,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("snapshots/snap_%04d.hdf5" % i, "r")
 
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
@@ -151,28 +132,28 @@ a_evol = np.logspace(-3, 0, 60)
 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.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(
+plt.figure()
+plt.subplot(111, xscale="log", yscale="log")
+
+plt.fill_between(a, T_mean - T_std, T_mean + T_std, color="C0", alpha=0.1)
+plt.plot(a, T_max, ls="-.", color="C0", lw=1.0, label="${\\rm max}~T$")
+plt.plot(a, T_min, ls=":", color="C0", lw=1.0, label="${\\rm min}~T$")
+plt.plot(a, T_mean, color="C0", label="${\\rm mean}~T$", lw=1.5)
+plt.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)
+plt.plot(a, 10 ** T_log_mean, color="C1", label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plt.plot(a, T_median, color="C2", label="${\\rm median}~T$", lw=1.5)
 
-legend(loc="upper left", frameon=False, handlelength=1.5)
+plt.legend(loc="upper left", frameon=False, handlelength=1.5)
 
 # Cooling model
 if cooling_model == "Constant Lambda":
-    text(
+    plt.text(
         1e-2,
         6e4,
         "$\Lambda_{\\rm const}/n_{\\rm H}^2 = %.1f\\times10^{%d}~[\\rm{cgs}]$"
@@ -180,15 +161,17 @@ if cooling_model == "Constant Lambda":
         fontsize=7,
     )
 elif cooling_model == "EAGLE":
-    text(1e-2, 6e4, "EAGLE (Wiersma et al. 2009)")
+    plt.text(1e-2, 6e4, "EAGLE (Wiersma et al. 2009)")
 elif cooling_model == b"Grackle":
-    text(1e-2, 6e4, "Grackle (Smith et al. 2016)")
+    plt.text(1e-2, 6e4, "Grackle (Smith et al. 2016)")
 else:
-    text(1e-2, 6e4, "No cooling")
+    plt.text(1e-2, 6e4, "No cooling")
 
 # Expected lines
-plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7)
-text(
+plt.plot(
+    [1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7
+)
+plt.text(
     2.5e-2,
     H_transition_temp * 1.07,
     "$T_{\\rm HII\\rightarrow HI}$",
@@ -196,10 +179,10 @@ text(
     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(
+plt.plot([1e-10, 1e10], [T_minimal, T_minimal], "k--", lw=0.5, alpha=0.7)
+plt.text(1e-2, T_minimal * 0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plt.plot(a_evol, T_cmb, "k--", lw=0.5, alpha=0.7)
+plt.text(
     a_evol[20],
     T_cmb[20] * 0.55,
     "$(1+z)^2\\times T_{\\rm CMB,0}$",
@@ -215,12 +198,14 @@ 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.0 / (redshift_ticks + 1.0)
 
-xticks(a_ticks, redshift_labels)
-minorticks_off()
+plt.xticks(a_ticks, redshift_labels)
+plt.minorticks_off()
+
+plt.xlabel("${\\rm Redshift}~z$", labelpad=0)
+plt.ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+plt.xlim(9e-3, 1.1)
+plt.ylim(20, 2.5e7)
 
-xlabel("${\\rm Redshift}~z$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
-xlim(9e-3, 1.1)
-ylim(5, 2.5e7)
+plt.tight_layout()
 
-savefig("Temperature_evolution.png", dpi=200)
+plt.savefig("Temperature_evolution.png", dpi=200)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index 1f93830629811b557705427204520b8c8db9ce53..c308cbd5382db94d36fef8e6d0d9b89594d08099 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -45,6 +45,7 @@ Snapshots:
   basename:            snap
   delta_time:          1.02
   scale_factor_first:  0.02
+  compression:         4
   
 # Parameters governing the conserved quantities statistics
 Statistics:
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
index 3eb80e98e8dfef70cf5dedaad876f37a54a587f5..f42dd927d45384e65f7ee2fb8b41d93a3df7b1cc 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
@@ -29,33 +29,14 @@ n_snapshots = 200
 import matplotlib
 
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 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.0,
-}
-rcParams.update(params)
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Read the simulation data
-sim = h5py.File("snap_0000.hdf5", "r")
+sim = h5py.File("snapshots/snap_0000.hdf5", "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"][0]
@@ -70,6 +51,10 @@ H_transition_temp = sim["/HydroScheme"].attrs[
 T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
 T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
 git = sim["Code"].attrs["Git Revision"]
+cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"].decode("utf-8")
+
+if cooling_model == "Constant Lambda":
+    Lambda = sim["/SubgridScheme"].attrs["Lambda/n_H^2 [cgs]"][0]
 
 # Cosmological parameters
 H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
@@ -83,7 +68,7 @@ unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
 
-# Primoridal ean molecular weight as a function of temperature
+# Primoridal mean 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.0 / (8.0 - 5.0 * (1.0 - H_frac))
@@ -121,7 +106,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("snapshots/snap_%04d.hdf5" % i, "r")
 
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
@@ -147,28 +132,46 @@ a_evol = np.logspace(-3, 0, 60)
 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.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(
+plt.figure()
+plt.subplot(111, xscale="log", yscale="log")
+
+plt.fill_between(a, T_mean - T_std, T_mean + T_std, color="C0", alpha=0.1)
+plt.plot(a, T_max, ls="-.", color="C0", lw=1.0, label="${\\rm max}~T$")
+plt.plot(a, T_min, ls=":", color="C0", lw=1.0, label="${\\rm min}~T$")
+plt.plot(a, T_mean, color="C0", label="${\\rm mean}~T$", lw=1.5)
+plt.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)
+plt.plot(a, 10 ** T_log_mean, color="C1", label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plt.plot(a, T_median, color="C2", label="${\\rm median}~T$", lw=1.5)
+
+plt.legend(loc="upper left", frameon=False, handlelength=1.5)
+
+# Cooling model
+if cooling_model == "Constant Lambda":
+    plt.text(
+        1e-2,
+        6e4,
+        "$\Lambda_{\\rm const}/n_{\\rm H}^2 = %.1f\\times10^{%d}~[\\rm{cgs}]$"
+        % (Lambda / 10.0 ** (int(log10(Lambda))), log10(Lambda)),
+        fontsize=7,
+    )
+elif cooling_model == "EAGLE":
+    plt.text(1e-2, 6e4, "EAGLE (Wiersma et al. 2009)")
+elif cooling_model == b"Grackle":
+    plt.text(1e-2, 6e4, "Grackle (Smith et al. 2016)")
+else:
+    plt.text(1e-2, 6e4, "No cooling")
 
 # Expected lines
-plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7)
-text(
+plt.plot(
+    [1e-10, 1e10], [H_transition_temp, H_transition_temp], "k--", lw=0.5, alpha=0.7
+)
+plt.text(
     2.5e-2,
     H_transition_temp * 1.07,
     "$T_{\\rm HII\\rightarrow HI}$",
@@ -176,10 +179,10 @@ text(
     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(
+plt.plot([1e-10, 1e10], [T_minimal, T_minimal], "k--", lw=0.5, alpha=0.7)
+plt.text(1e-2, T_minimal * 0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plt.plot(a_evol, T_cmb, "k--", lw=0.5, alpha=0.7)
+plt.text(
     a_evol[20],
     T_cmb[20] * 0.55,
     "$(1+z)^2\\times T_{\\rm CMB,0}$",
@@ -195,12 +198,14 @@ 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.0 / (redshift_ticks + 1.0)
 
-xticks(a_ticks, redshift_labels)
-minorticks_off()
+plt.xticks(a_ticks, redshift_labels)
+plt.minorticks_off()
+
+plt.xlabel("${\\rm Redshift}~z$", labelpad=0)
+plt.ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+plt.xlim(9e-3, 1.1)
+plt.ylim(20, 2.5e7)
 
-xlabel("${\\rm Redshift}~z$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
-xlim(9e-3, 1.1)
-ylim(20, 2.5e7)
+plt.tight_layout()
 
-savefig("Temperature_evolution.png", dpi=200)
+plt.savefig("Temperature_evolution.png", dpi=200)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml
index da0a2c7686c9b2957043c77cd9f9cadb6aa7cc09..518800b6150593caa40d580f1dd85d32feee761d 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml
@@ -41,9 +41,10 @@ SPH:
 
 # Parameters governing the snapshots
 Snapshots:
+  subdir:              snapshots
   basename:            snap
-  delta_time:          1.05
-  scale_factor_first:  0.05
+  delta_time:          1.02
+  scale_factor_first:  0.02
   compression:         4
   
 # Parameters governing the conserved quantities statistics
@@ -61,8 +62,8 @@ InitialConditions:
   periodic:                    1
   cleanup_h_factors:           1    
   cleanup_velocity_factors:    1  
-  generate_gas_in_ics:         1      # Generate gas particles from the DM-only ICs
-  cleanup_smoothing_lengths:   1      # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  generate_gas_in_ics:         1    # Generate gas particles from the DM-only ICs
+  cleanup_smoothing_lengths:   1    # Since we generate gas, make use of the (expensive) cleaning-up procedure.
 
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotRhoT.py b/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotRhoT.py
deleted file mode 100644
index e5db3348cd0f8c389643cbc12e79e58ccbaf8cbe..0000000000000000000000000000000000000000
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotRhoT.py
+++ /dev/null
@@ -1,161 +0,0 @@
-################################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-################################################################################
-
-# Computes the temperature evolution of the gas in a cosmological box
-
-# Physical constants needed for internal energy to temperature conversion
-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.0,
-}
-rcParams.update(params)
-
-snap = int(sys.argv[1])
-
-# Read the simulation data
-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]
-a = sim["/Cosmology"].attrs["Scale-factor"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-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.0 / (8.0 - 5.0 * (1.0 - H_frac))
-    else:
-        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.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)
-
-    # Neutral gas?
-    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
-
-
-rho = sim["/PartType0/Densities"][:]
-u = sim["/PartType0/InternalEnergies"][:]
-
-# Compute the temperature
-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 /= mH_in_kg
-
-# Life is better in log-space
-log_T = np.log10(Temp)
-log_rho = np.log10(rho)
-
-
-# Make a 2D histogram
-log_rho_min = -6
-log_rho_max = 3
-log_T_min = 1
-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)
-
-
-# Plot the interesting quantities
-figure()
-
-pcolormesh(bins_x, bins_y, np.log10(H).T)
-
-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$"]
-)
-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)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotTempEvolution.py
deleted file mode 100644
index 006f7742f23a8a0400b40ae9618e1d9d70868301..0000000000000000000000000000000000000000
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/plotTempEvolution.py
+++ /dev/null
@@ -1,226 +0,0 @@
-################################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-################################################################################
-
-# Computes the temperature evolution of the gas in a cosmological box
-
-# Physical constants needed for internal energy to temperature conversion
-k_in_J_K = 1.38064852e-23
-mH_in_kg = 1.6737236e-27
-
-# Number of snapshots generated
-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.0,
-}
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("snap_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-time = sim["/Header"].attrs["Time"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"]
-
-if cooling_model == "Constant Lambda":
-    Lambda = sim["/SubgridScheme"].attrs["Lambda/n_H^2 [cgs]"][0]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-unit_time_in_si = unit_time_in_cgs
-
-# Primoridal mean 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.0 / (8.0 - 5.0 * (1.0 - H_frac))
-    else:
-        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.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)
-
-    # Neutral gas?
-    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
-
-
-z = np.zeros(n_snapshots)
-a = np.zeros(n_snapshots)
-T_mean = np.zeros(n_snapshots)
-T_std = np.zeros(n_snapshots)
-T_log_mean = np.zeros(n_snapshots)
-T_log_std = np.zeros(n_snapshots)
-T_median = np.zeros(n_snapshots)
-T_min = np.zeros(n_snapshots)
-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")
-
-    z[i] = sim["/Cosmology"].attrs["Redshift"][0]
-    a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
-
-    u = sim["/PartType0/InternalEnergies"][:]
-
-    # Compute the temperature
-    u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
-    u /= a[i] ** (3 * (gas_gamma - 1.0))
-    Temp = T(u)
-
-    # Gather statistics
-    T_median[i] = np.median(Temp)
-    T_mean[i] = Temp.mean()
-    T_std[i] = Temp.std()
-    T_log_mean[i] = np.log10(Temp).mean()
-    T_log_std[i] = np.log10(Temp).std()
-    T_min[i] = Temp.min()
-    T_max[i] = Temp.max()
-
-# CMB evolution
-a_evol = np.logspace(-3, 0, 60)
-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.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)
-
-# Cooling model
-if cooling_model == "Constant Lambda":
-    text(
-        1e-2,
-        6e4,
-        "$\Lambda_{\\rm const}/n_{\\rm H}^2 = %.1f\\times10^{%d}~[\\rm{cgs}]$"
-        % (Lambda / 10.0 ** (int(log10(Lambda))), log10(Lambda)),
-        fontsize=7,
-    )
-elif cooling_model == "EAGLE":
-    text(1e-2, 6e4, "EAGLE (Wiersma et al. 2009)")
-elif cooling_model == b"Grackle":
-    text(1e-2, 6e4, "Grackle (Smith et al. 2016)")
-else:
-    text(1e-2, 6e4, "No cooling")
-
-# 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.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.0 / (redshift_ticks + 1.0)
-
-xticks(a_ticks, redshift_labels)
-minorticks_off()
-
-xlabel("${\\rm Redshift}~z$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
-xlim(9e-3, 1.1)
-ylim(5, 2.5e7)
-
-savefig("Temperature_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py b/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py
deleted file mode 100644
index 3d75a4069081b4d76dfd61ad72c581fae7dca823..0000000000000000000000000000000000000000
--- a/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py
+++ /dev/null
@@ -1,195 +0,0 @@
-# Script for plotting energy evolution of uniform box of gas with single star in the
-# centre when running with stochastic energy injection. It also checks that the change
-# in total energy of the gas particles is within a specified tolerance from what is
-# expected based on the mass of the star particle (Note that this tolerance could be
-# somewhat high because of Poisson noise and the relatively small number of injection
-# events)
-
-import matplotlib
-
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-import os.path
-import numpy as np
-import glob
-
-# Number of snapshots and elements
-newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
-n_snapshots = (
-    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
-)
-n_elements = 9
-
-# 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.3,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.18,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.19,
-    "lines.markersize": 6,
-    "lines.linewidth": 2.0,
-}
-
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("stellar_evolution_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-star_initial_mass = sim["/PartType4/Masses"][0]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-# Units
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
-unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-unit_time_in_si = unit_time_in_cgs
-
-# Calculate solar mass in internal units
-const_solar_mass = 1.98848e33 / unit_mass_in_cgs
-
-# Define Gyr
-Gyr_in_cgs = 1e9 * 365 * 24 * 3600.0
-
-# Find out how many particles (gas and star) we have
-n_parts = sim["/Header"].attrs["NumPart_Total"][0]
-n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
-
-# Declare arrays for data
-masses = zeros((n_parts, n_snapshots))
-star_masses = zeros((n_sparts, n_snapshots))
-internal_energy = zeros((n_parts, n_snapshots))
-velocity_parts = zeros((n_parts, 3, n_snapshots))
-time = zeros(n_snapshots)
-
-# Read fields we are checking from snapshots
-# for i in [0,n_snapshots-1]:
-for i in range(n_snapshots):
-    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
-    print("reading snapshot " + str(i))
-    masses[:, i] = sim["/PartType0/Masses"]
-    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
-    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
-    time[i] = sim["/Header"].attrs["Time"][0]
-
-# Check that the total amount of enrichment is as expected.
-# Define tolerance. Note, relatively high value used due to
-# Poisson noise associated with stochastic energy injection.
-eps = 0.15
-
-# Stochastic heating
-vel2 = zeros((n_parts, n_snapshots))
-vel2[:, :] = (
-    velocity_parts[:, 0, :] * velocity_parts[:, 0, :]
-    + velocity_parts[:, 1, :] * velocity_parts[:, 1, :]
-    + velocity_parts[:, 2, :] * velocity_parts[:, 2, :]
-)
-total_kinetic_energy_cgs = (
-    np.sum(np.multiply(vel2, masses) * 0.5, axis=0) * unit_energy_in_cgs
-)
-total_energy_cgs = (
-    np.sum(np.multiply(internal_energy, masses), axis=0) * unit_energy_in_cgs
-)
-total_energy_released_cgs = (
-    total_energy_cgs[n_snapshots - 1]
-    - total_energy_cgs[0]
-    + total_kinetic_energy_cgs[n_snapshots - 1]
-    - total_kinetic_energy_cgs[0]
-)
-
-# Calculate energy released
-energy_per_sn = 1.0e51 / unit_energy_in_cgs
-SNIa_efficiency = 2.0e-3
-SNIa_timescale_Gyr = 2.0
-expected_energy_released_cgs = np.zeros(n_snapshots)
-for i in range(n_snapshots):
-    age_Gyr = time[i] * unit_time_in_cgs / Gyr_in_cgs
-    total_sn = (
-        SNIa_efficiency
-        * (1.0 - np.exp(-age_Gyr / SNIa_timescale_Gyr))
-        / const_solar_mass
-    )
-    expected_energy_released_cgs[i] = total_sn * energy_per_sn * unit_energy_in_cgs
-
-# Did we get it right?
-if (
-    abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots - 1])
-    / expected_energy_released_cgs[n_snapshots - 1]
-    < eps
-):
-    print(
-        "total stochastic energy release consistent with expectation. total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-    )
-else:
-    print(
-        "total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        + " energy change fraction of total "
-        + str(
-            total_energy_released_cgs
-            / (total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        )
-    )
-
-# Plot the energy evolution
-figure()
-subplot(111)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    total_energy_cgs
-    + total_kinetic_energy_cgs
-    - total_energy_cgs[0]
-    - total_kinetic_energy_cgs[0],
-    color="k",
-    linewidth=0.5,
-    label="SWIFT",
-)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    expected_energy_released_cgs,
-    color="r",
-    linewidth=0.5,
-    label="expected",
-)
-xlabel("Time (Gyr)")
-ylabel("Total energy (erg)")
-legend()
-savefig("continuous_energy_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py
deleted file mode 100644
index baaaca241a5bb290197763f385431b03f3e1c41a..0000000000000000000000000000000000000000
--- a/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py
+++ /dev/null
@@ -1,197 +0,0 @@
-# Script for plotting energy evolution of uniform box of gas with single star in the
-# centre when running with stochastic energy injection. It also checks that the change
-# in total energy of the gas particles is within a specified tolerance from what is
-# expected based on the mass of the star particle (Note that this tolerance could be
-# somewhat high because of Poisson noise and the relatively small number of injection
-# events)
-
-import matplotlib
-
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-import os.path
-import numpy as np
-import glob
-
-# Number of snapshots and elements
-newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
-n_snapshots = (
-    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
-)
-n_elements = 9
-
-# 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.3,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.18,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.19,
-    "lines.markersize": 6,
-    "lines.linewidth": 2.0,
-}
-
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("stellar_evolution_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-star_initial_mass = sim["/PartType4/Masses"][0]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-# Units
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
-unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-unit_time_in_si = unit_time_in_cgs
-
-# Calculate solar mass in internal units
-const_solar_mass = 1.98848e33 / unit_mass_in_cgs
-
-# Define Gyr
-Gyr_in_cgs = 1e9 * 365 * 24 * 3600.0
-
-# Find out how many particles (gas and star) we have
-n_parts = sim["/Header"].attrs["NumPart_Total"][0]
-n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
-
-# Declare arrays for data
-masses = zeros((n_parts, n_snapshots))
-star_masses = zeros((n_sparts, n_snapshots))
-internal_energy = zeros((n_parts, n_snapshots))
-velocity_parts = zeros((n_parts, 3, n_snapshots))
-time = zeros(n_snapshots)
-
-# Read fields we are checking from snapshots
-# for i in [0,n_snapshots-1]:
-for i in range(n_snapshots):
-    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
-    print("reading snapshot " + str(i))
-    masses[:, i] = sim["/PartType0/Masses"]
-    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
-    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
-    time[i] = sim["/Header"].attrs["Time"][0]
-
-# Check that the total amount of enrichment is as expected.
-# Define tolerance. Note, relatively high value used due to
-# Poisson noise associated with stochastic energy injection.
-eps = 0.15
-
-# Stochastic heating
-vel2 = zeros((n_parts, n_snapshots))
-vel2[:, :] = (
-    velocity_parts[:, 0, :] * velocity_parts[:, 0, :]
-    + velocity_parts[:, 1, :] * velocity_parts[:, 1, :]
-    + velocity_parts[:, 2, :] * velocity_parts[:, 2, :]
-)
-total_kinetic_energy_cgs = (
-    np.sum(np.multiply(vel2, masses) * 0.5, axis=0) * unit_energy_in_cgs
-)
-total_energy_cgs = (
-    np.sum(np.multiply(internal_energy, masses), axis=0) * unit_energy_in_cgs
-)
-total_energy_released_cgs = (
-    total_energy_cgs[n_snapshots - 1]
-    - total_energy_cgs[0]
-    + total_kinetic_energy_cgs[n_snapshots - 1]
-    - total_kinetic_energy_cgs[0]
-)
-
-# Calculate energy released
-num_SNII_per_msun = 1.73621e-02
-energy_per_sn = 1.0e51 / unit_energy_in_cgs
-expected_energy_released_cgs = np.zeros(n_snapshots)
-for i in range(n_snapshots):
-    if time[i] * unit_time_in_cgs / Gyr_in_cgs < 0.03:
-        expected_energy_released_cgs[i] = 0
-    else:
-        expected_energy_released_cgs[i] = (
-            num_SNII_per_msun
-            * star_initial_mass
-            / const_solar_mass
-            * energy_per_sn
-            * unit_energy_in_cgs
-        )
-
-# Did we get it right?
-if (
-    abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots - 1])
-    / expected_energy_released_cgs[n_snapshots - 1]
-    < eps
-):
-    print(
-        "total stochastic energy release consistent with expectation. total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-    )
-else:
-    print(
-        "total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        + " energy change fraction of total "
-        + str(
-            total_energy_released_cgs
-            / (total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        )
-    )
-
-# Plot the energy evolution
-figure()
-subplot(111)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    total_energy_cgs
-    + total_kinetic_energy_cgs
-    - total_energy_cgs[0]
-    - total_kinetic_energy_cgs[0],
-    color="k",
-    linewidth=0.5,
-    label="SWIFT",
-)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    expected_energy_released_cgs,
-    color="r",
-    linewidth=0.5,
-    label="expected",
-)
-xlabel("Time (Gyr)")
-ylabel("Total energy (erg)")
-legend()
-savefig("stochastic_energy_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
index 74c0a7fcdb884cf0bccd6838918855298845da3e..2ad291563f4d39d83b743b0b9262b5955c84731a 100644
--- a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
+++ b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
@@ -203,7 +203,7 @@ plt.subplot(222, frameon=False)
 
 plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time, fontsize=10)
 plt.plot([-0.49, 0.1], [0.82, 0.82], "k-", lw=1)
-plt.text(-0.49, 0.7, "$\\textsc{Swift}$ %s" % git, fontsize=10)
+plt.text(-0.49, 0.7, "$SWIFT$ %s" % git, fontsize=10)
 plt.text(-0.49, 0.6, scheme, fontsize=10)
 plt.text(-0.49, 0.5, kernel, fontsize=10)
 plt.text(-0.49, 0.4, chemistry + "'s Chemistry", fontsize=10)
diff --git a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
deleted file mode 100644
index 3d75a4069081b4d76dfd61ad72c581fae7dca823..0000000000000000000000000000000000000000
--- a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
+++ /dev/null
@@ -1,195 +0,0 @@
-# Script for plotting energy evolution of uniform box of gas with single star in the
-# centre when running with stochastic energy injection. It also checks that the change
-# in total energy of the gas particles is within a specified tolerance from what is
-# expected based on the mass of the star particle (Note that this tolerance could be
-# somewhat high because of Poisson noise and the relatively small number of injection
-# events)
-
-import matplotlib
-
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-import os.path
-import numpy as np
-import glob
-
-# Number of snapshots and elements
-newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
-n_snapshots = (
-    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
-)
-n_elements = 9
-
-# 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.3,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.18,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.19,
-    "lines.markersize": 6,
-    "lines.linewidth": 2.0,
-}
-
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("stellar_evolution_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-star_initial_mass = sim["/PartType4/Masses"][0]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-# Units
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
-unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-unit_time_in_si = unit_time_in_cgs
-
-# Calculate solar mass in internal units
-const_solar_mass = 1.98848e33 / unit_mass_in_cgs
-
-# Define Gyr
-Gyr_in_cgs = 1e9 * 365 * 24 * 3600.0
-
-# Find out how many particles (gas and star) we have
-n_parts = sim["/Header"].attrs["NumPart_Total"][0]
-n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
-
-# Declare arrays for data
-masses = zeros((n_parts, n_snapshots))
-star_masses = zeros((n_sparts, n_snapshots))
-internal_energy = zeros((n_parts, n_snapshots))
-velocity_parts = zeros((n_parts, 3, n_snapshots))
-time = zeros(n_snapshots)
-
-# Read fields we are checking from snapshots
-# for i in [0,n_snapshots-1]:
-for i in range(n_snapshots):
-    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
-    print("reading snapshot " + str(i))
-    masses[:, i] = sim["/PartType0/Masses"]
-    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
-    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
-    time[i] = sim["/Header"].attrs["Time"][0]
-
-# Check that the total amount of enrichment is as expected.
-# Define tolerance. Note, relatively high value used due to
-# Poisson noise associated with stochastic energy injection.
-eps = 0.15
-
-# Stochastic heating
-vel2 = zeros((n_parts, n_snapshots))
-vel2[:, :] = (
-    velocity_parts[:, 0, :] * velocity_parts[:, 0, :]
-    + velocity_parts[:, 1, :] * velocity_parts[:, 1, :]
-    + velocity_parts[:, 2, :] * velocity_parts[:, 2, :]
-)
-total_kinetic_energy_cgs = (
-    np.sum(np.multiply(vel2, masses) * 0.5, axis=0) * unit_energy_in_cgs
-)
-total_energy_cgs = (
-    np.sum(np.multiply(internal_energy, masses), axis=0) * unit_energy_in_cgs
-)
-total_energy_released_cgs = (
-    total_energy_cgs[n_snapshots - 1]
-    - total_energy_cgs[0]
-    + total_kinetic_energy_cgs[n_snapshots - 1]
-    - total_kinetic_energy_cgs[0]
-)
-
-# Calculate energy released
-energy_per_sn = 1.0e51 / unit_energy_in_cgs
-SNIa_efficiency = 2.0e-3
-SNIa_timescale_Gyr = 2.0
-expected_energy_released_cgs = np.zeros(n_snapshots)
-for i in range(n_snapshots):
-    age_Gyr = time[i] * unit_time_in_cgs / Gyr_in_cgs
-    total_sn = (
-        SNIa_efficiency
-        * (1.0 - np.exp(-age_Gyr / SNIa_timescale_Gyr))
-        / const_solar_mass
-    )
-    expected_energy_released_cgs[i] = total_sn * energy_per_sn * unit_energy_in_cgs
-
-# Did we get it right?
-if (
-    abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots - 1])
-    / expected_energy_released_cgs[n_snapshots - 1]
-    < eps
-):
-    print(
-        "total stochastic energy release consistent with expectation. total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-    )
-else:
-    print(
-        "total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        + " energy change fraction of total "
-        + str(
-            total_energy_released_cgs
-            / (total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        )
-    )
-
-# Plot the energy evolution
-figure()
-subplot(111)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    total_energy_cgs
-    + total_kinetic_energy_cgs
-    - total_energy_cgs[0]
-    - total_kinetic_energy_cgs[0],
-    color="k",
-    linewidth=0.5,
-    label="SWIFT",
-)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    expected_energy_released_cgs,
-    color="r",
-    linewidth=0.5,
-    label="expected",
-)
-xlabel("Time (Gyr)")
-ylabel("Total energy (erg)")
-legend()
-savefig("continuous_energy_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
index 2e6cced06ab2346efe368f52e4e90a045e855713..a16101a1104191252b5a13312a6f5f43cdf0610f 100644
--- a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
@@ -1,6 +1,22 @@
-import matplotlib
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
-matplotlib.use("Agg")
 from pylab import *
 import h5py
 import os.path
@@ -14,28 +30,6 @@ n_snapshots = (
 )
 n_elements = 9
 
-# 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.3,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.18,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.19,
-    "lines.markersize": 6,
-    "lines.linewidth": 2.0,
-}
-
-rcParams.update(params)
-
 # Read the simulation data
 sim = h5py.File("stellar_evolution_0000.hdf5", "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
@@ -94,16 +88,16 @@ for i in range(n_snapshots):
     sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
     print("reading snapshot " + str(i))
     abundances[:, :, i] = sim["/PartType0/ElementMassFractions"]
-    metallicity[:, i] = sim["/PartType0/Metallicities"]
+    metallicity[:, i] = sim["/PartType0/MetalMassFractions"]
     masses[:, i] = sim["/PartType0/Masses"]
     star_masses[:, i] = sim["/PartType4/Masses"]
-    mass_from_AGB[:, i] = sim["/PartType0/TotalMassFromAGB"]
-    metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFracFromAGB"]
-    mass_from_SNII[:, i] = sim["/PartType0/TotalMassFromSNII"]
-    metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFracFromSNII"]
-    mass_from_SNIa[:, i] = sim["/PartType0/TotalMassFromSNIa"]
-    metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFracFromSNIa"]
-    iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFracFromSNIa"]
+    mass_from_AGB[:, i] = sim["/PartType0/MassesFromAGB"]
+    metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFractionsFromAGB"]
+    mass_from_SNII[:, i] = sim["/PartType0/MassesFromSNII"]
+    metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFractionsFromSNII"]
+    mass_from_SNIa[:, i] = sim["/PartType0/MassesFromSNIa"]
+    metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFractionsFromSNIa"]
+    iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFractionsFromSNIa"]
     internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
     velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
     time[i] = sim["/Header"].attrs["Time"][0]
diff --git a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
deleted file mode 100644
index baaaca241a5bb290197763f385431b03f3e1c41a..0000000000000000000000000000000000000000
--- a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
+++ /dev/null
@@ -1,197 +0,0 @@
-# Script for plotting energy evolution of uniform box of gas with single star in the
-# centre when running with stochastic energy injection. It also checks that the change
-# in total energy of the gas particles is within a specified tolerance from what is
-# expected based on the mass of the star particle (Note that this tolerance could be
-# somewhat high because of Poisson noise and the relatively small number of injection
-# events)
-
-import matplotlib
-
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-import os.path
-import numpy as np
-import glob
-
-# Number of snapshots and elements
-newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
-n_snapshots = (
-    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
-)
-n_elements = 9
-
-# 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.3,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.18,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.21,
-    "figure.subplot.hspace": 0.19,
-    "lines.markersize": 6,
-    "lines.linewidth": 2.0,
-}
-
-rcParams.update(params)
-
-# Read the simulation data
-sim = h5py.File("stellar_evolution_0000.hdf5", "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-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]
-T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
-T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
-git = sim["Code"].attrs["Git Revision"]
-star_initial_mass = sim["/PartType4/Masses"][0]
-
-# Cosmological parameters
-H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
-gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
-
-# Units
-unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
-unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
-unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
-unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
-unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
-unit_length_in_si = 0.01 * unit_length_in_cgs
-unit_mass_in_si = 0.001 * unit_mass_in_cgs
-unit_time_in_si = unit_time_in_cgs
-
-# Calculate solar mass in internal units
-const_solar_mass = 1.98848e33 / unit_mass_in_cgs
-
-# Define Gyr
-Gyr_in_cgs = 1e9 * 365 * 24 * 3600.0
-
-# Find out how many particles (gas and star) we have
-n_parts = sim["/Header"].attrs["NumPart_Total"][0]
-n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
-
-# Declare arrays for data
-masses = zeros((n_parts, n_snapshots))
-star_masses = zeros((n_sparts, n_snapshots))
-internal_energy = zeros((n_parts, n_snapshots))
-velocity_parts = zeros((n_parts, 3, n_snapshots))
-time = zeros(n_snapshots)
-
-# Read fields we are checking from snapshots
-# for i in [0,n_snapshots-1]:
-for i in range(n_snapshots):
-    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
-    print("reading snapshot " + str(i))
-    masses[:, i] = sim["/PartType0/Masses"]
-    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
-    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
-    time[i] = sim["/Header"].attrs["Time"][0]
-
-# Check that the total amount of enrichment is as expected.
-# Define tolerance. Note, relatively high value used due to
-# Poisson noise associated with stochastic energy injection.
-eps = 0.15
-
-# Stochastic heating
-vel2 = zeros((n_parts, n_snapshots))
-vel2[:, :] = (
-    velocity_parts[:, 0, :] * velocity_parts[:, 0, :]
-    + velocity_parts[:, 1, :] * velocity_parts[:, 1, :]
-    + velocity_parts[:, 2, :] * velocity_parts[:, 2, :]
-)
-total_kinetic_energy_cgs = (
-    np.sum(np.multiply(vel2, masses) * 0.5, axis=0) * unit_energy_in_cgs
-)
-total_energy_cgs = (
-    np.sum(np.multiply(internal_energy, masses), axis=0) * unit_energy_in_cgs
-)
-total_energy_released_cgs = (
-    total_energy_cgs[n_snapshots - 1]
-    - total_energy_cgs[0]
-    + total_kinetic_energy_cgs[n_snapshots - 1]
-    - total_kinetic_energy_cgs[0]
-)
-
-# Calculate energy released
-num_SNII_per_msun = 1.73621e-02
-energy_per_sn = 1.0e51 / unit_energy_in_cgs
-expected_energy_released_cgs = np.zeros(n_snapshots)
-for i in range(n_snapshots):
-    if time[i] * unit_time_in_cgs / Gyr_in_cgs < 0.03:
-        expected_energy_released_cgs[i] = 0
-    else:
-        expected_energy_released_cgs[i] = (
-            num_SNII_per_msun
-            * star_initial_mass
-            / const_solar_mass
-            * energy_per_sn
-            * unit_energy_in_cgs
-        )
-
-# Did we get it right?
-if (
-    abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots - 1])
-    / expected_energy_released_cgs[n_snapshots - 1]
-    < eps
-):
-    print(
-        "total stochastic energy release consistent with expectation. total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-    )
-else:
-    print(
-        "total stochastic energy release "
-        + str(total_energy_released_cgs)
-        + " expected "
-        + str(expected_energy_released_cgs[n_snapshots - 1])
-        + " initial total internal energy "
-        + str(total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        + " energy change fraction of total "
-        + str(
-            total_energy_released_cgs
-            / (total_energy_cgs[0] + total_kinetic_energy_cgs[0])
-        )
-    )
-
-# Plot the energy evolution
-figure()
-subplot(111)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    total_energy_cgs
-    + total_kinetic_energy_cgs
-    - total_energy_cgs[0]
-    - total_kinetic_energy_cgs[0],
-    color="k",
-    linewidth=0.5,
-    label="SWIFT",
-)
-plot(
-    time * unit_time_in_cgs / Gyr_in_cgs,
-    expected_energy_released_cgs,
-    color="r",
-    linewidth=0.5,
-    label="expected",
-)
-xlabel("Time (Gyr)")
-ylabel("Total energy (erg)")
-legend()
-savefig("stochastic_energy_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index bdffbe180cc9ebbe990fe0e01ea4c2beda4a1c74..9c7268b87397c783f380dcd604ec2277dea3d489 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -33,27 +33,7 @@ import numpy as np
 import glob
 import os.path
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.05,
-    "figure.subplot.right": 0.995,
-    "figure.subplot.bottom": 0.06,
-    "figure.subplot.top": 0.92,
-    "figure.subplot.wspace": 0.25,
-    "figure.subplot.hspace": 0.2,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
+style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Number of snapshots and elements
 newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"))  # , key=os.path.getctime)
@@ -193,7 +173,7 @@ eagle_total_metal_mass_SNIa = data[:, 2] * stellar_mass / Msun_in_cgs * unit_mas
 
 
 # Plot the interesting quantities
-figure()
+figure(figsize=(7, 7 / 1.6))
 
 suptitle("Star metallicity Z = %.4f" % Z_star)
 
@@ -345,4 +325,6 @@ xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
 ylabel("Change in total metal mass of gas particles ${[\\rm M_\\odot]}$", labelpad=2)
 ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
+tight_layout()
+
 savefig("box_evolution_Z_%.4f.png" % (Z_star), dpi=200)
diff --git a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
index 4b7c6e005547073859b065fafb52cafd302cc514..5e72bc477022ec5aa56aa8202d6c36afbba9fffa 100644
--- a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
@@ -33,6 +33,8 @@ import numpy as np
 import glob
 import os.path
 
+style.use("../../../tools/stylesheets/mnras.mplstyle")
+
 # Function to find index in array a for each element in array b
 def find_indices(a, b):
     result = np.zeros(len(b))
@@ -41,28 +43,6 @@ def find_indices(a, b):
     return result
 
 
-# Plot parameters
-params = {
-    "axes.labelsize": 10,
-    "axes.titlesize": 10,
-    "font.size": 12,
-    "legend.fontsize": 12,
-    "xtick.labelsize": 10,
-    "ytick.labelsize": 10,
-    "text.usetex": True,
-    "figure.figsize": (9.90, 6.45),
-    "figure.subplot.left": 0.1,
-    "figure.subplot.right": 0.99,
-    "figure.subplot.bottom": 0.1,
-    "figure.subplot.top": 0.95,
-    "figure.subplot.wspace": 0.2,
-    "figure.subplot.hspace": 0.2,
-    "lines.markersize": 6,
-    "lines.linewidth": 3.0,
-}
-rcParams.update(params)
-
-
 # Number of snapshots and elements
 newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
 n_snapshots = (
@@ -143,7 +123,7 @@ for i in range(n_snapshots):
     P = sim["/PartType0/Pressures"][:]
     rho = sim["/PartType0/Densities"][:]
     mass = sim["/PartType0/Masses"][:]
-    metallicity = sim["/PartType0/Metallicities"][:]
+    metallicity = sim["/PartType0/MetalMassFractions"][:]
     internal_energy = sim["/PartType0/InternalEnergies"][:]
     IDs = sim["/PartType0/ParticleIDs"][:]
 
@@ -156,7 +136,7 @@ for i in range(n_snapshots):
 
 
 # Plot the interesting quantities
-figure()
+figure(figsize=(7, 7 / 1.6))
 
 # Radial velocity --------------------------------
 subplot(221)
@@ -218,4 +198,6 @@ xlabel("Time (Myr)", labelpad=0)
 ylabel("Metallicity", labelpad=2)
 ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
+tight_layout()
+
 savefig("particle_evolution.png", dpi=200)