From d1eb8dc367548d19c6efba7c0fe19e6aad2dfce0 Mon Sep 17 00:00:00 2001
From: Folkert Nobels <nobels@strw.leidenuniv.nl>
Date: Wed, 30 Jan 2019 09:32:19 +0100
Subject: [PATCH] Update plotSolutions of star formation example and remove the
 two other examples

---
 .../plotSolution.py                           |  61 +-
 .../KS_law.py                                 | 134 -----
 .../StarFormationHistory.py                   |  36 --
 .../gasdensity.py                             | 563 ------------------
 .../IsolatedGalaxy_starformation_25p/getIC.sh |   2 -
 .../isolated_galaxy.yml                       | 115 ----
 .../IsolatedGalaxy_starformation_25p/run.sh   |   3 -
 .../gasdensity.py                             | 563 ------------------
 .../getIC.sh                                  |   2 -
 .../isolated_galaxy.yml                       | 114 ----
 .../run.sh                                    |   8 -
 11 files changed, 58 insertions(+), 1543 deletions(-)
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p/KS_law.py
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p/StarFormationHistory.py
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p/gasdensity.py
 delete mode 100755 examples/IsolatedGalaxy_starformation_25p/getIC.sh
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p/isolated_galaxy.yml
 delete mode 100755 examples/IsolatedGalaxy_starformation_25p/run.sh
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p_intermediate/gasdensity.py
 delete mode 100755 examples/IsolatedGalaxy_starformation_25p_intermediate/getIC.sh
 delete mode 100644 examples/IsolatedGalaxy_starformation_25p_intermediate/isolated_galaxy.yml
 delete mode 100755 examples/IsolatedGalaxy_starformation_25p_intermediate/run.sh

diff --git a/examples/IsolatedGalaxy_starformation/plotSolution.py b/examples/IsolatedGalaxy_starformation/plotSolution.py
index 4a0ad38f2a..fc3613b7bd 100644
--- a/examples/IsolatedGalaxy_starformation/plotSolution.py
+++ b/examples/IsolatedGalaxy_starformation/plotSolution.py
@@ -39,6 +39,9 @@ k_in_cgs = 1.38064852e-16
 mH_in_cgs = 1.6737236e-24
 year_in_cgs = 3600.0 * 24 * 365.0
 Msun_in_cgs = 1.98848e33
+G_in_cgs = 6.67259e-8
+pc_in_cgs = 3.08567758e18
+Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
 
 # Gemoetry info
 boxsize = f["/Header"].attrs["BoxSize"]
@@ -49,6 +52,9 @@ unit_length_in_cgs = f["/Units"].attrs["Unit length in cgs (U_L)"]
 unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
 unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
 
+# Calculate Gravitational constant in internal units
+G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
+
 # Read parameters of the model
 KS_law_slope = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawExponent"])
 KS_law_norm = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawCoeff_MSUNpYRpKPC2"])
@@ -60,7 +66,18 @@ KS_thresh_max_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_max_norm_HpCM
 KS_gamma_effective = float(
     f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
 )
+KS_high_den_thresh = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawHighDens_thresh_HpCM3"])
+KS_law_slope_high_den = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawHighDensExponent"])
 EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
+EAGLEfloor_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
+EAGLEfloor_T = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
+EAGLEfloor_cool_rho = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
+EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
+EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
+
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
+gamma = 5./3.
+KS_press_norm = k_in_cgs * EAGLEfloor_T * EAGLEfloor_rho_norm
 
 # Read gas properties
 gas_pos = f["/PartType0/Coordinates"][:, :]
@@ -71,6 +88,11 @@ gas_SFR = f["/PartType0/SFR"][:]
 gas_XH = f["/PartType0/ElementAbundance"][:, 0]
 gas_Z = f["/PartType0/Metallicity"][:]
 gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_sSFR = f["/PartType0/sSFR"][:]
+
+# Read the Star properties
+stars_pos = f["/PartType4/Coordinates"][:, :]
+stars_BirthDensity = f["/PartType4/BirthDensity"][:]
 
 # Centre the box
 gas_pos[:, 0] -= centre[0]
@@ -90,12 +112,15 @@ gas_nH = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
 gas_nH /= mH_in_cgs
 gas_nH *= gas_XH
 
+stars_BirthDensity = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
+stars_BirthDensity /= mH_in_cgs
+stars_BirthDensity *= gas_XH
+
 # Equations of state
 eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = eos_cool_rho ** 0.0 * 8000.0
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * eos_cool_rho **( EAGLEfloor_cool_gamma_effective - 1.0 )
 eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = (eos_Jeans_rho / 10 ** (-1)) ** (1.0 / 3.0) * 8000.0
-
+eos_Jeans_T = EAGLEfloor_T * (eos_Jeans_rho / EAGLEfloor_rho_norm) ** (KS_gamma_effective - 1.0 ) 
 # Plot the phase space diagram
 figure()
 subplot(111, xscale="log", yscale="log")
@@ -135,6 +160,36 @@ savefig("rho_SFR.png", dpi=200)
 
 ########################################################################3
 
+# Histogram of the birth density
+figure()
+subplot(111, xscale="linear", yscale="linear")
+hist(np.log10(stars_BirthDensity),density=True,bins=100)
+xlabel("${\\rm Birth Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("Probability", labelpad=-7)
+savefig("BirthDensity.png", dpi=200)
+
+# Plot of the specific star formation rate in the galaxy
+rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
+rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
+
+P_effective = KS_press_norm * ( rhos / EAGLEfloor_rho_norm)**(KS_gamma_effective)
+P_norm_high = KS_press_norm * (KS_high_den_thresh  / EAGLEfloor_rho_norm)**(KS_gamma_effective)
+sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
+KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
+sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**KS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+
+
+figure()
+subplot(111)
+hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
+plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
+plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
+xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=-7)
+savefig("density-sSFR.png", dpi=200)
+
+########################################################################3
+
 # Select gas in a pillow box around the galaxy
 mask = (
     (gas_pos[:, 0] > -15)
diff --git a/examples/IsolatedGalaxy_starformation_25p/KS_law.py b/examples/IsolatedGalaxy_starformation_25p/KS_law.py
deleted file mode 100644
index 296a3eed41..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/KS_law.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python
-import numpy as np
-import h5py
-import matplotlib.pyplot as plt
-
-# for the moment only use a single snapshot
-i = 29 
-
-# Load the snapshot
-#f = h5py.File('output_%04dcopy.hdf5'%i, 'r')
-f = h5py.File('output_%04d.hdf5'%i, 'r')
-
-# Load the stars
-starparticle = f['PartType4']
-
-# Load the properties of the stars
-coordinates = starparticle['Coordinates'][:,:]
-newstarflag = starparticle['NewStarFlag'][:]
-boxsize = f['Header'].attrs['BoxSize']
-masses_star = starparticle['Masses'][:]
-
-# Load the positions of the star particles and only take
-# the new star particles
-x = coordinates[:,0] - boxsize[0]/2.
-y = coordinates[:,1] - boxsize[1]/2.
-r = (x**2 + y**2)**.5
-birthdensity = starparticle['BirthDensity'][:]
-birthdensity=birthdensity[newstarflag==1]
-r_new = r[newstarflag==1]
-masses_new = masses_star[newstarflag==1]
-
-# Load the gas particles
-gasparticle = f['PartType0']
-
-# Load the properties of the gas particles
-coordinates_gas = gasparticle['Coordinates'][:,:]
-masses_gas = gasparticle['Masses'][:]
-
-# Load the positions of the gas particles
-x_gas = coordinates_gas[:,0] - boxsize[0]/2.
-y_gas = coordinates_gas[:,0] - boxsize[1]/2.
-r_gas = (x_gas**2 + y_gas**2)**.5
-
-r_array = np.linspace(0,15,100)
-area = np.zeros(len(r_array)-1)
-SigmaSFR = np.zeros(len(r_array)-1)
-Sigmagas = np.zeros(len(r_array)-1)
-massSFR_array = np.zeros(len(r_array)-1)
-massGAS_array = np.zeros(len(r_array)-1)
-for i in range(1,len(r_array)):
-    area_loop = np.pi*(r_array[i]**2 - r_array[i-1]**2)
-    mask = (r_new > r_array[i-1]) & (r_new < r_array[i])
-    mask2 = (r_gas > r_array[i-1]) & (r_gas < r_array[i])
-    massSFR = np.sum(masses_new[mask])
-    massSFR_array[i-1] = massSFR
-    massGAS = np.sum(masses_gas[mask2])
-    massGAS_array[i-1] = massGAS
-    #print('%d  area=%1.4f,  M_SFR=%1.4f,  M_gas=%1.4f'%(i,area_loop,massSFR,massGAS))
-    area[i-1] = area_loop
-    SigmaSFR[i-1] = massSFR/area_loop
-    Sigmagas[i-1] = massGAS/area_loop
-
-# Put SF law in correct units
-SigmaSFR *= 1e10/0.29e8
-Sigmagas *= 1e10/(1e3**2)
-SigmaSFR[SigmaSFR==0] = 1e-6
-print(Sigmagas)
-
-Sgas_array =np.logspace(-1,3,40) 
-Delta_Sgas = np.log10(Sgas_array[1]) - np.log10(Sgas_array[0])
-SigmaSFR_avg = np.zeros(len(Sgas_array))
-Sigmagas_avg = np.zeros(len(Sgas_array))
-
-for i in range(0,len(Sgas_array)):
-    totgasmass = 0.0
-    totnewstarmass = 0.0
-    area_loop = 0.0
-    for j in range(0,len(Sigmagas)):
-        #print('%d %1.4e %1.4e %1.4e %1.4e'%(j,Sigmagas[j],Sgas_array[i],10**(-Delta_Sgas),10**(Delta_Sgas)))
-        if (Sigmagas[j]>Sgas_array[i]*10**(-Delta_Sgas)) and (Sigmagas[j]<Sgas_array[i]*10**(Delta_Sgas)):
-            totgasmass += massGAS_array[j]
-            totnewstarmass += massSFR_array[j]
-            area_loop += area[j]
-    if area_loop==0.0:
-        SigmaSFR_avg[i] = 0.0
-        Sigmagas_avg[i] = 0.0
-    else:
-        SigmaSFR_avg[i] = totnewstarmass/area_loop
-        Sigmagas_avg[i] = totgasmass/area_loop
-
-print(SigmaSFR_avg)
-print(Sigmagas_avg)
-
-# Put SF law in correct units
-SigmaSFR_avg *= 1e10/0.29e8
-Sigmagas_avg *= 1e10/(1e6)
-SigmaSFR_avg[SigmaSFR_avg==0]=1e-6
-
-def KS_law(Sgas,n=1.4,A=1.515):
-    return A*1e-4*(Sgas)**n
-
-Sigmagaslaw = np.logspace(-1,3,100) 
-SigmaSFRlaw = KS_law(Sigmagaslaw)
-
-Z = 0
-ncrit = 10 # g/cm^3
-ncrit = 10 # M_sun /pc^2
-
-# plot SF law
-plt.loglog(Sigmagas,SigmaSFR,'o',label='Data')
-#plt.loglog(Sigmagas_avg,SigmaSFR_avg,'o',label='binned Data')
-plt.loglog(Sigmagaslaw, SigmaSFRlaw,'r-.',label='KS law')
-#plt.loglog(np.logspace(-1,3,40),1e-5*np.ones(40),'o')
-plt.axvline(x=10,linestyle='--', color='r',label='$\Sigma_c = 10$')
-plt.xlabel('$\Sigma_{gas}$')
-plt.ylabel('$\Sigma_{SFR}$')
-plt.legend()
-plt.show()
-   
-
-
-'''
-plt.hist(birthdensity,bins=100,range=[0,2.5])
-plt.show()
-
-plt.hist(x_gas,bins=100)
-plt.show()
-
-plt.hist(r,bins=100)
-plt.show()
-
-plt.hist(r_gas,bins=100)
-plt.show()
-'''
diff --git a/examples/IsolatedGalaxy_starformation_25p/StarFormationHistory.py b/examples/IsolatedGalaxy_starformation_25p/StarFormationHistory.py
deleted file mode 100644
index b08a7db365..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/StarFormationHistory.py
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/usr/bin/env python3
-import numpy as np
-import h5py
-import matplotlib.pyplot as plt
-
-snaps = 30 
-
-SFH = np.zeros(snaps)
-
-for i in range(0,snaps):
-    f = h5py.File('output_%04d.hdf5'%i, 'r')
-    header = f["Header"]
-    time = header.attrs['Time'][0]
-    print(time)
-    particles = f["PartType4"]
-    newstar = particles["NewStarFlag"][:]
-    SFH[i] = np.sum(newstar)
-
-plt.plot(SFH)
-plt.xlabel('Snapshot number')
-plt.ylabel('Total number of stars formed')
-plt.show()
-
-
-
-
-f = h5py.File('output_%04d.hdf5'%i, 'r')
-particles = f["PartType4"]
-newstar = particles["NewStarFlag"][:]
-Coordinates = particles["Coordinates"][:,:]
-x = Coordinates[:,0][newstar==1]
-y = Coordinates[:,1][newstar==1]
-
-print(len(x), len(y))
-plt.plot(x,y,'.')
-plt.show()
diff --git a/examples/IsolatedGalaxy_starformation_25p/gasdensity.py b/examples/IsolatedGalaxy_starformation_25p/gasdensity.py
deleted file mode 100644
index 7b7d7e190a..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/gasdensity.py
+++ /dev/null
@@ -1,563 +0,0 @@
-"""
-Makes a movie using sphviewer and ffmpeg.
-
-Edit low_frac and up_frac to focus on a certain view of the box.
-The colour map can also be changed via colour_map.
-
-Usage: python3 makeMovie.py CoolingHalo_
-
-Written by James Willis (james.s.willis@durham.ac.uk)
-"""
-
-import glob
-import h5py as h5
-import numpy as np
-import matplotlib.pyplot as plt
-
-from multiprocessing import Pool
-from tqdm import tqdm
-from pylab import sys
-from sphviewer.tools import QuickView
-from matplotlib.animation import FuncAnimation
-from matplotlib.colors import LogNorm
-from scipy.optimize import curve_fit
-
-def gen_data(filename,pixels,case):
-    """
-    Data generation function. Reads a snapshot, filters the particle data and
-    creates an image using QuickView from sphviewer.
-    """
-
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates = f["/PartType0/Coordinates"][:,:]
-        masses = f["/PartType0/Masses"][:]
-        smoothing_lengths = f["/PartType0/SmoothingLength"][:]
-        coordinates_star= f["/PartType4/Coordinates"][:,:]
-        star_flag = f["/PartType4/NewStarFlag"][:]
-        SFR = f["/PartType0/SFR"][:]
-    
-    # Filter data  
-    low_frac=0.48
-    up_frac=0.52
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz)
-
-    sfr_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (SFR > 0 )
-    
-    masses = masses[part_mask]
-    SFR = SFR[sfr_mask]
-   
-    if case ==0:
-        coordinates = coordinates[part_mask]
-        smoothing_lengths = smoothing_lengths[part_mask]
-        masvar = masses
-    elif case==1:
-        coordinates = coordinates[part_mask]
-        smoothing_lengths = smoothing_lengths[part_mask]
-        masvar = masses**2
-    elif case==2:
-        coordinates = coordinates[sfr_mask]
-        smoothing_lengths = smoothing_lengths[sfr_mask]
-        masvar = SFR
-    # Use sphviewer QuickView to generate an image
-    qv = QuickView(
-        coordinates,
-        mass=masvar,
-        r="infinity",
-        xsize=pixels,
-        ysize=pixels,
-        p=0,  # Viewing angle theta
-        roll=0,  # Viewing angle phi
-        plot=False,
-        logscale=False,
-        hsml=smoothing_lengths
-    )
-
-    image_data = qv.get_image()
-    extent = qv.get_extent()
-   
-    return image_data, extent
-    
-def getnewstars(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        star_flag = f["/PartType4/NewStarFlag"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (star_flag == 1)
-
-    return coordinates[:,0][part_mask]-box_size/2., coordinates[:,1][part_mask]-box_size/2.
-
-def getdensfrtemp(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        SFR = f["/PartType0/SFR"][:]
-        density = f["/PartType0/Density"][:]
-        temperature = f["/PartType0/Temperature"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) 
-    SFR = SFR[part_mask]
-    density = density[part_mask]
-    temperature = temperature[part_mask]
-
-    return density, SFR, temperature
-
-def getbirthdensity(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        birthdensity = f["/PartType4/BirthDensity"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (birthdensity>0)
-
-    birthdensity = birthdensity[part_mask]
-
-    return birthdensity 
-    
-def getdensityrate(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        sfrrate = f["/PartType0/sSFR"][:]
-        density = f["/PartType0/Density"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (sfrrate>0)
-
-    sfrrate = sfrrate[part_mask]
-    density = density[part_mask]
-
-    return sfrrate, density 
-
-def getgasmass(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        mass = f["/PartType0/Masses"][:][0]
-
-    return mass*1e10
-
-def getSFH(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        mass = f["/PartType4/Masses"][:]
-        flag = f["/PartType4/NewStarFlag"][:]
-        birth_time = f["/PartType4/Birth_time"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (flag==1)
-
-    birth_time = birth_time[part_mask]
-    mass = mass[part_mask]
-    
-    histogram = np.histogram(birth_time,bins=200,weights=mass*2e4,range=[0,.1])
-    values = histogram[0]
-    xvalues =( histogram[1][:-1] +histogram[1][1:])/2.
-    return xvalues, values
-
-def KS_law(Sigmagas,n=1.4,A=1.515e-4):
-    return A * Sigmagas**n
-
-def specificSFR(rho,fg=.25,gamma_eos=4/3):
-    return .363 * (10 * fg * (rho/0.1)**gamma_eos )**.2 
-
-def specificSFRhighden(rho,gamma_eos=4/3):
-    return specificSFR(1e3)* (( rho/1e3 )**gamma_eos)**.5
-
-def linearfunc(x,a,b):
-    return a*x + b
-    
-def makefullplot(filename,binsize=200):
-    # Get the Surface density from the snapshot
-    image_data,extent = gen_data(filename,binsize,case=0)
-    image_data2 = image_data[:,:]*1e10/(1e2**2)
-    surface_density = image_data[:,:]*1e10/(1e2**2)
-    if binsize==20:
-        image_data2 /= 1e2
-        surface_density /= 1e2
-    minvalue = np.min(image_data2[image_data2!=0])
-    image_data2[image_data2==0] = minvalue*.9
-    
-    # Get the SFR surface density from the snapshot
-    image_data,extent = gen_data(filename,binsize,case=2)
-    image_data3 = image_data[:,:] * 1e10/1e9/(.1**2) 
-    SFR_smoothed = image_data3/image_data2
-    SFR_surface_density = image_data[:,:] * 1e10/1e9/(.1**2)/surface_density
-    if binsize==20:
-        image_data3 /= 1e2
-        SFR_surface_density /= 1e2
-        SFR_smoothed /= 1e2
-    plot_SFR_smoothed = np.ones(np.shape(image_data3))*np.min(SFR_smoothed[SFR_smoothed!=0])*.9
-    plot_SFR_smoothed[image_data3!=0] = SFR_smoothed[image_data3!=0]
-    
-    # Define the font used in the plots 
-    font = {'color':'white', 'size':12}
-
-    # Make the first plot of the Gas density
-    fig = plt.figure(1,figsize=(16.5,16.5))
-    plt.subplot(3,3,1)
-    plt.imshow(np.log10(image_data2), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,      # ticks along the bottom edge are off
-        right=False,         # ticks along the top edge are off
-        labelleft=False) # labels along the bottom edge are off
-    cbar = plt.colorbar() 
-    plt.clim(-1.5,3.0)
-    #cbar.set_label('Surface density ($M_\odot pc^{-2}$)', rotation=270)
-    
-    # Make the second plot of the SFR density
-    plt.subplot(3,3,2)
-    plt.imshow(np.log10(plot_SFR_smoothed), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,      # ticks along the bottom edge are off
-        right=False,         # ticks along the top edge are off
-        labelleft=False) # labels along the bottom edge are off
-    cbar = plt.colorbar() 
-    plt.clim(-5,-1)
-
-    # Make the third plot of the stars that are formed
-    plt.subplot(3,3,3)
-    xnewstar, ynewstar = getnewstars(filename)
-    newstarmatrix = np.histogram2d(xnewstar,ynewstar,bins=binsize)[0]
-    newstarmatrix[newstarmatrix==0] = 1e-1
-    plt.imshow(np.log10(newstarmatrix), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('New stars')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,        # ticks along the bottom edge are off
-        right=False,       # ticks along the top edge are off
-        labelleft=False)   # labels along the bottom edge are off
-    cbar = plt.colorbar()
-    plt.clim(-1,3)
-    
-    # Make the fourth plot of the KS law
-    plt.subplot(3,3,4)
-    sigma_gas = image_data2.flatten()
-    sigma_star = SFR_smoothed.flatten()
-    sigma_gas=sigma_gas[sigma_star>0]
-    sigma_star=sigma_star[sigma_star>0]
-    sigma_gas2 = surface_density.flatten()
-    sigma_star2 = SFR_surface_density.flatten()
-
-    #KShistogram = np.histogram2d(np.log10(sigma_gas),np.log10(sigma_star),bins=50,range=[[-1,4],[-5,2]])
-    #plt.imshow(np.transpose(KShistogram[0]), extent=np.array([-1,4,-5,2]))
-    plt.title('KS law')
-
-    fitresult = curve_fit(linearfunc,np.log10(sigma_gas),np.log10(sigma_star),p0=(1.4,np.log10(1.515e-4)),bounds=((0,-6),(3,-2))) 
-    #print(fitresult[0], fitresult[1])
-    
-    sigma_gas_range = 10**np.linspace(-1,2,100)
-    sigma_star_range = KS_law(sigma_gas_range)
-    sigma_star_one_range = KS_law(sigma_gas_range,n=0.7)
-
-    if binsize==200:
-        plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=50)
-    else:
-        plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=25)
-    #plt.hist2d(np.log10(sigma_gas2),np.log10(sigma_star2),range=[[-1,2],[-5,-2]],bins=50)
-    plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_range),'k--',label='EAGLE SF')
-    plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_one_range),'k-.',label='n=.7')
-    plt.plot(np.log10(sigma_gas_range), np.log10(KS_law(sigma_gas_range,n=0.6)),'k:',label='n=.6')
-    #plt.plot(np.log10(sigma_gas_range), linearfunc(np.log10(sigma_gas_range), *fitresult[0]),'r--',label='fit')
-    plt.xlabel('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
-    plt.ylabel('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
-    cbar = plt.colorbar()
-    if binsize==200:
-        plt.clim(0,700)
-    elif binsize==20:
-        plt.clim(0,40)
-    #plt.scatter(np.log10(sigma_gas),np.log10(sigma_star))
-    plt.legend()
-
-    # make a plot of the gas density vs SFR
-    plt.subplot(3,3,5)
-    critden = .1*((0.0129)/.002)**(-.64)
-    mu = 1.3
-    density, SFR, temperature = getdensfrtemp(filename)
-    density *= 1e10/(1e3)**3*40.4759/mu
-    SFR *= 1e10/1e9
-    SFR_plot = SFR[:]
-    plt.hist2d(np.log10(density[SFR_plot>0]),np.log10(SFR_plot[SFR_plot>0]), bins=50,range=[[np.log10(critden),5],[-4.4,-2.5]])
-    #plt.loglog(density*1e10/(1e3)**3*40.4759/mu,SFR_plot*1e10/1e9,'.')
-    xx = 10**np.linspace(-1,2,100)
-    def sfr_model(density,critden=.1,gammaeff=4./3.,mgas=1.4995e4):
-        return 6.9e-11 * mgas * (density/critden)**(gammaeff/5.)
-    rhos = 10**np.linspace(-1,4,100)
-    rhoshigh = 10**np.linspace(2,5,100)
-    massparticle =getgasmass(filename) 
-    #plt.plot(np.log10(xx),np.log10(1e-4*xx**(4./15.)),'k--',label='Slope of EAGLE SF')
-    #plt.plot(np.log10(xx),np.log10(sfr_model(xx)),'r--',label='Slope of EAGLE SF')
-    plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)*massparticle/1e9),'k--',label='SFR low density EAGLE')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.ylabel('$\log$ Star Formation rate (SFR) [$ M_\odot \\rm yr^{-1} $]')
-    plt.xlim(np.log10(critden),5)
-    cbar = plt.colorbar()
-    plt.legend()
-
-    SFR_tracer = np.log10(SFR_plot[SFR_plot>0])
-    
-    plt.subplot(3,3,6)
-    plt.scatter(np.log10(density[SFR_plot>0]),np.log10(temperature[SFR_plot>0]),s=1,c=SFR_tracer,cmap='viridis')
-    plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
-    plt.axvline(x=0,linestyle='--',color='k',linewidth=3,label='End of Wiersma Tables')
-    cbar = plt.colorbar()
-    plt.clim(-4.4,-2.5)
-    plt.scatter(np.log10(density[SFR_plot<=0]),np.log10(temperature[SFR_plot<=0]),s=1,c='gray')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.ylabel('$\log$ temperature [$\\rm K$]')
-    plt.ylim(2,4)
-    plt.xlim(-3.5,5)
-    plt.legend()
-
-    plt.subplot(3,3,7)
-    birthdensity = getbirthdensity(filename)*404.759
-    plt.hist(np.log10(birthdensity),bins=50,density=True,range=(np.log10(critden),5))
-    plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.xlim(np.log10(critden),5)
-    plt.legend()
-
-    plt.subplot(3,3,8)
-    sfrrate, gasdensity = getdensityrate(filename)
-    gasdensity *= 1e10/(1e3)**3*40.4759/mu
-    plt.hist2d(np.log10(gasdensity),np.log10(sfrrate), bins=50,range=[[-1.5,5],[-.5,2.5]])
-    cbar = plt.colorbar()
-    rhos = 10**np.linspace(-1,4,100)
-    rhoshigh = 10**np.linspace(2,5,100)
-    plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)),'k--',label='sSFR low density EAGLE')
-    plt.plot(np.log10(rhoshigh),np.log10(specificSFRhighden(rhoshigh)),'k--',label='sSFR high density EAGLE')
-    plt.xlabel('density [\\rm cm^{-3}]')
-    plt.ylabel('sSFR ($\\rm Gyr^{-1}$)')
-    plt.legend()
-    
-    plt.subplot(3,3,9)
-    timearray, SFH = getSFH(filename)
-    plt.plot(timearray*1e3,SFH, label='spart birth_time')
-    plt.xlabel('Time (Myr)')
-    plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
-    plt.legend()
-
-    if binsize==200:
-        plt.savefig('./ksplot/'+filename+'_0.1kpc.png')
-    elif binsize==20:
-        plt.savefig('./ksplot_1kpc/'+filename+'_0.1kpc.png')
-    plt.close()
-
-
-def getmassandsfr(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        SFR = f["/PartType0/SFR"][:]
-        coordinates_star = f["/PartType4/Coordinates"][:,:]
-        masses_star = f["/PartType4/Masses"][:]
-        
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (SFR>0) 
-    mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) < 
-        absmaxz)  
-
-    SFR_snap = np.sum(SFR[part_mask])
-    mass_star = np.sum(masses_star[mask])
-    return SFR_snap, mass_star
-
-def getsfrsnapwide():
-    
-    time = np.arange(1,101,1)
-    SFR_sparticles = np.zeros(100)
-    SFR_gparticles = np.zeros(100)
-    new_sparticles = np.zeros(100)
-    previous_mass = 0
-    previous_numb = 0
-    for i in range(1,43):
-        print(i)
-        # Read the data
-        filename= 'output_%04d.hdf5'%i
-        with h5.File(filename, "r") as f:
-            box_size = f["/Header"].attrs["BoxSize"][0]
-            coordinates= f["/PartType0/Coordinates"][:,:]
-            SFR = f["/PartType0/SFR"][:]
-            coordinates_star = f["/PartType4/Coordinates"][:,:]
-            masses_star = f["/PartType4/Masses"][:]
-            newstar = f["/PartType4/NewStarFlag"][:]
-
-        absmaxz = 2 #kpc
-        absmaxxy = 10 #kpc
-
-        part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-            absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-            absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-            absmaxz) & (SFR>0) 
-        mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) < 
-            absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) < 
-            absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) < 
-            absmaxz) & (newstar==1) 
-
-        SFR = SFR[part_mask]
-        masses_star = masses_star[mask]*1e4
-        newstar = newstar[mask]
-
-        total_SFR = np.sum(SFR)
-        total_mass = np.sum(masses_star)
-        new_spart = np.sum(newstar)
-        new_sparticles[i] = new_spart - previous_numb
-        SFR_sparticles[i] = total_mass - previous_mass 
-        SFR_gparticles[i] = total_SFR*10
-        previous_numb = new_spart
-        previous_mass = total_mass
-
-    factor = masses_star[0]
-    print(new_sparticles)
-    print(len(time),len(SFR_sparticles),len(SFR_gparticles))
-    return time[:-1], SFR_sparticles[1:], SFR_gparticles[1:], np.sqrt(new_sparticles[1:])*factor
-
-
-
-if __name__ == "__main__":
-
-    # Some global variables to define the whole run
-    input_file_name = sys.argv[1]
-    pixels = 1000 
-    num_procs = 4 # No. of processors to use
-    colour_map = "inferno"
-    
-    # Find the list of files and sort them in numerical order
-    file_list = sorted(glob.glob(input_file_name + "*.hdf5"), key=lambda x:
-            int(x[len(input_file_name):len(input_file_name)+4]))
-    
-    # get the image data 
-    #for i in [1,5,10,33,66,90]:
-    #    print('output_%04d.hdf5'%i)
-    #    makefullplot('output_%04d.hdf5'%i)
-    #    makefullplot('output_%04d.hdf5'%i,binsize=20)
-
-    time, SFR1, SFR2, SFR_error = getsfrsnapwide()
-    time2, SFR3 = getSFH('output_%04d.hdf5'%42)
-    plt.plot(time,SFR1,label='Using created star particles')
-    print(SFR_error)
-    plt.fill_between(time, SFR1+SFR_error,SFR1-SFR_error,label='Uncertainty',alpha=.3)
-    plt.plot(time,SFR2,label='Using SFR of gas particles')
-    plt.plot(time2*1e3,SFR3,label='Using birth_time of star particles')
-    plt.xlabel('Time (Myr)')
-    plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
-    plt.legend()
-    plt.show()
-    '''
-    snapshots = 52 
-    mass = np.zeros(snapshots)
-    sfr = np.zeros(snapshots)
-    for i in range(1,snapshots):
-        if i!=6:
-            print(i)
-            sfr[i], mass[i] = getmassandsfr('output_%04d.hdf5'%i)
-        
-    plt.plot(np.log10(mass*1e10))
-    plt.xlabel('Snapshot number (Myr)')
-    plt.ylabel('Mass')
-    plt.savefig('Mass increase')
-    plt.close()
-
-    plt.plot(np.log10(sfr*10))
-    plt.xlabel('Snapshot number (Myr)')
-    plt.ylabel('SFR ($M_\odot / \\rm yr$)')
-    plt.savefig('SFR')
-    plt.close()
-    '''
-
-    for i in [1,5,10,33,66,90]:
-        if i!=6:
-            print('output_%04d.hdf5'%i)
-            makefullplot('output_%04d.hdf5'%i)
-            makefullplot('output_%04d.hdf5'%i,binsize=20)
-
-
diff --git a/examples/IsolatedGalaxy_starformation_25p/getIC.sh b/examples/IsolatedGalaxy_starformation_25p/getIC.sh
deleted file mode 100755
index 25cee8dfc0..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/getIC.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-wget https://home.strw.leidenuniv.nl/~nobels/data/starforming_25p.hdf5
diff --git a/examples/IsolatedGalaxy_starformation_25p/isolated_galaxy.yml b/examples/IsolatedGalaxy_starformation_25p/isolated_galaxy.yml
deleted file mode 100644
index 9082855933..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/isolated_galaxy.yml
+++ /dev/null
@@ -1,115 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
-  UnitMass_in_cgs:     1.9891E43   # 10^10 solar masses 
-  UnitLength_in_cgs:   3.08567758E21   # 1 kpc 
-  UnitVelocity_in_cgs: 1E5   # km/s
-  UnitCurrent_in_cgs:  1   # Amperes
-  UnitTemp_in_cgs:     1   # Kelvin
-
-# Parameters for the self-gravity scheme
-Gravity:
-  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
-  eta:          0.025               # Constant dimensionless multiplier for time integration.
-  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.01 # Comoving softening length (in internal units).
-  max_physical_softening: 0.01    # Physical softening length (in internal units).
-
-# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
-TimeIntegration:
-  time_begin:        0.    # The starting time of the simulation (in internal units).
-  time_end:          2.    # The end time of the simulation (in internal units).
-  dt_min:            1e-9  # The minimal time-step size of the simulation (in internal units).
-  dt_max:            1e-2  # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
-  basename:   output      # Common part of the name of output files
-  time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
-  delta_time: 0.001        # Time difference between consecutive outputs (in internal units)
-
-  
-# Parameters governing the conserved quantities statistics
-Statistics:
-  delta_time:           1e-2     # Time between statistics output
-  time_first:             0.     # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
-
-# Parameters related to the initial conditions
-InitialConditions:
-  file_name:  starforming_25p.hdf5 # The file to read
-  periodic:                    0    # Are we running with periodic ICs?
-
-# Parameters for the hydrodynamics scheme
-SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  h_tolerance:           1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
-  h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
-  max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
-  max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
-  minimal_temperature:   100        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
-  H_mass_fraction:       0.755    # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
-  H_ionization_temperature: 1e4   # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
-
-EAGLECooling:
-  dir_name:                ./coolingtables/  # Location of the Wiersma+08 cooling tables
-  H_reion_z:               11.5              # Redshift of Hydrogen re-ionization
-  He_reion_z_centre:       3.5               # Redshift of the centre of the Helium re-ionization Gaussian
-  He_reion_z_sigma:        0.5               # Spread in redshift of the  Helium re-ionization Gaussian
-  He_reion_eV_p_H:         2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
-
-# Primordial abundances
-EAGLEChemistry:
-  init_abundance_metal:     0.0129          # Inital fraction of particle mass in *all* metals
-  init_abundance_Hydrogen:  0.7065       # Inital fraction of particle mass in Hydrogen
-  init_abundance_Helium:    0.2806        # Inital fraction of particle mass in Helium
-  init_abundance_Carbon:    0.00207        # Inital fraction of particle mass in Carbon
-  init_abundance_Nitrogen:  0.000836        # Inital fraction of particle mass in Nitrogen
-  init_abundance_Oxygen:    0.00549        # Inital fraction of particle mass in Oxygen
-  init_abundance_Neon:      0.00141        # Inital fraction of particle mass in Neon
-  init_abundance_Magnesium: 0.000591        # Inital fraction of particle mass in Magnesium
-  init_abundance_Silicon:   0.000683        # Inital fraction of particle mass in Silicon
-  init_abundance_Iron:      0.0011        # Inital fraction of particle mass in Iron
-
-# Hernquist potential parameters
-HernquistPotential:
-  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
-  position:        [0.,0.,0.]    # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
-  idealizeddisk:   1        # Run with an idealized galaxy disk
-  M200:            30.0   # M200 of the galaxy disk
-  h:               0.704    # reduced Hubble constant (value does not specify the used units!)
-  concentration:   7.1      # concentration of the Halo
-  diskfraction:              0.040   # Disk mass fraction
-  bulgefraction:              0.010   # Bulge mass fraction
-  timestep_mult:   0.01     # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
-  epsilon:         0.1      # Softening size (internal units)
- 
-# Schaye and Dalla Vecchia 2008 star formation
-SchayeSF:
-  thresh_MinOverDens:               57.7      # The critical density contrast to form stars
-  thresh_temp:                      1e5     # The critical temperature to form stars
-  fg:                               0.25    # The mass fraction in gas
-  SchmidtLawCoeff_MSUNpYRpKPC2:     1.515e-4  # The normalization of the Kennicutt-Schmidt law
-  SchmidtLawExponent:               1.4     # The power law of the Kennicutt-Schmidt law
-  SchmidtLawHighDensExponent:       2.0     # The high density exponent for the Kennicutt-Schmidt law
-  SchmidtLawHighDens_thresh_HpCM3:  1e3
-  Schaye2004:                       1       # whether to use the metallicity dependent critical star formation of Schaye (2004) (1) or not (0).
-  thresh_norm_HpCM3:                0.1      # Critical sf normalization to use (is not a normalization when Schaye2004=0, than it is the value.
-  thresh_max_norm_HpCM3:            10.0    # Maximum norm of the critical sf density
-  MetDep_Z0:                        0.002   # Scale metallicity to use for the equation (not used for Schaye2004=0)
-  MetDep_SFthresh_Slope:            -0.64   # Scaling of the critical density with the metallicity (not used for Schaye2004=0)
-  thresh_MaxPhysDensOn:             0       # Default is 0.
-  thresh_MaxOverDens_HpCM3:         1e5     # Density at which the SF law changes
-  EOS_Jeans_GammaEffective:         1.33333 # The polytropic index
-  EOS_Jeans_TemperatureNorm_K:      1e5     # No idea how this works
-  EOS_JEANS_DensityNorm_HpCM3:      0.1      # No idea what the value is.
-
-# Parameters for the EAGLE "equation of state"
-EAGLEEntropyFloor:
-  Jeans_density_threshold_H_p_cm3: 0.1       # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
-  Jeans_over_density_threshold:    10.       # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
-  Jeans_temperature_norm_K:        8000      # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
-  Jeans_gamma_effective:           1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
-  Cool_density_threshold_H_p_cm3: 1e-5       # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
-  Cool_over_density_threshold:    10.        # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
-  Cool_temperature_norm_K:        8000       # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
-  Cool_gamma_effective:           1.         # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/IsolatedGalaxy_starformation_25p/run.sh b/examples/IsolatedGalaxy_starformation_25p/run.sh
deleted file mode 100755
index d20476b86f..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p/run.sh
+++ /dev/null
@@ -1,3 +0,0 @@
-#!/bin/bash
-
-../swift --threads=32 --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
diff --git a/examples/IsolatedGalaxy_starformation_25p_intermediate/gasdensity.py b/examples/IsolatedGalaxy_starformation_25p_intermediate/gasdensity.py
deleted file mode 100644
index f9c5dce2cf..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p_intermediate/gasdensity.py
+++ /dev/null
@@ -1,563 +0,0 @@
-"""
-Makes a movie using sphviewer and ffmpeg.
-
-Edit low_frac and up_frac to focus on a certain view of the box.
-The colour map can also be changed via colour_map.
-
-Usage: python3 makeMovie.py CoolingHalo_
-
-Written by James Willis (james.s.willis@durham.ac.uk)
-"""
-
-import glob
-import h5py as h5
-import numpy as np
-import matplotlib.pyplot as plt
-
-from multiprocessing import Pool
-from tqdm import tqdm
-from pylab import sys
-from sphviewer.tools import QuickView
-from matplotlib.animation import FuncAnimation
-from matplotlib.colors import LogNorm
-from scipy.optimize import curve_fit
-
-def gen_data(filename,pixels,case):
-    """
-    Data generation function. Reads a snapshot, filters the particle data and
-    creates an image using QuickView from sphviewer.
-    """
-
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates = f["/PartType0/Coordinates"][:,:]
-        masses = f["/PartType0/Masses"][:]
-        smoothing_lengths = f["/PartType0/SmoothingLength"][:]
-        coordinates_star= f["/PartType4/Coordinates"][:,:]
-        star_flag = f["/PartType4/NewStarFlag"][:]
-        SFR = f["/PartType0/SFR"][:]
-    
-    # Filter data  
-    low_frac=0.48
-    up_frac=0.52
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz)
-
-    sfr_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (SFR > 0 )
-    
-    masses = masses[part_mask]
-    SFR = SFR[sfr_mask]
-   
-    if case ==0:
-        coordinates = coordinates[part_mask]
-        smoothing_lengths = smoothing_lengths[part_mask]
-        masvar = masses
-    elif case==1:
-        coordinates = coordinates[part_mask]
-        smoothing_lengths = smoothing_lengths[part_mask]
-        masvar = masses**2
-    elif case==2:
-        coordinates = coordinates[sfr_mask]
-        smoothing_lengths = smoothing_lengths[sfr_mask]
-        masvar = SFR
-    # Use sphviewer QuickView to generate an image
-    qv = QuickView(
-        coordinates,
-        mass=masvar,
-        r="infinity",
-        xsize=pixels,
-        ysize=pixels,
-        p=0,  # Viewing angle theta
-        roll=0,  # Viewing angle phi
-        plot=False,
-        logscale=False,
-        hsml=smoothing_lengths
-    )
-
-    image_data = qv.get_image()
-    extent = qv.get_extent()
-   
-    return image_data, extent
-    
-def getnewstars(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        star_flag = f["/PartType4/NewStarFlag"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (star_flag == 1)
-
-    return coordinates[:,0][part_mask]-box_size/2., coordinates[:,1][part_mask]-box_size/2.
-
-def getdensfrtemp(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        SFR = f["/PartType0/SFR"][:]
-        density = f["/PartType0/Density"][:]
-        temperature = f["/PartType0/Temperature"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) 
-    SFR = SFR[part_mask]
-    density = density[part_mask]
-    temperature = temperature[part_mask]
-
-    return density, SFR, temperature
-
-def getbirthdensity(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        birthdensity = f["/PartType4/BirthDensity"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (birthdensity>0)
-
-    birthdensity = birthdensity[part_mask]
-
-    return birthdensity 
-    
-def getdensityrate(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        sfrrate = f["/PartType0/sSFR"][:]
-        density = f["/PartType0/Density"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (sfrrate>0)
-
-    sfrrate = sfrrate[part_mask]
-    density = density[part_mask]
-
-    return sfrrate, density 
-
-def getgasmass(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        mass = f["/PartType0/Masses"][:][0]
-
-    return mass*1e10
-
-def getSFH(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType4/Coordinates"][:,:]
-        mass = f["/PartType4/Masses"][:]
-        flag = f["/PartType4/NewStarFlag"][:]
-        birth_time = f["/PartType4/Birth_time"][:]
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (flag==1)
-
-    birth_time = birth_time[part_mask]
-    mass = mass[part_mask]
-    
-    histogram = np.histogram(birth_time,bins=200,weights=mass*2e4,range=[0,.1])
-    values = histogram[0]
-    xvalues =( histogram[1][:-1] +histogram[1][1:])/2.
-    return xvalues, values
-
-def KS_law(Sigmagas,n=1.4,A=1.515e-4):
-    return A * Sigmagas**n
-
-def specificSFR(rho,fg=.25,gamma_eos=4/3):
-    return .363 * (10 * fg * (rho/0.1)**gamma_eos )**.2 
-
-def specificSFRhighden(rho,gamma_eos=4/3):
-    return specificSFR(1e3)* (( rho/1e3 )**gamma_eos)**.5
-
-def linearfunc(x,a,b):
-    return a*x + b
-    
-def makefullplot(filename,binsize=200):
-    # Get the Surface density from the snapshot
-    image_data,extent = gen_data(filename,binsize,case=0)
-    image_data2 = image_data[:,:]*1e10/(1e2**2)
-    surface_density = image_data[:,:]*1e10/(1e2**2)
-    if binsize==20:
-        image_data2 /= 1e2
-        surface_density /= 1e2
-    minvalue = np.min(image_data2[image_data2!=0])
-    image_data2[image_data2==0] = minvalue*.9
-    
-    # Get the SFR surface density from the snapshot
-    image_data,extent = gen_data(filename,binsize,case=2)
-    image_data3 = image_data[:,:] * 1e10/1e9/(.1**2) 
-    SFR_smoothed = image_data3/image_data2
-    SFR_surface_density = image_data[:,:] * 1e10/1e9/(.1**2)/surface_density
-    if binsize==20:
-        image_data3 /= 1e2
-        SFR_surface_density /= 1e2
-        SFR_smoothed /= 1e2
-    plot_SFR_smoothed = np.ones(np.shape(image_data3))*np.min(SFR_smoothed[SFR_smoothed!=0])*.9
-    plot_SFR_smoothed[image_data3!=0] = SFR_smoothed[image_data3!=0]
-    
-    # Define the font used in the plots 
-    font = {'color':'white', 'size':12}
-
-    # Make the first plot of the Gas density
-    fig = plt.figure(1,figsize=(16.5,16.5))
-    plt.subplot(3,3,1)
-    plt.imshow(np.log10(image_data2), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,      # ticks along the bottom edge are off
-        right=False,         # ticks along the top edge are off
-        labelleft=False) # labels along the bottom edge are off
-    cbar = plt.colorbar() 
-    plt.clim(-1.5,3.0)
-    #cbar.set_label('Surface density ($M_\odot pc^{-2}$)', rotation=270)
-    
-    # Make the second plot of the SFR density
-    plt.subplot(3,3,2)
-    plt.imshow(np.log10(plot_SFR_smoothed), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,      # ticks along the bottom edge are off
-        right=False,         # ticks along the top edge are off
-        labelleft=False) # labels along the bottom edge are off
-    cbar = plt.colorbar() 
-    plt.clim(-5,-1)
-
-    # Make the third plot of the stars that are formed
-    plt.subplot(3,3,3)
-    xnewstar, ynewstar = getnewstars(filename)
-    newstarmatrix = np.histogram2d(xnewstar,ynewstar,bins=binsize)[0]
-    newstarmatrix[newstarmatrix==0] = 1e-1
-    plt.imshow(np.log10(newstarmatrix), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
-    plt.title('New stars')
-    plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
-    plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
-    plt.tick_params(
-        axis='x',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        bottom=False,      # ticks along the bottom edge are off
-        top=False,         # ticks along the top edge are off
-        labelbottom=False) # labels along the bottom edge are off
-    plt.tick_params(
-        axis='y',          # changes apply to the x-axis
-        which='both',      # both major and minor ticks are affected
-        left=False,        # ticks along the bottom edge are off
-        right=False,       # ticks along the top edge are off
-        labelleft=False)   # labels along the bottom edge are off
-    cbar = plt.colorbar()
-    plt.clim(-1,3)
-    
-    # Make the fourth plot of the KS law
-    plt.subplot(3,3,4)
-    sigma_gas = image_data2.flatten()
-    sigma_star = SFR_smoothed.flatten()
-    sigma_gas=sigma_gas[sigma_star>0]
-    sigma_star=sigma_star[sigma_star>0]
-    sigma_gas2 = surface_density.flatten()
-    sigma_star2 = SFR_surface_density.flatten()
-
-    #KShistogram = np.histogram2d(np.log10(sigma_gas),np.log10(sigma_star),bins=50,range=[[-1,4],[-5,2]])
-    #plt.imshow(np.transpose(KShistogram[0]), extent=np.array([-1,4,-5,2]))
-    plt.title('KS law')
-
-    fitresult = curve_fit(linearfunc,np.log10(sigma_gas),np.log10(sigma_star),p0=(1.4,np.log10(1.515e-4)),bounds=((0,-6),(3,-2))) 
-    #print(fitresult[0], fitresult[1])
-    
-    sigma_gas_range = 10**np.linspace(-1,2,100)
-    sigma_star_range = KS_law(sigma_gas_range)
-    sigma_star_one_range = KS_law(sigma_gas_range,n=0.7)
-
-    if binsize==200:
-        plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=50)
-    else:
-        plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=25)
-    #plt.hist2d(np.log10(sigma_gas2),np.log10(sigma_star2),range=[[-1,2],[-5,-2]],bins=50)
-    plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_range),'k--',label='EAGLE SF')
-    plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_one_range),'k-.',label='n=.7')
-    plt.plot(np.log10(sigma_gas_range), np.log10(KS_law(sigma_gas_range,n=0.6)),'k:',label='n=.6')
-    #plt.plot(np.log10(sigma_gas_range), linearfunc(np.log10(sigma_gas_range), *fitresult[0]),'r--',label='fit')
-    plt.xlabel('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
-    plt.ylabel('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
-    cbar = plt.colorbar()
-    if binsize==200:
-        plt.clim(0,700)
-    elif binsize==20:
-        plt.clim(0,40)
-    #plt.scatter(np.log10(sigma_gas),np.log10(sigma_star))
-    plt.legend()
-
-    # make a plot of the gas density vs SFR
-    plt.subplot(3,3,5)
-    critden = .1*((0.0129)/.002)**(-.64)
-    mu = 1.3
-    density, SFR, temperature = getdensfrtemp(filename)
-    density *= 1e10/(1e3)**3*40.4759/mu
-    SFR *= 1e10/1e9
-    SFR_plot = SFR[:]
-    plt.hist2d(np.log10(density[SFR_plot>0]),np.log10(SFR_plot[SFR_plot>0]), bins=50,range=[[np.log10(critden),5],[-4.4,-2.5]])
-    #plt.loglog(density*1e10/(1e3)**3*40.4759/mu,SFR_plot*1e10/1e9,'.')
-    xx = 10**np.linspace(-1,2,100)
-    def sfr_model(density,critden=.1,gammaeff=4./3.,mgas=1.4995e4):
-        return 6.9e-11 * mgas * (density/critden)**(gammaeff/5.)
-    rhos = 10**np.linspace(-1,4,100)
-    rhoshigh = 10**np.linspace(2,5,100)
-    massparticle =getgasmass(filename) 
-    #plt.plot(np.log10(xx),np.log10(1e-4*xx**(4./15.)),'k--',label='Slope of EAGLE SF')
-    #plt.plot(np.log10(xx),np.log10(sfr_model(xx)),'r--',label='Slope of EAGLE SF')
-    plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)*massparticle/1e9),'k--',label='SFR low density EAGLE')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.ylabel('$\log$ Star Formation rate (SFR) [$ M_\odot \\rm yr^{-1} $]')
-    plt.xlim(np.log10(critden),5)
-    cbar = plt.colorbar()
-    plt.legend()
-
-    SFR_tracer = np.log10(SFR_plot[SFR_plot>0])
-    
-    plt.subplot(3,3,6)
-    plt.scatter(np.log10(density[SFR_plot>0]),np.log10(temperature[SFR_plot>0]),s=1,c=SFR_tracer,cmap='viridis')
-    plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
-    plt.axvline(x=0,linestyle='--',color='k',linewidth=3,label='End of Wiersma Tables')
-    cbar = plt.colorbar()
-    plt.clim(-4.4,-2.5)
-    plt.scatter(np.log10(density[SFR_plot<=0]),np.log10(temperature[SFR_plot<=0]),s=1,c='gray')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.ylabel('$\log$ temperature [$\\rm K$]')
-    plt.ylim(3.5,5.5)
-    plt.xlim(-3.5,5)
-    plt.legend()
-
-    plt.subplot(3,3,7)
-    birthdensity = getbirthdensity(filename)*404.759
-    plt.hist(np.log10(birthdensity),bins=50,density=True,range=(np.log10(critden),5))
-    plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
-    plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
-    plt.xlim(np.log10(critden),5)
-    plt.legend()
-
-    plt.subplot(3,3,8)
-    sfrrate, gasdensity = getdensityrate(filename)
-    gasdensity *= 1e10/(1e3)**3*40.4759/mu
-    plt.hist2d(np.log10(gasdensity),np.log10(sfrrate), bins=50,range=[[-1.5,5],[-.5,2.5]])
-    cbar = plt.colorbar()
-    rhos = 10**np.linspace(-1,4,100)
-    rhoshigh = 10**np.linspace(2,5,100)
-    plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)),'k--',label='sSFR low density EAGLE')
-    plt.plot(np.log10(rhoshigh),np.log10(specificSFRhighden(rhoshigh)),'k--',label='sSFR high density EAGLE')
-    plt.xlabel('density [\\rm cm^{-3}]')
-    plt.ylabel('sSFR ($\\rm Gyr^{-1}$)')
-    plt.legend()
-    
-    plt.subplot(3,3,9)
-    timearray, SFH = getSFH(filename)
-    plt.plot(timearray*1e3,SFH, label='spart birth_time')
-    plt.xlabel('Time (Myr)')
-    plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
-    plt.legend()
-
-    if binsize==200:
-        plt.savefig('./ksplot/'+filename+'_0.1kpc.png')
-    elif binsize==20:
-        plt.savefig('./ksplot_1kpc/'+filename+'_0.1kpc.png')
-    plt.close()
-
-
-def getmassandsfr(filename):
-    
-    # Read the data
-    with h5.File(filename, "r") as f:
-        box_size = f["/Header"].attrs["BoxSize"][0]
-        coordinates= f["/PartType0/Coordinates"][:,:]
-        SFR = f["/PartType0/SFR"][:]
-        coordinates_star = f["/PartType4/Coordinates"][:,:]
-        masses_star = f["/PartType4/Masses"][:]
-        
-
-    absmaxz = 2 #kpc
-    absmaxxy = 10 #kpc
-
-    part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-        absmaxz) & (SFR>0) 
-    mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) < 
-        absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) < 
-        absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) < 
-        absmaxz)  
-
-    SFR_snap = np.sum(SFR[part_mask])
-    mass_star = np.sum(masses_star[mask])
-    return SFR_snap, mass_star
-
-def getsfrsnapwide():
-    
-    time = np.arange(1,101,1)
-    SFR_sparticles = np.zeros(100)
-    SFR_gparticles = np.zeros(100)
-    new_sparticles = np.zeros(100)
-    previous_mass = 0
-    previous_numb = 0
-    for i in range(1,100):
-        print(i)
-        # Read the data
-        filename= 'output_%04d.hdf5'%i
-        with h5.File(filename, "r") as f:
-            box_size = f["/Header"].attrs["BoxSize"][0]
-            coordinates= f["/PartType0/Coordinates"][:,:]
-            SFR = f["/PartType0/SFR"][:]
-            coordinates_star = f["/PartType4/Coordinates"][:,:]
-            masses_star = f["/PartType4/Masses"][:]
-            newstar = f["/PartType4/NewStarFlag"][:]
-
-        absmaxz = 2 #kpc
-        absmaxxy = 10 #kpc
-
-        part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) < 
-            absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) < 
-            absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) < 
-            absmaxz) & (SFR>0) 
-        mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) < 
-            absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) < 
-            absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) < 
-            absmaxz) & (newstar==1) 
-
-        SFR = SFR[part_mask]
-        masses_star = masses_star[mask]*1e4
-        newstar = newstar[mask]
-
-        total_SFR = np.sum(SFR)
-        total_mass = np.sum(masses_star)
-        new_spart = np.sum(newstar)
-        new_sparticles[i] = new_spart - previous_numb
-        SFR_sparticles[i] = total_mass - previous_mass 
-        SFR_gparticles[i] = total_SFR*10
-        previous_numb = new_spart
-        previous_mass = total_mass
-
-    factor = masses_star[0]
-    print(new_sparticles)
-    print(len(time),len(SFR_sparticles),len(SFR_gparticles))
-    return time[:-1], SFR_sparticles[1:], SFR_gparticles[1:], np.sqrt(new_sparticles[1:])*factor
-
-
-
-if __name__ == "__main__":
-
-    # Some global variables to define the whole run
-    input_file_name = sys.argv[1]
-    pixels = 1000 
-    num_procs = 4 # No. of processors to use
-    colour_map = "inferno"
-    
-    # Find the list of files and sort them in numerical order
-    file_list = sorted(glob.glob(input_file_name + "*.hdf5"), key=lambda x:
-            int(x[len(input_file_name):len(input_file_name)+4]))
-    
-    # get the image data 
-    #for i in [1,5,10,33,66,90]:
-    #    print('output_%04d.hdf5'%i)
-    #    makefullplot('output_%04d.hdf5'%i)
-    #    makefullplot('output_%04d.hdf5'%i,binsize=20)
-
-    '''
-    snapshots = 52 
-    mass = np.zeros(snapshots)
-    sfr = np.zeros(snapshots)
-    for i in range(1,snapshots):
-        if i!=6:
-            print(i)
-            sfr[i], mass[i] = getmassandsfr('output_%04d.hdf5'%i)
-        
-    plt.plot(np.log10(mass*1e10))
-    plt.xlabel('Snapshot number (Myr)')
-    plt.ylabel('Mass')
-    plt.savefig('Mass increase')
-    plt.close()
-
-    plt.plot(np.log10(sfr*10))
-    plt.xlabel('Snapshot number (Myr)')
-    plt.ylabel('SFR ($M_\odot / \\rm yr$)')
-    plt.savefig('SFR')
-    plt.close()
-    '''
-
-    for i in [1,5,10,33,66,90]:
-        if i!=6:
-            print('output_%04d.hdf5'%i)
-            makefullplot('output_%04d.hdf5'%i)
-            makefullplot('output_%04d.hdf5'%i,binsize=20)
-
-
-    time, SFR1, SFR2, SFR_error = getsfrsnapwide()
-    time2, SFR3 = getSFH('output_%04d.hdf5'%99)
-    plt.plot(time,SFR1,label='Using created star particles')
-    print(SFR_error)
-    plt.fill_between(time, SFR1+SFR_error,SFR1-SFR_error,label='Uncertainty',alpha=.3)
-    plt.plot(time,SFR2,label='Using SFR of gas particles')
-    plt.plot(time2*1e3,SFR3,label='Using birth_time of star particles')
-    plt.xlabel('Time (Myr)')
-    plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
-    plt.legend()
-    plt.show()
diff --git a/examples/IsolatedGalaxy_starformation_25p_intermediate/getIC.sh b/examples/IsolatedGalaxy_starformation_25p_intermediate/getIC.sh
deleted file mode 100755
index 1d5bcb30dc..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p_intermediate/getIC.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-wget https://home.strw.leidenuniv.nl/~nobels/data/starforming_intermediate.hdf5
diff --git a/examples/IsolatedGalaxy_starformation_25p_intermediate/isolated_galaxy.yml b/examples/IsolatedGalaxy_starformation_25p_intermediate/isolated_galaxy.yml
deleted file mode 100644
index 44ce055182..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p_intermediate/isolated_galaxy.yml
+++ /dev/null
@@ -1,114 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
-  UnitMass_in_cgs:     1.9891E43   # 10^10 solar masses 
-  UnitLength_in_cgs:   3.08567758E21   # 1 kpc 
-  UnitVelocity_in_cgs: 1E5   # km/s
-  UnitCurrent_in_cgs:  1   # Amperes
-  UnitTemp_in_cgs:     1   # Kelvin
-
-# Parameters for the self-gravity scheme
-Gravity:
-  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
-  eta:          0.025               # Constant dimensionless multiplier for time integration.
-  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.01 # Comoving softening length (in internal units).
-  max_physical_softening: 0.01    # Physical softening length (in internal units).
-
-# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
-TimeIntegration:
-  time_begin:        0.    # The starting time of the simulation (in internal units).
-  time_end:          0.1    # The end time of the simulation (in internal units).
-  dt_min:            1e-9  # The minimal time-step size of the simulation (in internal units).
-  dt_max:            1e-2  # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
-  basename:   output      # Common part of the name of output files
-  time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
-  delta_time: 0.001        # Time difference between consecutive outputs (in internal units)
-
-  
-# Parameters governing the conserved quantities statistics
-Statistics:
-  delta_time:           1e-2     # Time between statistics output
-  time_first:             0.     # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
-
-# Parameters related to the initial conditions
-InitialConditions:
-  file_name:  starforming_intermediate.hdf5 # The file to read
-  periodic:                    0    # Are we running with periodic ICs?
-
-# Parameters for the hydrodynamics scheme
-SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  h_tolerance:           1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
-  h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
-  max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
-  max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
-  minimal_temperature:   100        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
-  H_ionization_temperature: 1e4   # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
-
-EAGLECooling:
-  dir_name:                ../coolingtables/  # Location of the Wiersma+08 cooling tables
-  H_reion_z:               11.5              # Redshift of Hydrogen re-ionization
-  He_reion_z_centre:       3.5               # Redshift of the centre of the Helium re-ionization Gaussian
-  He_reion_z_sigma:        0.5               # Spread in redshift of the  Helium re-ionization Gaussian
-  He_reion_eV_p_H:         2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
-
-# Primordial abundances
-EAGLEChemistry:
-  init_abundance_metal:     0.0129          # Inital fraction of particle mass in *all* metals
-  init_abundance_Hydrogen:  0.7065       # Inital fraction of particle mass in Hydrogen
-  init_abundance_Helium:    0.2806        # Inital fraction of particle mass in Helium
-  init_abundance_Carbon:    0.00207        # Inital fraction of particle mass in Carbon
-  init_abundance_Nitrogen:  0.000836        # Inital fraction of particle mass in Nitrogen
-  init_abundance_Oxygen:    0.00549        # Inital fraction of particle mass in Oxygen
-  init_abundance_Neon:      0.00141        # Inital fraction of particle mass in Neon
-  init_abundance_Magnesium: 0.000591        # Inital fraction of particle mass in Magnesium
-  init_abundance_Silicon:   0.000683        # Inital fraction of particle mass in Silicon
-  init_abundance_Iron:      0.0011        # Inital fraction of particle mass in Iron
-
-# Hernquist potential parameters
-HernquistPotential:
-  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
-  position:        [0.,0.,0.]    # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
-  idealizeddisk:   1        # Run with an idealized galaxy disk
-  M200:            30.0   # M200 of the galaxy disk
-  h:               0.704    # reduced Hubble constant (value does not specify the used units!)
-  concentration:   7.1      # concentration of the Halo
-  diskfraction:              0.040   # Disk mass fraction
-  bulgefraction:              0.010   # Bulge mass fraction
-  timestep_mult:   0.01     # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
-  epsilon:         0.01      # Softening size (internal units)
- 
-# Schaye and Dalla Vecchia 2008 star formation
-SchayeSF:
-  thresh_MinOverDens:               57.7      # The critical density contrast to form stars
-  thresh_temp:                      1e5     # The critical temperature to form stars
-  fg:                               0.25    # The mass fraction in gas
-  SchmidtLawCoeff_MSUNpYRpKPC2:     1.515e-4  # The normalization of the Kennicutt-Schmidt law
-  SchmidtLawExponent:               1.4     # The power law of the Kennicutt-Schmidt law
-  SchmidtLawHighDensExponent:       2.0     # The high density exponent for the Kennicutt-Schmidt law
-  SchmidtLawHighDens_thresh_HpCM3:  1e3
-  Schaye2004:                       1       # whether to use the metallicity dependent critical star formation of Schaye (2004) (1) or not (0).
-  thresh_norm_HpCM3:                0.1      # Critical sf normalization to use (is not a normalization when Schaye2004=0, than it is the value.
-  thresh_max_norm_HpCM3:            10.0    # Maximum norm of the critical sf density
-  MetDep_Z0:                        0.002   # Scale metallicity to use for the equation (not used for Schaye2004=0)
-  MetDep_SFthresh_Slope:            -0.64   # Scaling of the critical density with the metallicity (not used for Schaye2004=0)
-  thresh_MaxPhysDensOn:             0       # Default is 0.
-  thresh_MaxOverDens_HpCM3:         1e5     # Density at which the SF law changes
-  EOS_Jeans_GammaEffective:         1.3333333 # The polytropic index
-  EOS_Jeans_TemperatureNorm_K:      1e5     # No idea how this works
-  EOS_JEANS_DensityNorm_HpCM3:      0.1      # No idea what the value is.
-
-# Parameters for the EAGLE "equation of state"
-EAGLEEntropyFloor:
-  Jeans_density_threshold_H_p_cm3: 0.1       # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
-  Jeans_over_density_threshold:    10.       # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
-  Jeans_temperature_norm_K:        8000      # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
-  Jeans_gamma_effective:           1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
-  Cool_density_threshold_H_p_cm3: 1e-5       # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
-  Cool_over_density_threshold:    10.        # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
-  Cool_temperature_norm_K:        8000       # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
-  Cool_gamma_effective:           1.         # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/IsolatedGalaxy_starformation_25p_intermediate/run.sh b/examples/IsolatedGalaxy_starformation_25p_intermediate/run.sh
deleted file mode 100755
index f16f8c170c..0000000000
--- a/examples/IsolatedGalaxy_starformation_25p_intermediate/run.sh
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/bin/bash
-
-if [ ! -e starforming_intermediate.hdf5 ]
-then 
-  ./getIC.sh
-fi 
-
-../swift --threads=32 --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
-- 
GitLab