diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/README b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/README
new file mode 100644
index 0000000000000000000000000000000000000000..719ed3356701983e60e0f032791e4bb9a0524978
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/README
@@ -0,0 +1,43 @@
+Isolated Galaxy generated by the MakeNewDisk code from Springel, Di Matteo &
+Hernquist (2005). The done analysis in this example is similar to the work done
+by Schaye and Dalla Vecchia (2008) (After this SD08). The default example runs
+the simulation for a galaxy with similar mass of their fiducial model and should
+produce plots similar to their middle pannel Figure 4. The code needs to be
+configured to run with the Hernquist external potential as well as the cooling &
+star-formation model of interest. Using the EAGLE model allows to reproduce the
+results of SD08.
+
+The code can also be run for other situations to check to verify the law using
+different parameters, changes that were done in SD08 are given by:
+ - gas fraction of 10% instead of 30%, change the IC to f10.hdf5, see getIC.sh,
+ should reproduce something similar to Figure 4 left hand pannel. Requires
+ change of fg=.1
+ - gas fraction of 90% instead of 30%, change the IC to f90.hdf5, see getIC.sh,
+ should reproduce something similar to Figure 4 right hand pannel. Requires
+ change of fg=.9
+ - Changing the effective equation of state to adiabatic, Jeans_gamma_effective
+ = 1.666667. Should result in something similar to Figure 5 left hand pannel
+ of SD08.
+ - Changing the effective equation of state to isothermal, Jeans_gamma_effective
+ = 1.0000. Should result in something similar to Figure 5 middle hand pannel
+ of SD08.
+ - Changing the slope of the Kennicutt-Schmidt law to 1.7, SchmidtLawExponent =
+ 1.7, this should result in a plot similar to Figure 6 of SD08.
+ - Increasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 1.0,
+ should reproduce plot similar to Figure 7.
+ - Decreasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 0.01,
+ should reproduce plot similar to Figure 7.
+ - Running with a lower resultion of a factor 8, change the IC to lowres8.hdf5,
+ see getIC.sh.
+ - Running with a lower resultion of a factor 64, change the IC to lowres64.hdf5,
+ see getIC.sh.
+ - Running with a lower resultion of a factor 512, change the IC to lowres512.hdf5,
+ see getIC.sh.
+
+Other options to verify the correctness of the code is by chaning the following
+parameters:
+ - Changing the normalization to A/2 or 2A.
+ - Running the code with zero metallicity.
+ - Running the code with a factor 6 higher resolution idealized disks, use
+ highres6.hdf5, see getIC.sh.
+ - Running with different SPH schemes like Anarchy-PU.
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/SFH.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/SFH.py
new file mode 100644
index 0000000000000000000000000000000000000000..fa9d9258530396fb7f95237a45af5db9c0da4603
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/SFH.py
@@ -0,0 +1,100 @@
+"""
+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 tqdm import tqdm
+
+
+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.0) > -absmaxxy)
+ & ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
+ & ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
+ & ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
+ & ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
+ & ((coordinates[:, 2] - box_size / 2.0) < 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, 0.1])
+ values = histogram[0]
+ xvalues = (histogram[1][:-1] + histogram[1][1:]) / 2.0
+ return xvalues, values
+
+
+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 tqdm(range(1, 100)):
+ # 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"][:]
+
+ absmaxz = 2 # kpc
+ absmaxxy = 10 # kpc
+
+ part_mask = (
+ ((coordinates[:, 0] - box_size / 2.0) > -absmaxxy)
+ & ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
+ & ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
+ & ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
+ & ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
+ & ((coordinates[:, 2] - box_size / 2.0) < absmaxz)
+ & (SFR > 0)
+ )
+
+ SFR = SFR[part_mask]
+
+ total_SFR = np.sum(SFR)
+ SFR_gparticles[i] = total_SFR * 10
+
+ return time[:-1], SFR_gparticles[1:]
+
+
+if __name__ == "__main__":
+
+ time, SFR1 = getsfrsnapwide() # , SFR2, SFR_error = getsfrsnapwide()
+ time2, SFR3 = getSFH("output_%04d.hdf5" % 100)
+ plt.plot(time2[1:] * 1e3, SFR3[1:], label="Using birth_time of star particles")
+ plt.plot(time, SFR1, label="Using SFR of gas particles", color="g")
+ plt.xlabel("Time (Myr)")
+ plt.ylabel("SFH ($\\rm M_\odot \\rm yr^{-1}$)")
+ plt.ylim(0, 20)
+ plt.legend()
+ plt.savefig("SFH.png")
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getEagleCoolingTable.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getEagleCoolingTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5cfd93ef0f4603e40b7675f3f2c254b2250f699f
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getEagleCoolingTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/EAGLE/coolingtables.tar.gz
+tar -xf coolingtables.tar.gz
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cc868b42a4755397308d9bb7a1bba30d9a15e169
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/fid.hdf5
+# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/f10.hdf5
+# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/f90.hdf5
+ wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres8.hdf5
+# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres64.hdf5
+#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres512.hdf5
+# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/highres6.hdf5
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2bdf1b9cab8b362426631f9e7aaba93bf6a6fb02
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
@@ -0,0 +1,117 @@
+# 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: lowres8.hdf5 # The file to read
+ periodic: 0 # Are we running with periodic ICs?
+ stars_smoothing_length: 0.5
+
+# 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_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 10.
+
+# Standard EAGLE cooling options
+EAGLECooling:
+ dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
+ H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
+ H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
+ 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: 137.0 # M200 of the galaxy disk
+ h: 0.704 # reduced Hubble constant (value does not specify the used units!)
+ concentration: 9.0 # concentration of the Halo
+ diskfraction: 0.040 # Disk mass fraction
+ bulgefraction: 0.014 # 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)
+
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# 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
+
+EagleStellarEvolution:
+ filename: /cosma5/data/Eagle/BG_Tables/YieldTables/
+ imf_model: Chabrier
+
+EAGLEFeedback:
+ lifetime_flag: 2
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..89a87923148cb6872ab17b6d7229aef597ef3287
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
@@ -0,0 +1,368 @@
+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py as h5
+
+# Plot parameters
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 9,
+ "legend.fontsize": 9,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (3.15, 3.15),
+ "figure.subplot.left": 0.15,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.13,
+ "figure.subplot.top": 0.99,
+ "figure.subplot.wspace": 0.15,
+ "figure.subplot.hspace": 0.12,
+ "lines.markersize": 6,
+ "lines.linewidth": 2.0,
+ "text.latex.unicode": True,
+}
+rcParams.update(params)
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
+
+snap = int(sys.argv[1])
+filename = "output_%.4d.hdf5"%snap
+
+f = h5.File(filename, "r")
+
+# Physical constants
+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"]
+centre = boxsize / 2.0
+
+# Read units
+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 SF model
+KS_law_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_exponent"])
+KS_law_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_normalisation"])
+KS_thresh_Z0 = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_Z0"])
+KS_thresh_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_slope"])
+KS_thresh_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"])
+KS_gas_fraction = float(f["/Parameters"].attrs["EAGLEStarFormation:gas_fraction"])
+KS_thresh_max_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"])
+KS_high_den_thresh = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"])
+KS_law_slope_high_den = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"])
+EOS_gamma_effective = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"])
+EOS_density_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"])
+EOS_temp_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"])
+
+# Read reference metallicity
+EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
+
+# Read parameters of the entropy floor
+EAGLEfloor_Jeans_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
+EAGLEfloor_Jeans_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
+EAGLEfloor_Jeans_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
+EAGLEfloor_cool_rho_norm = 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"])
+
+# Properties of the KS law
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
+gamma = 5./3.
+EOS_press_norm = k_in_cgs * EOS_temp_norm * EOS_density_norm
+
+# Star formation threshold
+SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0)**(KS_thresh_slope)
+
+# Read gas properties
+gas_pos = f["/PartType0/Coordinates"][:, :]
+gas_mass = f["/PartType0/Masses"][:]
+gas_rho = f["/PartType0/Density"][:]
+gas_T = f["/PartType0/Temperature"][:]
+gas_SFR = f["/PartType0/SFR"][:]
+gas_XH = f["/PartType0/ElementAbundance"][:, 0]
+gas_Z = f["/PartType0/Metallicity"][:]
+gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_sSFR = gas_SFR / gas_mass
+
+# Read the Star properties
+stars_pos = f["/PartType4/Coordinates"][:, :]
+stars_BirthDensity = f["/PartType4/BirthDensity"][:]
+stars_BirthTime = f["/PartType4/BirthTime"][:]
+stars_XH = f["/PartType4/ElementAbundance"][:,0]
+
+# Centre the box
+gas_pos[:, 0] -= centre[0]
+gas_pos[:, 1] -= centre[1]
+gas_pos[:, 2] -= centre[2]
+
+stars_pos[:,0] -= centre[0]
+stars_pos[:,1] -= centre[1]
+stars_pos[:,2] -= centre[2]
+
+# Turn the mass into better units
+gas_mass *= unit_mass_in_cgs / Msun_in_cgs
+
+# Turn the SFR into better units
+gas_SFR = np.maximum(gas_SFR, np.zeros(np.size(gas_SFR)))
+gas_SFR /= unit_time_in_cgs / year_in_cgs
+gas_SFR *= unit_mass_in_cgs / Msun_in_cgs
+
+# Make it a Hydrogen number density
+gas_nH = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
+gas_nH /= mH_in_cgs
+gas_nH *= gas_XH
+
+stars_BirthDensity *= unit_mass_in_cgs / unit_length_in_cgs ** 3
+stars_BirthDensity /= mH_in_cgs
+stars_BirthDensity *= stars_XH
+
+# Equations of state
+eos_cool_rho = np.logspace(-5, 5, 1000)
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (eos_cool_rho / EAGLEfloor_cool_rho_norm) ** ( EAGLEfloor_cool_gamma_effective - 1.0 )
+eos_Jeans_rho = np.logspace(-1, 5, 1000)
+eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm) ** (EAGLEfloor_Jeans_gamma_effective - 1.0 )
+
+########################################################################3
+
+# Plot the phase space diagram
+figure()
+subplot(111, xscale="log", yscale="log")
+plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
+plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
+scatter(gas_nH, gas_T, s=0.2)
+xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
+xlim(3e-6, 3e3)
+ylim(500.0, 2e5)
+savefig("rhoT.png", dpi=200)
+
+# Plot the phase space diagram for SF gas
+figure()
+subplot(111, xscale="log", yscale="log")
+plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
+plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
+plot([SF_thresh, SF_thresh], [1, 1e10], "k:", lw=0.6)
+text(SF_thresh*0.9, 2e4, "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$"%SF_thresh, fontsize=8, rotation=90, ha="right", va="bottom")
+scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
+xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
+xlim(3e-6, 3e3)
+ylim(500.0, 2e5)
+savefig("rhoT_SF.png", dpi=200)
+
+########################################################################3
+
+# 3D Density vs SFR
+figure()
+subplot(111, xscale="log", yscale="log")
+scatter(gas_nH, gas_SFR, s=0.2)
+plot([1, 100], 2e-5 * np.array([1, 100]) ** 0.266667, "k--", lw=1)
+xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm SFR}~[{\\rm M_\\odot~\\cdot~yr^{-1}}]$", labelpad=-7)
+xlim(1e-4, 3e3)
+ylim(8e-6, 2.5e-4)
+savefig("rho_SFR.png", dpi=200)
+
+########################################################################3
+
+star_mask = (
+ (stars_pos[:, 0] > -15)
+ & (stars_pos[:, 0] < 15)
+ & (stars_pos[:, 1] > -15)
+ & (stars_pos[:, 1] < 15)
+ & (stars_pos[:, 2] < 1.0)
+ & (stars_pos[:, 2] > -1.0)
+)
+
+stars_BirthDensity = stars_BirthDensity[star_mask]
+#stars_BirthFlag = stars_BirthFlag[star_mask]
+stars_BirthTime = stars_BirthTime[star_mask]
+
+# Histogram of the birth density
+figure()
+subplot(111, xscale="linear", yscale="linear")
+hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5])
+xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm 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 = EOS_press_norm * ( rhos / EOS_density_norm)**(EOS_gamma_effective)
+P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm)**(EOS_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)**EOS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+
+# density - sSFR plane
+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=2)
+ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=0)
+xticks([-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"])
+yticks([0, 1, 2], ["$10^0$", "$10^1$", "$10^2$"])
+xlim(-1.4, 4.9)
+ylim(-0.5, 2.2)
+savefig("density-sSFR.png", dpi=200)
+
+########################################################################3
+
+# Select gas in a pillow box around the galaxy
+mask = (
+ (gas_pos[:, 0] > -15)
+ & (gas_pos[:, 0] < 15)
+ & (gas_pos[:, 1] > -15)
+ & (gas_pos[:, 1] < 15)
+ & (gas_pos[:, 2] < 1.0)
+ & (gas_pos[:, 2] > -1.0)
+)
+gas_pos = gas_pos[mask, :]
+gas_SFR = gas_SFR[mask]
+gas_nH = gas_nH[mask]
+gas_rho = gas_rho[mask]
+gas_T = gas_T[mask]
+gas_mass = gas_mass[mask]
+gas_Z = gas_Z[mask]
+gas_hsml = gas_hsml[mask]
+
+
+# Make a crude map of the gas
+figure()
+subplot(111)
+scatter(gas_pos[:, 0], gas_pos[:, 1], s=0.1)
+xlabel("${\\rm Pos}~x~[{\\rm kpc}]$", labelpad=0)
+ylabel("${\\rm Pos}~y~[{\\rm kpc}]$", labelpad=-3)
+xlim(-12, 12)
+ylim(-12, 12)
+savefig("face_on.png", dpi=200)
+
+figure()
+subplot(111)
+scatter(gas_pos[:, 0], gas_pos[:, 2], s=0.1)
+xlabel("${\\rm Pos}~x~[{\\rm kpc}]$", labelpad=0)
+ylabel("${\\rm Pos}~z~[{\\rm kpc}]$", labelpad=-3)
+xlim(-12, 12)
+ylim(-12, 12)
+savefig("edge_on.png", dpi=200)
+
+# Now a SF map
+rcParams.update({"figure.figsize": (4.15, 3.15)})
+figure()
+subplot(111)
+scatter(gas_pos[:, 0], gas_pos[:, 1], s=0.1, c=gas_SFR)
+xlabel("${\\rm Pos}~x~[{\\rm kpc}]$", labelpad=0)
+ylabel("${\\rm Pos}~y~[{\\rm kpc}]$", labelpad=-3)
+colorbar()
+xlim(-12, 12)
+ylim(-12, 12)
+savefig("SF_face_on.png", dpi=200)
+
+
+########################################################################3
+
+# Bin the data in kpc-size patches
+
+x_edges = np.linspace(-15, 15, 31)
+y_edges = np.linspace(-15, 15, 31)
+
+map_mass, _, _, _ = stats.binned_statistic_2d(
+ gas_pos[:, 0], gas_pos[:, 1], gas_mass, statistic="sum", bins=(x_edges, y_edges)
+)
+map_SFR, _, _, _ = stats.binned_statistic_2d(
+ gas_pos[:, 0], gas_pos[:, 1], gas_SFR, statistic="sum", bins=(x_edges, y_edges)
+)
+
+# Mass map
+figure()
+subplot(111)
+pcolormesh(x_edges, y_edges, np.log10(map_mass))
+colorbar()
+xlim(-12, 12)
+ylim(-12, 12)
+xlabel("${\\rm Pos}~x~[{\\rm kpc}]$", labelpad=0)
+ylabel("${\\rm Pos}~y~[{\\rm kpc}]$", labelpad=-3)
+savefig("Map_mass.png", dpi=200)
+
+# SF map
+figure()
+subplot(111)
+pcolormesh(x_edges, y_edges, np.log10(map_SFR), vmax=-0.5, vmin=-4.5)
+colorbar()
+xlim(-12, 12)
+ylim(-12, 12)
+xlabel("${\\rm Pos}~x~[{\\rm kpc}]$", labelpad=0)
+ylabel("${\\rm Pos}~y~[{\\rm kpc}]$", labelpad=-3)
+savefig("Map_SFR.png", dpi=200)
+
+#########################################################################
+
+# Give a minimum SF surface density for the plots
+map_SFR[map_SFR < 1e-6] = 1e-6
+
+# Theoretical threshold (assumes all gas has the same Z)
+KS_n_thresh = KS_thresh_norm * (gas_Z[0] / KS_thresh_Z0) ** KS_thresh_slope
+if np.isfinite(KS_n_thresh) == False:
+ KS_n_thresh = KS_thresh_max_norm
+KS_sigma_thresh = 29.0 * np.sqrt(KS_gas_fraction) * np.sqrt(KS_n_thresh)
+
+# Theoretical KS law
+KS_sigma_mass = np.logspace(-1, 3, 100)
+KS_sigma_SFR = KS_law_norm * KS_sigma_mass ** KS_law_slope
+
+# KS relation
+rcParams.update({"figure.figsize": (3.15, 3.15), "figure.subplot.left": 0.18})
+figure()
+subplot(111, xscale="log", yscale="log")
+plot(KS_sigma_mass, KS_sigma_SFR, "k--", lw=0.6)
+plot([KS_sigma_thresh, KS_sigma_thresh], [1e-8, 1e8], "k--", lw=0.6)
+text(
+ KS_sigma_thresh * 0.95,
+ 2.2,
+ "$\\Sigma_{\\rm c} = %.2f~{\\rm M_\\odot\\cdot pc^{-2}}$" % KS_sigma_thresh,
+ va="top",
+ ha="right",
+ rotation=90,
+ fontsize=7,
+)
+text(16, 10 ** (-3.5), "$n_{\\rm H,c} = %.3f~{\\rm cm^{-3}}$" % KS_n_thresh, fontsize=7)
+text(
+ 16,
+ 2e-6,
+ "${\\rm K\\textendash S~law}$:\n$\\Sigma_{\\rm SFR} = A \\times \\Sigma_g^n$\n$n=%.1f$\n$A=%.3f\\times10^{-4}~{\\rm M_\\odot / yr^{1} / kpc^{2}}$\n$f_{\\rm g} = %.1f$\n$\gamma_{\\rm eos} = %.3f$\n$Z=%1.4f$"
+ % (
+ KS_law_slope,
+ KS_law_norm * 10 ** 4,
+ KS_gas_fraction,
+ EOS_gamma_effective,
+ EAGLE_Z,
+ ),
+ fontsize=7,
+)
+scatter(map_mass.flatten() / 1e6, map_SFR.flatten(), s=0.4)
+xlim(0.3, 900)
+ylim(3e-7, 3)
+xlabel("$\\Sigma_g~[{\\rm M_\\odot\\cdot pc^{-2}}]$", labelpad=0)
+ylabel(
+ "$\\Sigma_{\\rm SFR}~[{\\rm M_\\odot \\cdot yr^{-1} \\cdot kpc^{-2}}]$", labelpad=0
+)
+savefig("KS_law.png", dpi=200)
+close()
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..67da3c390be1240323941b835e056dcd6e27feed
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
@@ -0,0 +1,155 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py
+import numpy as np
+import glob
+import os.path
+
+def find_indices(a,b):
+ result = np.zeros(len(b))
+ for i in range(len(b)):
+ result[i] = ((np.where(a == b[i]))[0])[0]
+
+ return result
+
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.1,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.1,
+'figure.subplot.top' : 0.95,
+'figure.subplot.wspace' : 0.2,
+'figure.subplot.hspace' : 0.2,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+# Number of snapshots and elements
+newest_snap_name = max(glob.glob('output_*.hdf5'), key=os.path.getctime)
+n_snapshots = int(newest_snap_name.replace('output_','').replace('.hdf5','')) + 1
+n_elements = 9
+
+# Read the simulation data
+sim = h5py.File("output_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+stellar_mass = sim["/PartType4/Masses"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+unit_temp_in_cgs = sim["/Units"].attrs["Unit temperature in cgs (U_T)"]
+unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
+unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
+unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
+unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
+Gyr_in_cgs = 3.155e16
+Msun_in_cgs = 1.989e33
+
+box_energy = zeros(n_snapshots)
+box_mass = zeros(n_snapshots)
+box_star_mass = zeros(n_snapshots)
+box_metal_mass = zeros(n_snapshots)
+element_mass = zeros((n_snapshots,n_elements))
+t = zeros(n_snapshots)
+
+# Read data from snapshots
+for i in range(n_snapshots):
+ print("reading snapshot "+str(i))
+ # Read the simulation data
+ sim = h5py.File("output_%04d.hdf5"%i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
+ #ids = sim["/PartType0/ParticleIDs"][:]
+
+ masses = sim["/PartType0/Masses"][:]
+ box_mass[i] = np.sum(masses)
+
+ star_masses = sim["/PartType4/Masses"][:]
+ box_star_mass[i] = np.sum(star_masses)
+
+ metallicities = sim["/PartType0/Metallicity"][:]
+ box_metal_mass[i] = np.sum(metallicities * masses)
+
+ internal_energies = sim["/PartType0/InternalEnergy"][:]
+ box_energy[i] = np.sum(masses * internal_energies)
+
+# Plot the interesting quantities
+figure()
+
+# Box mass --------------------------------
+subplot(221)
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
+xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
+ylabel("Change in total gas particle mass (Msun)", labelpad=2)
+ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+# Box metal mass --------------------------------
+subplot(222)
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_metal_mass[1:] - box_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
+ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
+ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+# Box energy --------------------------------
+subplot(223)
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_energy[1:] - box_energy[0])* unit_energy_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
+ylabel("Change in total energy of gas particles (erg)", labelpad=2)
+ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+# Box mass --------------------------------
+subplot(224)
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='gas')
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='r', marker = "*", ms=0.5, label='stars')
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1] + box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='g', marker = "*", ms=0.5, label='total')
+xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
+ylabel("Change in total gas particle mass (Msun)", labelpad=2)
+ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+legend()
+
+savefig("box_evolution.png", dpi=200)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..fd1bbf56ed99a0fb0a4bfcdd9a43e80cbe7ecee5
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+if [ ! -e lowres8.hdf5 ]
+then
+ echo "Fetching initial conditions for the isolated galaxy example..."
+ ./getIC.sh
+fi
+
+if [ ! -e cooling_tables/ ]
+then
+ echo "Fetching EAGLE cooling tables for the isolated galaxy example..."
+ ./getEagleCoolingTable.sh
+fi
+
+../../swift --threads=32 --feedback --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
+
+# Kennicutt-Schmidt law plot
+python3 plotSolution.py
+
+# Plot that the random star formation matches the expected SFH
+python3 SFH.py