diff --git a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
index 1f9ae6557d1a55bd75d3440679ae42d653da3399..4d24529b4590601325105964c094fd2e834cd6c6 100644
--- a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
+++ b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
@@ -1,184 +1,208 @@
 import matplotlib
+
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
-import h5py as h5
 import swiftsimio
 import sys
 import glob
 import unyt
 import numpy as np
 
+try:
+    plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+except:
+    print("Can't find Matplotlib stylesheet.")
+
 ## read command line arguments
 snapshot_name = sys.argv[1]
 
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 9,
-'legend.fontsize': 9,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': False,
-'figure.figsize' : (4.15,3.15),
-'figure.subplot.left'    : 0.12,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.12,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 2.,
-'text.latex.unicode': True
-}
-plt.rcParams.update(params)
-
 ################# Read in observational data
 
-data_schaye = np.genfromtxt("./datasets/schaye_et_al_2000_thermal_history.dat",skip_header = 2)
-data_walther = np.genfromtxt("./datasets/walther_et_al_2019_thermal_history.dat",skip_header = 2)
+data_schaye = np.genfromtxt(
+    "./datasets/schaye_et_al_2000_thermal_history.dat", skip_header=2
+)
+data_walther = np.genfromtxt(
+    "./datasets/walther_et_al_2019_thermal_history.dat", skip_header=2
+)
 
 data_schaye = data_schaye.T
 data_walther = data_walther.T
 
 schaye_z_lower_error = data_schaye[0] - data_schaye[1]
 schaye_z_upper_error = data_schaye[2] - data_schaye[0]
-schaye_T_lower_error = np.log10(data_schaye[3]*1.0e4) - np.log10(data_schaye[4]*1.0e4)
-schaye_T_upper_error = np.log10(data_schaye[5]*1.0e4) - np.log10(data_schaye[3]*1.0e4)
-walther_T_lower_error = np.log10(data_walther[1]*1.0e4) - np.log10(data_walther[2]*1.0e4)
-walther_T_upper_error = np.log10(data_walther[3]*1.0e4) - np.log10(data_walther[1]*1.0e4)
+
+schaye_T_lower_error = np.log10(data_schaye[3] * 1.0e4) - np.log10(
+    data_schaye[4] * 1.0e4
+)
+schaye_T_upper_error = np.log10(data_schaye[5] * 1.0e4) - np.log10(
+    data_schaye[3] * 1.0e4
+)
+walther_T_lower_error = np.log10(data_walther[1] * 1.0e4) - np.log10(
+    data_walther[2] * 1.0e4
+)
+walther_T_upper_error = np.log10(data_walther[3] * 1.0e4) - np.log10(
+    data_walther[1] * 1.0e4
+)
 
 ############### Read in simulation data
+density_units = "Solar_Mass / Mpc**3"
+temperature_units = "K"
 
 ## First, get list of all snapshots
-reg_exp = "%s*.hdf5" %snapshot_name
+reg_exp = "%s*.hdf5" % snapshot_name
 snap_list = glob.glob(reg_exp)
+snap_data = sorted(
+    [swiftsimio.load(snap) for snap in snap_list], key=lambda x: x.metadata.z
+)
 
-z = []
-T_mean = []
-T_std = []
-rho_mean = []
-rho_std = []
+z = np.empty(len(snap_data), dtype=np.float32)
+T_mean = unyt.unyt_array(np.empty_like(z), units=temperature_units)
+T_std = unyt.unyt_array(np.empty_like(z), units=temperature_units)
+rho_mean = unyt.unyt_array(np.empty_like(z), units=density_units)
+rho_std = unyt.unyt_array(np.empty_like(z), units=density_units)
 
 ## loop through list
-for snap in snap_list:
-    
-    # This loads all metadata but explicitly does _not_ read any particle data yet
-    data = swiftsimio.load(snap)
+for index, data in enumerate(snap_data):
+    # Stick redshift in a list
+    z[index] = data.metadata.z
 
-    # Get the redshift
-    z = np.append(z, data.metadata.z)
-        
-    # Convert gas temperatures and densities to right units
-    data.gas.temperatures.convert_to_cgs()
+    # Convert units to something sensible so `np.mean` doesn't freak out
+    data.gas.temperatures.convert_to_units(temperature_units)
+    data.gas.densities.convert_to_units(density_units)
 
     # Get mean and standard deviation of temperature
-    T_mean.append(np.mean(data.gas.temperatures) * data.gas.temperatures.units)
-    T_std.append(np.std(data.gas.temperatures) * data.gas.temperatures.units)
+    T_mean[index] = np.mean(data.gas.temperatures)
+    T_std[index] = np.std(data.gas.temperatures)
 
     # Get mean and standard deviation of density
-    rho_mean.append(np.mean(data.gas.densities) * data.gas.densities.units)
-    rho_std.append(np.std(data.gas.densities) * data.gas.densities.units)
-    
-## Turn into numpy arrays
-T_mean = np.array(T_mean) * data.gas.temperatures.units
-T_std = np.array(T_std) * data.gas.temperatures.units
-rho_mean = np.array(rho_mean) * data.gas.densities.units
-rho_std = np.array(rho_std) * data.gas.densities.units
+    rho_mean[index] = np.mean(data.gas.densities)
+    rho_std[index] = np.std(data.gas.densities)
+
 
 ## Put Density into units of mean baryon density today
 
 # first calculate rho_bar_0 from snapshot metadata
-### from the first snapshot, get cosmology information
-d = swiftsimio.load(snap_list[0])
-cosmology = d.metadata.cosmology
-H0 = cosmology["H0 [internal units]"] / (d.units.time)
+data = snap_data[0]
+cosmology = data.metadata.cosmology
+H0 = cosmology["H0 [internal units]"] / (data.units.time)
 Omega_bar = cosmology["Omega_b"]
 
 ### now calculate rho_bar_0 and divide through
-rho_bar_0 = 3.0 * H0**2 / (8. * np.pi * unyt.G) * Omega_bar 
+rho_bar_0 = 3.0 * H0 ** 2 / (8.0 * np.pi * unyt.G) * Omega_bar
 rho_mean /= rho_bar_0
 rho_std /= rho_bar_0
 
-### sort arrays into redshift order
-ind_sorted = np.argsort(z)
-z = z[ind_sorted]
-T_mean = T_mean[ind_sorted]
-T_std = T_std[ind_sorted]
-rho_mean = rho_mean[ind_sorted]
-rho_std = rho_std[ind_sorted]
-
 ### from the first snapshot, get code information
-d = swiftsimio.load(snap_list[0])
-code_info = d.metadata.code
-git_branch = code_info["Git Branch"].decode('UTF-8')
-git_revision = code_info["Git Revision"].decode('UTF-8')
-hydro_metadata = d.metadata.hydro_scheme
-scheme = hydro_metadata["Scheme"].decode('UTF-8')
-params = d.metadata.parameters
-
-subgrid_metadata = d.metadata.subgrid_scheme
-cooling_model = subgrid_metadata["Cooling Model"].decode('UTF-8')
-chemistry_model = subgrid_metadata["Chemistry Model"].decode('UTF-8')
-
-if cooling_model == 'EAGLE':
-    z_r_H = float(params['EAGLECooling:H_reion_z'])
-    H_heat_input = float(params['EAGLECooling:H_reion_eV_p_H'])
-    z_r_He_centre = float(params['EAGLECooling:He_reion_z_centre'])
-    z_r_He_sigma = float(params['EAGLECooling:He_reion_z_sigma'])
-    He_heat_input = float(params['EAGLECooling:He_reion_eV_p_H'])
-
-metallicity = "Unknown"
-if chemistry_model == 'EAGLE':
-    metallicity = float(params['EAGLEChemistry:init_abundance_metal'])
-    
+code_info = data.metadata.code
+git_branch = code_info["Git Branch"].decode("UTF-8")
+git_revision = code_info["Git Revision"].decode("UTF-8")
+
+hydro_metadata = data.metadata.hydro_scheme
+scheme = hydro_metadata["Scheme"].decode("UTF-8")
+params = data.metadata.parameters
+
+subgrid_metadata = data.metadata.subgrid_scheme
+cooling_model = subgrid_metadata["Cooling Model"].decode("UTF-8")
+chemistry_model = subgrid_metadata["Chemistry Model"].decode("UTF-8")
+
+if cooling_model == "EAGLE":
+    z_r_H = float(params["EAGLECooling:H_reion_z"])
+    H_heat_input = float(params["EAGLECooling:H_reion_eV_p_H"])
+    z_r_He_centre = float(params["EAGLECooling:He_reion_z_centre"])
+    z_r_He_sigma = float(params["EAGLECooling:He_reion_z_sigma"])
+    He_heat_input = float(params["EAGLECooling:He_reion_eV_p_H"])
+
+if chemistry_model == "EAGLE":
+    metallicity = float(params["EAGLEChemistry:init_abundance_metal"])
+else:
+    metallicity = "Unknown"
+
 # Make plot of temperature evolution  --------------------------------
-fig = plt.figure()
+fig, ax = plt.subplots(constrained_layout=True)
 
 # Plot sim properties
-if cooling_model == 'EAGLE':
-    plt.plot([z_r_H, z_r_H], [3.4, 4.4], 'k--', alpha=0.5, lw=0.7)
-    plt.text(z_r_H + 0.1, 3.55, "H reion.", rotation=90, alpha=0.5, fontsize=7, va="bottom")
-    plt.plot([z_r_He_centre, z_r_He_centre], [3.4, 4.4], 'k--', alpha=0.5, lw=0.7)
-    plt.text(z_r_He_centre + 0.1, 3.55, "HeII reion.", rotation=90, alpha=0.5, fontsize=7, va="bottom")
-    
+if cooling_model == "EAGLE":
+    ax.axvline(z_r_H, color="k", linestyle="--", alpha=0.5, lw=0.7, zorder=-1)
+    ax.text(
+        z_r_H + 0.1, 3.55, "H reion.", rotation=90, alpha=0.5, fontsize=7, va="bottom"
+    )
+    ax.axvline(z_r_He_centre, color="k", linestyle="--", alpha=0.5, lw=0.7, zorder=-1)
+    ax.text(
+        z_r_He_centre + 0.1,
+        3.55,
+        "HeII reion.",
+        rotation=90,
+        alpha=0.5,
+        fontsize=7,
+        va="bottom",
+    )
+
 # Plot observational data
-plt.errorbar(data_schaye[0],
-             np.log10(data_schaye[3]*1.0e4),
-             xerr = [schaye_z_lower_error,schaye_z_upper_error],
-             yerr = [schaye_T_lower_error,schaye_T_upper_error],
-             fmt='s', mec='0.3', color='0.3', markersize=4, markeredgewidth=0.5, linewidth=0.5, mfc='w', label="Schaye et al. (2000)")
-plt.errorbar(data_walther[0],
-             np.log10(data_walther[1]*1.0e4),
-             yerr = [walther_T_lower_error,walther_T_upper_error],
-              fmt='.', mec='0.3', color='0.3', markersize=7, markeredgewidth=0.5, linewidth=0.5, mfc='w', label = "Walther et al. (2019)")
+ax.errorbar(
+    data_schaye[0],
+    np.log10(data_schaye[3] * 1.0e4),
+    xerr=[schaye_z_lower_error, schaye_z_upper_error],
+    yerr=[schaye_T_lower_error, schaye_T_upper_error],
+    fmt="s",
+    mec="0.3",
+    color="0.3",
+    markersize=4,
+    markeredgewidth=0.5,
+    linewidth=0.5,
+    mfc="w",
+    label="Schaye et al. (2000)",
+)
+ax.errorbar(
+    data_walther[0],
+    np.log10(data_walther[1] * 1.0e4),
+    yerr=[walther_T_lower_error, walther_T_upper_error],
+    fmt=".",
+    mec="0.3",
+    color="0.3",
+    markersize=7,
+    markeredgewidth=0.5,
+    linewidth=0.5,
+    mfc="w",
+    label="Walther et al. (2019)",
+)
 
 # Plot simulation
-plt.plot(z, np.log10(T_mean))
+ax.plot(z, np.log10(T_mean))
 
 # Legend
-plt.legend(loc="upper right", frameon=True, fontsize=8, handletextpad=0.1, facecolor="w", edgecolor="w", framealpha=1.)
-plt.text(0.2, 4.8, "SWIFT %s \nCooling model: %s \n$\\rho=%.3f$ $\\Omega_{b}\\rho_{crit,0}$\n$Z=%s$"%(git_revision, cooling_model, rho_mean[-1], metallicity), va="top", ha="left", fontsize=8)
-
-plt.xlim(0, 12.2)
-plt.ylim(3.5,4.85)
-plt.xlabel("Redshift", labelpad=-0.5)
-plt.ylabel(r"$\log_{10}(T/K)$", labelpad=0)
-plt.savefig("Temperature_evolution.png", dpi=200)
-
-
-# Make plot of denisty evolution  --------------------------------
-plt.rcParams.update({'figure.subplot.left'    : 0.14})
-fig = plt.figure()
-
-plt.text(0.2, 1.011, "SWIFT %s"%git_revision, va="top", ha="left", fontsize=8)
-
-plt.fill_between(z,rho_mean - rho_std,rho_mean + rho_std,alpha = 0.5)
-plt.plot(z,rho_mean)
-
-plt.axhline(y = 1.0, linestyle = '--', color='k', alpha=0.5, lw=1.)
-
-plt.xlim(0.0,12.2)
-plt.ylabel(r"$\delta_b = \rho / \Omega_b\rho_{crit,0}$", labelpad=0.)
-plt.ylim(0.988,1.012)
-plt.yticks([0.99, 1., 1.01])
-plt.xlabel("Redshift", labelpad=-0.5)
-plt.savefig("Density_evolution.png", dpi=200)
+ax.legend(loc="upper right", frameon=True)
+ax.text(
+    0.025,
+    0.975,
+    "SWIFT %s \nCooling model: %s \n$\\rho=%.3f$ $\\Omega_{b}\\rho_{crit,0}$\n$Z=%s$"
+    % (git_revision, cooling_model, rho_mean[-1], metallicity),
+    va="top",
+    ha="left",
+    bbox={"fc": "white", "ec": "none"},
+    transform=plt.gca().transAxes,
+    zorder=1,
+)
+
+ax.set_xlim(0, 12.2)
+ax.set_ylim(3.5, 4.85)
+ax.set_xlabel("Redshift", labelpad=-0.5)
+ax.set_ylabel(r"$\log_{10}(T/K)$", labelpad=0)
+fig.savefig("Temperature_evolution.png")
+
+
+# Make plot of density evolution  --------------------------------
+fig, ax = plt.subplots(constrained_layout=True)
+
+ax.text(0.2, 1.011, "SWIFT %s" % git_revision, va="top", ha="left", fontsize=8)
+
+ax.fill_between(z, rho_mean - rho_std, rho_mean + rho_std, alpha=0.5)
+ax.plot(z, rho_mean)
+
+ax.axhline(y=1.0, linestyle="--", color="k", alpha=0.5, lw=1.0)
+
+ax.set_xlim(0.0, 12.2)
+ax.set_ylabel(r"$\delta_b = \rho / \Omega_b\rho_{crit,0}$", labelpad=0.0)
+ax.set_ylim(0.988, 1.012)
+ax.set_yticks([0.99, 1.0, 1.01])
+ax.set_xlabel("Redshift", labelpad=-0.5)
+fig.savefig("Density_evolution.png")
diff --git a/examples/HydroTests/SedovBlast_1D/plotSolution.py b/examples/HydroTests/SedovBlast_1D/plotSolution.py
index d82ad9e94610a916d900ea93863b2881757e73b3..231d32ab02adac8a2b785651a149a7f00101ebf2 100644
--- a/examples/HydroTests/SedovBlast_1D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_1D/plotSolution.py
@@ -1,31 +1,31 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # 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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# 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 2D Sedov blast wave.
 # The script works for a given initial box and dumped energy and computes the solution at a later time t.
 
 # Parameters
-rho_0 = 1.          # Background Density
-P_0 = 1.e-6         # Background Pressure
-E_0 = 1.            # Energy of the explosion
-gas_gamma = 5./3.   # Gas polytropic index
+rho_0 = 1.0  # Background Density
+P_0 = 1.0e-6  # Background Pressure
+E_0 = 1.0  # Energy of the explosion
+gas_gamma = 5.0 / 3.0  # Gas polytropic index
 
 
 # ---------------------------------------------------------------
@@ -33,38 +33,19 @@ gas_gamma = 5./3.   # Gas polytropic index
 # ---------------------------------------------------------------
 
 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.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
 
 # Read the simulation data
-sim = h5py.File("sedov_%04d.hdf5"%snap, "r")
+sim = h5py.File("sedov_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -73,21 +54,28 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
 r = abs(x)
-v_r = x * vel[:,0] / r
+v_r = x * vel[:, 0] / r
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
+
+try:
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
+    plot_diffusion = True
+except:
+    plot_diffusion = False
+
 try:
-    alpha = sim["/PartType0/ViscosityParameters"][:]
-    plot_alpha = True 
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
+    plot_viscosity = True
 except:
-    plot_alpha = False
+    plot_viscosity = False
 
 
 # Now, work our the solution....
@@ -95,25 +83,29 @@ except:
 from scipy.special import gamma as Gamma
 from numpy import *
 
-def calc_a(g,nu=3):
+
+def calc_a(g, nu=3):
     """ 
     exponents of the polynomials of the sedov solution
     g - the polytropic gamma
     nu - the dimension
     """
-    a = [0]*8
-   
+    a = [0] * 8
+
     a[0] = 2.0 / (nu + 2)
-    a[2] = (1-g) / (2*(g-1) + nu)
-    a[3] = nu / (2*(g-1) + nu)
-    a[5] = 2 / (g-2)
-    a[6] = g / (2*(g-1) + nu)
-   
-    a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
-    a[4] = a[1]*(nu+2) / (2-g)
-    a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+    a[2] = (1 - g) / (2 * (g - 1) + nu)
+    a[3] = nu / (2 * (g - 1) + nu)
+    a[5] = 2 / (g - 2)
+    a[6] = g / (2 * (g - 1) + nu)
+
+    a[1] = (((nu + 2) * g) / (2.0 + nu * (g - 1.0))) * (
+        (2.0 * nu * (2.0 - g)) / (g * (nu + 2.0) ** 2) - a[2]
+    )
+    a[4] = a[1] * (nu + 2) / (2 - g)
+    a[7] = (2 + nu * (g - 1)) * a[1] / (nu * (2 - g))
     return a
 
+
 def calc_beta(v, g, nu=3):
     """ 
     beta values for the sedov solution (coefficients of the polynomials of the similarity variables) 
@@ -122,15 +114,33 @@ def calc_beta(v, g, nu=3):
     nu- the dimension
     """
 
-    beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
-            -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
-     -0.5/(g-1)), dtype=float64)
+    beta = (
+        (nu + 2)
+        * (g + 1)
+        * array(
+            (
+                0.25,
+                (g / (g - 1)) * 0.5,
+                -(2 + nu * (g - 1))
+                / 2.0
+                / ((nu + 2) * (g + 1) - 2 * (2 + nu * (g - 1))),
+                -0.5 / (g - 1),
+            ),
+            dtype=float64,
+        )
+    )
 
     beta = outer(beta, v)
 
-    beta += (g+1) * array((0.0,  -1.0/(g-1),
-                           (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
-                           1.0/(g-1)), dtype=float64).reshape((4,1))
+    beta += (g + 1) * array(
+        (
+            0.0,
+            -1.0 / (g - 1),
+            (nu + 2) / ((nu + 2) * (g + 1) - 2.0 * (2 + nu * (g - 1))),
+            1.0 / (g - 1),
+        ),
+        dtype=float64,
+    ).reshape((4, 1))
 
     return beta
 
@@ -154,9 +164,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     a = calc_a(g, nu)
     beta = calc_beta(v, g=g, nu=nu)
     lbeta = log(beta)
-    
+
     r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
-    rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+    rho = ((g + 1.0) / (g - 1.0)) * exp(
+        a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2]
+    )
     p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
     u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
     p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
@@ -165,14 +177,17 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     # It is not a singularity, however, the gradients of our variables (wrt v) are.
     # r -> 0, u -> 0, rho -> 0, p-> constant
 
-    u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+    u[0] = 0.0
+    rho[0] = 0.0
+    r[0] = 0.0
+    p[0] = p[1]
 
     # volume of an n-sphere
     vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
 
     # note we choose to evaluate the integral in this way because the
-    # volumes of the first few elements (i.e near v=vmin) are shrinking 
-    # very slowly, so we dramatically improve the error convergence by 
+    # volumes of the first few elements (i.e near v=vmin) are shrinking
+    # very slowly, so we dramatically improve the error convergence by
     # finding the volumes exactly. This is most important for the
     # pressure integral, as this is on the order of the volume.
 
@@ -186,11 +201,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 
     # shock speed
     shock_speed = fac * (2.0 / (nu + 2))
-    rho_s = ((g + 1) / (g - 1)) * rho0                                                                            
+    rho_s = ((g + 1) / (g - 1)) * rho0
     r_s = shock_speed * t * (nu + 2) / 2.0
     p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
     u_s = (2.0 * shock_speed) / (g + 1)
-    
+
     r *= fac * t
     u *= fac
     p *= fac * fac * rho0
@@ -202,93 +217,115 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 1)
 
 # Append points for after the shock
-r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock * 1.5])
 rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
 P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
 v_s = np.insert(v_s, np.size(v_s), [0, 0])
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
-
+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()
+figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=1.0,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
 
 # Velocity profile --------------------------------
 subplot(231)
-plot(r, v_r, '.', color='r', ms=2.)
-plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
+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)
 
 # Density profile --------------------------------
 subplot(232)
-plot(r, rho, '.', color='r', ms=2.)
-plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
+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)
 
 # Pressure profile --------------------------------
 subplot(233)
-plot(r, P, '.', color='r', ms=2.)
-plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+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)
 
 # Internal energy profile -------------------------
 subplot(234)
-plot(r, u, '.', color='r', ms=2.)
-plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+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)
 
-# Entropy/alpha profile ---------------------------------
+# Entropy profile ---------------------------------
 subplot(235)
+xlabel("Radius $r$")
+if plot_diffusion or plot_viscosity:
+    if plot_diffusion:
+        plot(r, diffusion, **scatter_props)
+
+    if plot_viscosity:
+        plot(r, viscosity, **scatter_props)
 
-if plot_alpha:
-    plot(r, alpha, '.', color='r', ms=2.0)
-    plot([r_shock, r_shock], [-1, 2], "--", color="k", alpha=0.8, lw=1.2)
-    ylabel(r"${\rm{Viscosity}}~\alpha$", labelpad=0)
-    # Show location of shock
-    ylim(0, 2)
+    ylabel(r"Rate Coefficient $\alpha$", labelpad=0)
+    legend()
 else:
-    plot(r, S, '.', color='r', ms=2.0)
-    plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    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)
 
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
 xlim(0, 1.3 * r_shock)
-
 # Information -------------------------------------
 subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Sedov blast with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
-text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), 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)
+text_fontsize = 5
+
+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(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
 xlim(-0.5, 0.5)
 ylim(0, 1)
 xticks([])
 yticks([])
 
+tight_layout()
 
-savefig("Sedov.png", dpi=200)
-
-
-
+savefig("Sedov.png")
 
diff --git a/examples/HydroTests/SedovBlast_2D/plotSolution.py b/examples/HydroTests/SedovBlast_2D/plotSolution.py
index 6f504a09c9432368ce141ec0d28c28699f5ba7f3..d18b9de84dcce969d203fdf34dd1e08406e12cd0 100644
--- a/examples/HydroTests/SedovBlast_2D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_2D/plotSolution.py
@@ -1,31 +1,31 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # 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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# 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 2D Sedov blast wave.
 # The script works for a given initial box and dumped energy and computes the solution at a later time t.
 
 # Parameters
-rho_0 = 1.          # Background Density
-P_0 = 1.e-6         # Background Pressure
-E_0 = 1.            # Energy of the explosion
-gas_gamma = 5./3.   # Gas polytropic index
+rho_0 = 1.0  # Background Density
+P_0 = 1.0e-6  # Background Pressure
+E_0 = 1.0  # Energy of the explosion
+gas_gamma = 5.0 / 3.0  # Gas polytropic index
 
 
 # ---------------------------------------------------------------
@@ -33,39 +33,18 @@ gas_gamma = 5./3.   # Gas polytropic index
 # ---------------------------------------------------------------
 
 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 from scipy import stats
 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.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
-sim = h5py.File("sedov_%04d.hdf5"%snap, "r")
+sim = h5py.File("sedov_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -74,60 +53,94 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-y = pos[:,1] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
-r = sqrt(x**2 + y**2)
-v_r = (x * vel[:,0] + y * vel[:,1]) / r
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+r = sqrt(x ** 2 + y ** 2)
+v_r = (x * vel[:, 0] + y * vel[:, 1]) / r
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
+try:
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
+    plot_diffusion = True
+except:
+    plot_diffusion = False
+
+try:
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
+    plot_viscosity = True
+except:
+    plot_viscosity = False
+
 # Bin te data
-r_bin_edge = np.arange(0., 0.5, 0.01)
-r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
-rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
-v_bin,_,_ = stats.binned_statistic(r, v_r, statistic='mean', bins=r_bin_edge)
-P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
-S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
-u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
-rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
-v2_bin,_,_ = stats.binned_statistic(r, v_r**2, statistic='mean', bins=r_bin_edge)
-P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
-S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
-u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
-rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
-v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
-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)
+r_bin_edge = np.arange(0.0, 0.5, 0.01)
+r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)
+v_bin, _, _ = stats.binned_statistic(r, v_r, statistic="mean", bins=r_bin_edge)
+P_bin, _, _ = stats.binned_statistic(r, P, statistic="mean", bins=r_bin_edge)
+S_bin, _, _ = stats.binned_statistic(r, S, statistic="mean", bins=r_bin_edge)
+u_bin, _, _ = stats.binned_statistic(r, u, statistic="mean", bins=r_bin_edge)
+rho2_bin, _, _ = stats.binned_statistic(r, rho ** 2, statistic="mean", bins=r_bin_edge)
+v2_bin, _, _ = stats.binned_statistic(r, v_r ** 2, statistic="mean", bins=r_bin_edge)
+P2_bin, _, _ = stats.binned_statistic(r, P ** 2, statistic="mean", bins=r_bin_edge)
+S2_bin, _, _ = stats.binned_statistic(r, S ** 2, statistic="mean", bins=r_bin_edge)
+u2_bin, _, _ = stats.binned_statistic(r, u ** 2, statistic="mean", bins=r_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin ** 2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin ** 2)
+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)
+
+if plot_diffusion:
+    alpha_diff_bin, _, _ = stats.binned_statistic(
+        r, diffusion, statistic="mean", bins=r_bin_edge
+    )
+    alpha2_diff_bin, _, _ = stats.binned_statistic(
+        r, diffusion ** 2, statistic="mean", bins=r_bin_edge
+    )
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin ** 2)
+
+if plot_viscosity:
+    alpha_visc_bin, _, _ = stats.binned_statistic(
+        r, viscosity, statistic="mean", bins=r_bin_edge
+    )
+    alpha2_visc_bin, _, _ = stats.binned_statistic(
+        r, viscosity ** 2, statistic="mean", bins=r_bin_edge
+    )
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin ** 2)
 
 # Now, work our the solution....
 
 from scipy.special import gamma as Gamma
 from numpy import *
 
-def calc_a(g,nu=3):
+
+def calc_a(g, nu=3):
     """ 
     exponents of the polynomials of the sedov solution
     g - the polytropic gamma
     nu - the dimension
     """
-    a = [0]*8
-   
+    a = [0] * 8
+
     a[0] = 2.0 / (nu + 2)
-    a[2] = (1-g) / (2*(g-1) + nu)
-    a[3] = nu / (2*(g-1) + nu)
-    a[5] = 2 / (g-2)
-    a[6] = g / (2*(g-1) + nu)
-   
-    a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
-    a[4] = a[1]*(nu+2) / (2-g)
-    a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+    a[2] = (1 - g) / (2 * (g - 1) + nu)
+    a[3] = nu / (2 * (g - 1) + nu)
+    a[5] = 2 / (g - 2)
+    a[6] = g / (2 * (g - 1) + nu)
+
+    a[1] = (((nu + 2) * g) / (2.0 + nu * (g - 1.0))) * (
+        (2.0 * nu * (2.0 - g)) / (g * (nu + 2.0) ** 2) - a[2]
+    )
+    a[4] = a[1] * (nu + 2) / (2 - g)
+    a[7] = (2 + nu * (g - 1)) * a[1] / (nu * (2 - g))
     return a
 
+
 def calc_beta(v, g, nu=3):
     """ 
     beta values for the sedov solution (coefficients of the polynomials of the similarity variables) 
@@ -136,15 +149,33 @@ def calc_beta(v, g, nu=3):
     nu- the dimension
     """
 
-    beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
-            -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
-     -0.5/(g-1)), dtype=float64)
+    beta = (
+        (nu + 2)
+        * (g + 1)
+        * array(
+            (
+                0.25,
+                (g / (g - 1)) * 0.5,
+                -(2 + nu * (g - 1))
+                / 2.0
+                / ((nu + 2) * (g + 1) - 2 * (2 + nu * (g - 1))),
+                -0.5 / (g - 1),
+            ),
+            dtype=float64,
+        )
+    )
 
     beta = outer(beta, v)
 
-    beta += (g+1) * array((0.0,  -1.0/(g-1),
-                           (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
-                           1.0/(g-1)), dtype=float64).reshape((4,1))
+    beta += (g + 1) * array(
+        (
+            0.0,
+            -1.0 / (g - 1),
+            (nu + 2) / ((nu + 2) * (g + 1) - 2.0 * (2 + nu * (g - 1))),
+            1.0 / (g - 1),
+        ),
+        dtype=float64,
+    ).reshape((4, 1))
 
     return beta
 
@@ -168,9 +199,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     a = calc_a(g, nu)
     beta = calc_beta(v, g=g, nu=nu)
     lbeta = log(beta)
-    
+
     r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
-    rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+    rho = ((g + 1.0) / (g - 1.0)) * exp(
+        a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2]
+    )
     p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
     u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
     p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
@@ -179,14 +212,17 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     # It is not a singularity, however, the gradients of our variables (wrt v) are.
     # r -> 0, u -> 0, rho -> 0, p-> constant
 
-    u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+    u[0] = 0.0
+    rho[0] = 0.0
+    r[0] = 0.0
+    p[0] = p[1]
 
     # volume of an n-sphere
     vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
 
     # note we choose to evaluate the integral in this way because the
-    # volumes of the first few elements (i.e near v=vmin) are shrinking 
-    # very slowly, so we dramatically improve the error convergence by 
+    # volumes of the first few elements (i.e near v=vmin) are shrinking
+    # very slowly, so we dramatically improve the error convergence by
     # finding the volumes exactly. This is most important for the
     # pressure integral, as this is on the order of the volume.
 
@@ -200,11 +236,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 
     # shock speed
     shock_speed = fac * (2.0 / (nu + 2))
-    rho_s = ((g + 1) / (g - 1)) * rho0                                                                            
+    rho_s = ((g + 1) / (g - 1)) * rho0
     r_s = shock_speed * t * (nu + 2) / 2.0
     p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
     u_s = (2.0 * shock_speed) / (g + 1)
-    
+
     r *= fac * t
     u *= fac
     p *= fac * fac * rho0
@@ -216,89 +252,140 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 2)
 
 # Append points for after the shock
-r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock * 1.5])
 rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
 P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
 v_s = np.insert(v_s, np.size(v_s), [0, 0])
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
-
+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()
+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(r, v_r, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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{Radial~velocity}}~v_r$", labelpad=0)
+plot(r, v_r, **scatter_props)
+plot(r_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Radialvelocity $v_r$")
 xlim(0, 1.3 * r_shock)
 ylim(-0.2, 3.8)
 
 # Density profile --------------------------------
 subplot(232)
-plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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=2)
+plot(r, rho, **scatter_props)
+plot(r_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Density $\\rho$")
 xlim(0, 1.3 * r_shock)
 ylim(-0.2, 5.2)
 
 # Pressure profile --------------------------------
 subplot(233)
-plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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)
+plot(r, P, **scatter_props)
+plot(r_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Pressure $P$")
 xlim(0, 1.3 * r_shock)
 ylim(-1, 12.5)
 
 # Internal energy profile -------------------------
 subplot(234)
-plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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)
+plot(r, u, **scatter_props)
+plot(r_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Internal Energy $u$")
 xlim(0, 1.3 * r_shock)
 ylim(-2, 22)
 
 # Entropy profile ---------------------------------
 subplot(235)
-plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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=0)
-xlim(0, 1.3 * r_shock)
-ylim(-5, 50)
+xlabel("Radius $r$")
+if plot_diffusion or plot_viscosity:
+    if plot_diffusion:
+        plot(r, diffusion, **scatter_props)
+        errorbar(
+            r_bin,
+            alpha_diff_bin,
+            yerr=alpha_diff_sigma_bin,
+            **errorbar_props,
+            label="Diffusion",
+        )
+
+    if plot_viscosity:
+        plot(r, viscosity, **scatter_props)
+        errorbar(
+            r_bin,
+            alpha_visc_bin,
+            yerr=alpha_visc_sigma_bin,
+            label="Viscosity",
+            color="C4",
+            ms=binned_marker_size,
+            fmt=".",
+            lw=1.2,
+        )
+
+    ylabel(r"Rate Coefficient $\alpha$", labelpad=0)
+    legend()
+else:
+    plot(r, S, **scatter_props)
+    plot(r_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+    errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+    ylabel("Entropy $S$", labelpad=0)
+    ylim(-5, 50)
 
+xlim(0, 1.3 * r_shock)
 # Information -------------------------------------
 subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Sedov blast with  $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
-text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), 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)
+text_fontsize = 5
+
+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(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
 xlim(-0.5, 0.5)
 ylim(0, 1)
 xticks([])
 yticks([])
 
+tight_layout()
 
-savefig("Sedov.png", dpi=200)
-
-
-
-
+savefig("Sedov.png")
diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py
index fec4f1101406be5803a3f1601812d1cb85275409..b1c2278f800c89fa6d7190cafc6b79cdefd57443 100644
--- a/examples/HydroTests/SedovBlast_3D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py
@@ -1,31 +1,31 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # 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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# 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 3D Sedov blast wave.
 # The script works for a given initial box and dumped energy and computes the solution at a later time t.
 
 # Parameters
-rho_0 = 1.          # Background Density
-P_0 = 1.e-6         # Background Pressure
-E_0 = 1.0           # Energy of the explosion
-gas_gamma = 5./3.   # Gas polytropic index
+rho_0 = 1.0  # Background Density
+P_0 = 1.0e-6  # Background Pressure
+E_0 = 1.0  # Energy of the explosion
+gas_gamma = 5.0 / 3.0  # Gas polytropic index
 
 
 # ---------------------------------------------------------------
@@ -33,39 +33,19 @@ gas_gamma = 5./3.   # Gas polytropic index
 # ---------------------------------------------------------------
 
 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
+from numpy import power
 from scipy import stats
 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.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
-sim = h5py.File("sedov_%04d.hdf5"%snap, "r")
+sim = h5py.File("sedov_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -74,13 +54,13 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-y = pos[:,1] - boxSize / 2
-z = pos[:,2] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
-r = sqrt(x**2 + y**2 + z**2)
-v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+z = pos[:, 2] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+r = sqrt(x ** 2 + y ** 2 + z ** 2)
+v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
@@ -99,59 +79,70 @@ except:
     plot_viscosity = False
 
 # Bin te data
-r_bin_edge = np.arange(0., 0.5, 0.01)
-r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
-rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
-v_bin,_,_ = stats.binned_statistic(r, v_r, statistic='mean', bins=r_bin_edge)
-P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
-S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
-u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
-rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
-v2_bin,_,_ = stats.binned_statistic(r, v_r**2, statistic='mean', bins=r_bin_edge)
-P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
-S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
-u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
-rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
-v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
-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)
+r_bin_edge = np.arange(0.0, 0.5, 0.01)
+r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)
+v_bin, _, _ = stats.binned_statistic(r, v_r, statistic="mean", bins=r_bin_edge)
+P_bin, _, _ = stats.binned_statistic(r, P, statistic="mean", bins=r_bin_edge)
+S_bin, _, _ = stats.binned_statistic(r, S, statistic="mean", bins=r_bin_edge)
+u_bin, _, _ = stats.binned_statistic(r, u, statistic="mean", bins=r_bin_edge)
+rho2_bin, _, _ = stats.binned_statistic(r, rho ** 2, statistic="mean", bins=r_bin_edge)
+v2_bin, _, _ = stats.binned_statistic(r, v_r ** 2, statistic="mean", bins=r_bin_edge)
+P2_bin, _, _ = stats.binned_statistic(r, P ** 2, statistic="mean", bins=r_bin_edge)
+S2_bin, _, _ = stats.binned_statistic(r, S ** 2, statistic="mean", bins=r_bin_edge)
+u2_bin, _, _ = stats.binned_statistic(r, u ** 2, statistic="mean", bins=r_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin ** 2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin ** 2)
+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)
 
 if plot_diffusion:
-    alpha_diff_bin,_,_ = stats.binned_statistic(r, diffusion, statistic='mean', bins=r_bin_edge)
-    alpha2_diff_bin,_,_ = stats.binned_statistic(r, diffusion**2, statistic='mean', bins=r_bin_edge)
-    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2)
+    alpha_diff_bin, _, _ = stats.binned_statistic(
+        r, diffusion, statistic="mean", bins=r_bin_edge
+    )
+    alpha2_diff_bin, _, _ = stats.binned_statistic(
+        r, diffusion ** 2, statistic="mean", bins=r_bin_edge
+    )
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin ** 2)
 
 if plot_viscosity:
-    alpha_visc_bin,_,_ = stats.binned_statistic(r, viscosity, statistic='mean', bins=r_bin_edge)
-    alpha2_visc_bin,_,_ = stats.binned_statistic(r, viscosity**2, statistic='mean', bins=r_bin_edge)
-    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2)
+    alpha_visc_bin, _, _ = stats.binned_statistic(
+        r, viscosity, statistic="mean", bins=r_bin_edge
+    )
+    alpha2_visc_bin, _, _ = stats.binned_statistic(
+        r, viscosity ** 2, statistic="mean", bins=r_bin_edge
+    )
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin ** 2)
 
 
 # Now, work our the solution....
 
 from scipy.special import gamma as Gamma
-from numpy import *
 
-def calc_a(g,nu=3):
+
+def calc_a(g, nu=3):
     """ 
     exponents of the polynomials of the sedov solution
     g - the polytropic gamma
     nu - the dimension
     """
-    a = [0]*8
-   
+    a = [0] * 8
+
     a[0] = 2.0 / (nu + 2)
-    a[2] = (1-g) / (2*(g-1) + nu)
-    a[3] = nu / (2*(g-1) + nu)
-    a[5] = 2 / (g-2)
-    a[6] = g / (2*(g-1) + nu)
-   
-    a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
-    a[4] = a[1]*(nu+2) / (2-g)
-    a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+    a[2] = (1 - g) / (2 * (g - 1) + nu)
+    a[3] = nu / (2 * (g - 1) + nu)
+    a[5] = 2 / (g - 2)
+    a[6] = g / (2 * (g - 1) + nu)
+
+    a[1] = (((nu + 2) * g) / (2.0 + nu * (g - 1.0))) * (
+        (2.0 * nu * (2.0 - g)) / (g * (nu + 2.0) ** 2) - a[2]
+    )
+    a[4] = a[1] * (nu + 2) / (2 - g)
+    a[7] = (2 + nu * (g - 1)) * a[1] / (nu * (2 - g))
     return a
 
+
 def calc_beta(v, g, nu=3):
     """ 
     beta values for the sedov solution (coefficients of the polynomials of the similarity variables) 
@@ -160,15 +151,33 @@ def calc_beta(v, g, nu=3):
     nu- the dimension
     """
 
-    beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
-            -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
-     -0.5/(g-1)), dtype=float64)
+    beta = (
+        (nu + 2)
+        * (g + 1)
+        * array(
+            (
+                0.25,
+                (g / (g - 1)) * 0.5,
+                -(2 + nu * (g - 1))
+                / 2.0
+                / ((nu + 2) * (g + 1) - 2 * (2 + nu * (g - 1))),
+                -0.5 / (g - 1),
+            ),
+            dtype=float64,
+        )
+    )
 
     beta = outer(beta, v)
 
-    beta += (g+1) * array((0.0,  -1.0/(g-1),
-                           (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
-                           1.0/(g-1)), dtype=float64).reshape((4,1))
+    beta += (g + 1) * array(
+        (
+            0.0,
+            -1.0 / (g - 1),
+            (nu + 2) / ((nu + 2) * (g + 1) - 2.0 * (2 + nu * (g - 1))),
+            1.0 / (g - 1),
+        ),
+        dtype=float64,
+    ).reshape((4, 1))
 
     return beta
 
@@ -192,9 +201,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     a = calc_a(g, nu)
     beta = calc_beta(v, g=g, nu=nu)
     lbeta = log(beta)
-    
+
     r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
-    rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+    rho = ((g + 1.0) / (g - 1.0)) * exp(
+        a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2]
+    )
     p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
     u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
     p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
@@ -203,14 +214,17 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
     # It is not a singularity, however, the gradients of our variables (wrt v) are.
     # r -> 0, u -> 0, rho -> 0, p-> constant
 
-    u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+    u[0] = 0.0
+    rho[0] = 0.0
+    r[0] = 0.0
+    p[0] = p[1]
 
     # volume of an n-sphere
     vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
 
     # note we choose to evaluate the integral in this way because the
-    # volumes of the first few elements (i.e near v=vmin) are shrinking 
-    # very slowly, so we dramatically improve the error convergence by 
+    # volumes of the first few elements (i.e near v=vmin) are shrinking
+    # very slowly, so we dramatically improve the error convergence by
     # finding the volumes exactly. This is most important for the
     # pressure integral, as this is on the order of the volume.
 
@@ -224,11 +238,11 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 
     # shock speed
     shock_speed = fac * (2.0 / (nu + 2))
-    rho_s = ((g + 1) / (g - 1)) * rho0                                                                            
+    rho_s = ((g + 1) / (g - 1)) * rho0
     r_s = shock_speed * t * (nu + 2) / 2.0
     p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
     u_s = (2.0 * shock_speed) / (g + 1)
-    
+
     r *= fac * t
     u *= fac
     p *= fac * fac * rho0
@@ -240,103 +254,141 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
 r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 3)
 
 # Append points for after the shock
-r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock * 1.5])
 rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
 P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
 v_s = np.insert(v_s, np.size(v_s), [0, 0])
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
-
+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()
+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_r, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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{Radial~velocity}}~v_r$", labelpad=0)
+plot(r, v_r, **scatter_props)
+plot(r_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Radialvelocity $v_r$")
 xlim(0, 1.3 * r_shock)
 ylim(-0.2, 3.8)
 
 # Density profile --------------------------------
 subplot(232)
-plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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=2)
+plot(r, rho, **scatter_props)
+plot(r_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Density $\\rho$")
 xlim(0, 1.3 * r_shock)
 ylim(-0.2, 5.2)
 
 # Pressure profile --------------------------------
 subplot(233)
-plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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)
+plot(r, P, **scatter_props)
+plot(r_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Pressure $P$")
 xlim(0, 1.3 * r_shock)
 ylim(-1, 12.5)
 
 # Internal energy profile -------------------------
 subplot(234)
-plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_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)
+plot(r, u, **scatter_props)
+plot(r_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+xlabel("Radius $r$")
+ylabel("Internal Energy $u$")
 xlim(0, 1.3 * r_shock)
 ylim(-2, 22)
 
 # Entropy profile ---------------------------------
 subplot(235)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-
-
+xlabel("Radius $r$")
 if plot_diffusion or plot_viscosity:
     if plot_diffusion:
-        plot(r, diffusion, ".", color='r', ms=0.5, alpha=0.2)
-        errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion")
+        plot(r, diffusion, **scatter_props)
+        errorbar(
+            r_bin,
+            alpha_diff_bin,
+            yerr=alpha_diff_sigma_bin,
+            **errorbar_props,
+            label="Diffusion",
+        )
 
     if plot_viscosity:
-        plot(r, viscosity, ".", color='g', ms=0.5, alpha=0.2)
-        errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity")
-
-    ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0)
+        plot(r, viscosity, **scatter_props)
+        errorbar(
+            r_bin,
+            alpha_visc_bin,
+            yerr=alpha_visc_sigma_bin,
+            label="Viscosity",
+            color="C4",
+            ms=binned_marker_size,
+            fmt=".",
+            lw=1.2,
+        )
+
+    ylabel(r"Rate Coefficient $\alpha$", labelpad=0)
     legend()
 else:
-    plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
-    plot(r_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)
-    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    plot(r, S, **scatter_props)
+    plot(r_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+    errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+    ylabel("Entropy $S$", labelpad=0)
     ylim(-5, 50)
 
 xlim(0, 1.3 * r_shock)
 # Information -------------------------------------
 subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Sedov blast with  $\\gamma=%.3f$ in 3D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
-text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), 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)
+text_fontsize = 5
+
+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(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
 xlim(-0.5, 0.5)
 ylim(0, 1)
 xticks([])
 yticks([])
 
+tight_layout()
 
-savefig("Sedov.png", dpi=200)
-
-
-
+savefig("Sedov.png")
 
diff --git a/tools/stylesheets/README b/tools/stylesheets/README
new file mode 100644
index 0000000000000000000000000000000000000000..4b72981e6d55667be54b8389001502b68b59dc99
--- /dev/null
+++ b/tools/stylesheets/README
@@ -0,0 +1,2 @@
+This folder contains matplotlib stylesheets for use in python
+plotting scripts.
diff --git a/tools/stylesheets/mnras.mplstyle b/tools/stylesheets/mnras.mplstyle
new file mode 100644
index 0000000000000000000000000000000000000000..2b1ff222b5fb3a20782896f23a91275957d9d269
--- /dev/null
+++ b/tools/stylesheets/mnras.mplstyle
@@ -0,0 +1,30 @@
+# Line Properties
+lines.linewidth: 2
+
+# Font options
+font.size: 8
+font.family: STIXGeneral
+mathtext.fontset: stix
+
+# LaTeX options
+text.usetex: False
+
+# Legend options
+legend.frameon: False
+legend.fontsize: 8
+
+# Figure options for publications
+figure.dpi: 300
+figure.figsize: 3.321, 3.0 # Correct width for MNRAS
+
+# Histogram options
+hist.bins: auto
+
+# Ticks inside plots; more space devoted to plot.
+xtick.direction: in
+ytick.direction: in
+# Always plot ticks on top of data
+axes.axisbelow: False
+
+# Setup colours
+axes.prop_cycle: cycler('color', ['00AEFF','FF8C40','CC4314',  '5766B3', '68246D'])