diff --git a/examples/HydrostaticHalo/density_profile.py b/examples/HydrostaticHalo/density_profile.py
index d0afd399f951cf3b727e869ca8571a3a802c2e8d..5248587ec343d3c0ffe2cef0cbd8716b9a1e055c 100644
--- a/examples/HydrostaticHalo/density_profile.py
+++ b/examples/HydrostaticHalo/density_profile.py
@@ -1,6 +1,27 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys
 
 #for the plotting
@@ -46,7 +67,8 @@ for i in range(n_snaps):
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
@@ -63,58 +85,46 @@ for i in range(n_snaps):
     bin_width = bin_edges[1] - bin_edges[0]
     hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin
 
-#find the mass in each radial bin
 
+    #find the mass in each radial bin
     mass_dset = f["PartType0/Masses"]
-#mass of each particles should be equal
+
+    #mass of each particles should be equal
     part_mass = np.array(mass_dset)[0]
     part_mass_cgs = part_mass * unit_mass_cgs
     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
 
     mass_hist = hist * part_mass_over_virial_mass
     radial_bin_mids = np.linspace(bin_width/2.,max_r - bin_width/2.,n_radial_bins)
-#volume in each radial bin
+
+    #volume in each radial bin
     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
 
-#now divide hist by the volume so we have a density in each bin
 
+    #now divide hist by the volume so we have a density in each bin
     density = mass_hist / volume
 
-    ##read the densities
-
-    # density_dset = f["PartType0/Density"]
-    # density = np.array(density_dset)
-    # density_cgs = density * unit_mass_cgs / unit_length_cgs**3
-    # rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
-
     t = np.linspace(10./n_radial_bins,10.0,1000)
     rho_analytic = t**(-2)/(4.*np.pi)
 
-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
 
-    #initial analytic density profile
-    
+    #initial analytic density profile    
     if (i == 0):
         r_0 = radial_bin_mids[0]
         rho_0 = density[0]
-
         rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
-    plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
-    #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$\rho / \rho_{init})$")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    #plt.ylim((1.e-2,1.e1))
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
-    plt.xlim((radial_bin_mids[0],max_r))
-    plt.ylim((0,20))
-    plt.plot((0,max_r),(1,1))
-    #plt.xscale('log')
-    #plt.yscale('log')
-    plt.legend(loc = "upper right")
+
+    figure()
+    plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
+    #plot(t,rho_analytic,label = "Initial analytic density profile")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$\rho / \rho_{init}$")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    xlim((radial_bin_mids[0],max_r))
+    ylim((0,2))
+    plot((0,max_r),(1,1))
+    legend(loc = "upper right")
     plot_filename = "./plots/density_profile/density_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()
 
diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py
index ea52cf8fc5fd098a46f05eaa58494529a868000c..f1be049adb8e972f89fd9ffe86106b1b9f3b19dc 100644
--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
@@ -1,6 +1,27 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys
 
 def do_binning(x,y,x_bin_edges):
@@ -48,8 +69,6 @@ unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
 v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
-#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
-#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
@@ -64,7 +83,8 @@ for i in range(n_snaps):
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
@@ -75,11 +95,11 @@ for i in range(n_snaps):
     radius_cgs = radius*unit_length_cgs
     radius_over_virial_radius = radius_cgs / r_vir_cgs
 
-#get the internal energies
+    #get the internal energies
     u_dset = f["PartType0/InternalEnergy"]
     u = np.array(u_dset)
 
-#make dimensionless
+    #make dimensionless
     u /= v_c**2/(2. * (gamma - 1.))
     r = radius_over_virial_radius
 
@@ -90,21 +110,16 @@ for i in range(n_snaps):
     radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
     binned_u = u_totals / hist
 
-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
-
-    plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    #plt.plot((0,1),(1,1),label = "Analytic Solution")
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
-    plt.legend(loc = "lower right")
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,2))
+    figure()
+    plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
+    legend(loc = "lower right")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    ylim((0,2))
     plot_filename = "./plots/internal_energy/internal_energy_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()
 
 
         
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index d23ead6a67f43c9d19d76a797e72d050a3978d61..82584282559c1fceb0492aada671ff83fb74c924 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -1,11 +1,14 @@
 #!/bin/bash
 
 # Generate the initial conditions if they are not present.
-echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 100000
+if [ ! -e Hydrostatic.hdf5 ]
+then
+    echo "Generating initial conditions for the isothermal potential box example..."
+    python makeIC.py 100000
+fi
 
 # Run for 10 dynamical times
-../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
+../swift -g -s -t 1 hydrostatic.yml 2>&1 | tee output.log
 
 echo "Plotting density profiles"
 mkdir plots
diff --git a/examples/HydrostaticHalo/test_energy_conservation.py b/examples/HydrostaticHalo/test_energy_conservation.py
index ca091050c4127d11a37a2cc7504e42d244031e25..8368d475813d248ca93c12e46737b062752ab779 100644
--- a/examples/HydrostaticHalo/test_energy_conservation.py
+++ b/examples/HydrostaticHalo/test_energy_conservation.py
@@ -1,6 +1,27 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys
 
 n_snaps = int(sys.argv[1])
@@ -24,7 +45,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
@@ -45,7 +66,8 @@ for i in range(n_snaps):
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
@@ -73,7 +95,6 @@ for i in range(n_snaps):
     internal_energy_array = np.append(internal_energy_array,total_internal_energy)
 
 #put energies in units of v_c^2 and rescale by number of particles
-
 pe = potential_energy_array / (N*v_c**2)
 ke = kinetic_energy_array / (N*v_c**2)
 ie = internal_energy_array / (N*v_c**2)
@@ -82,14 +103,15 @@ te = pe + ke + ie
 dyn_time_cgs = r_vir_cgs / v_c_cgs
 time_array = time_array_cgs / dyn_time_cgs
 
-plt.plot(time_array,ke,label = "Kinetic Energy")
-plt.plot(time_array,pe,label = "Potential Energy")
-plt.plot(time_array,ie,label = "Internal Energy")
-plt.plot(time_array,te,label = "Total Energy")
-plt.legend(loc = "lower right")
-plt.xlabel(r"$t / t_{dyn}$")
-plt.ylabel(r"$E / v_c^2$")
-plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
-plt.ylim((-2,2))
-plt.savefig("energy_conservation.png",format = 'png')
+figure()
+plot(time_array,ke,label = "Kinetic Energy")
+plot(time_array,pe,label = "Potential Energy")
+plot(time_array,ie,label = "Internal Energy")
+plot(time_array,te,label = "Total Energy")
+legend(loc = "lower right")
+xlabel(r"$t / t_{dyn}$")
+ylabel(r"$E / v_c^2$")
+title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
+ylim((-2,2))
+savefig("energy_conservation.png",format = 'png')
 
diff --git a/examples/HydrostaticHalo/velocity_profile.py b/examples/HydrostaticHalo/velocity_profile.py
index 9133195d942233514148aa419003ee0ab7923494..f8f607362846a323937a9203dab8bc228f52a149 100644
--- a/examples/HydrostaticHalo/velocity_profile.py
+++ b/examples/HydrostaticHalo/velocity_profile.py
@@ -1,6 +1,27 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+ # 
+ ##############################################################################
+
 import numpy as np
 import h5py as h5
-import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
 import sys
 
 def do_binning(x,y,x_bin_edges):
@@ -62,7 +83,8 @@ for i in range(n_snaps):
     f = h5.File(filename,'r')
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+
+    #translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
@@ -73,16 +95,15 @@ for i in range(n_snaps):
     radius_cgs = radius*unit_length_cgs
     radius_over_virial_radius = radius_cgs / r_vir_cgs
 
-#get the internal energies
+    #get the internal energies
     vel_dset = f["PartType0/Velocities"]
     vel = np.array(vel_dset)
 
-#make dimensionless
+    #make dimensionless
     vel /= v_c
     r = radius_over_virial_radius
 
     #find radial component of velocity
-
     v_r = np.zeros(r.size)
     for j in range(r.size):
         v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
@@ -94,18 +115,13 @@ for i in range(n_snaps):
     radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
     binned_v_r = v_r_totals / hist
 
-    #calculate cooling radius
-
-    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
-
-    plt.plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
-    #plt.plot((0,1),(1,1),label = "Analytic Solution")
-    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
-    plt.legend(loc = "upper right")
-    plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$v_r / v_c$")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,2))
+    figure()
+    plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
+    legend(loc = "upper right")
+    xlabel(r"$r / r_{vir}$")
+    ylabel(r"$v_r / v_c$")
+    title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    ylim((-1,1))
     plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
+    savefig(plot_filename,format = "png")
+    close()