diff --git a/examples/IsolatedGalaxy_starformation/README b/examples/IsolatedGalaxy_starformation/README
new file mode 100644
index 0000000000000000000000000000000000000000..16d2abde23b742e421bbd8c04a632edb0e35108a
--- /dev/null
+++ b/examples/IsolatedGalaxy_starformation/README
@@ -0,0 +1,40 @@
+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 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_starformation/getIC.sh b/examples/IsolatedGalaxy_starformation/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..197bc332fc1d2d140d572ecb4bfd8ca9b30e50b9
--- /dev/null
+++ b/examples/IsolatedGalaxy_starformation/getIC.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/fid.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/f10.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/f90.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres8.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres64.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres512.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/highres6.hdf5
diff --git a/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml
new file mode 100644
index 0000000000000000000000000000000000000000..fb2e582b756b7e74dac513957e5f2dd4f5ac4130
--- /dev/null
+++ b/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml
@@ -0,0 +1,111 @@
+# 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: fid.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: 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)
+
+# 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.3 # 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
+
+# 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/ks_plotter.py b/examples/IsolatedGalaxy_starformation/ks_plotter.py
new file mode 100644
index 0000000000000000000000000000000000000000..c046629c7069761f821b9b6a6ae87bfbb8b987bd
--- /dev/null
+++ b/examples/IsolatedGalaxy_starformation/ks_plotter.py
@@ -0,0 +1,316 @@
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py as h5
+from sphviewer.tools import QuickView
+
+# 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"]})
+
+
+filename = "output_0033.hdf5"
+
+f = h5.File(filename, "r")
+
+# Physical constants
+k_in_cgs = 1.38064852e-16
+mH_in_cgs = 1.6737236e-24
+year_in_cgs = 3600. * 24 * 365.
+Msun_in_cgs = 1.98848e33
+
+# Gemoetry info
+boxsize = f["/Header"].attrs["BoxSize"]
+centre = boxsize / 2.
+
+# 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)"]
+
+# Read parameters of the model
+KS_law_slope = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawExponent"])
+KS_law_norm = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawCoeff_MSUNpYRpKPC2"])
+KS_thresh_Z0 = float(f["/Parameters"].attrs["SchayeSF:MetDep_Z0"])
+KS_thresh_slope = float(f["/Parameters"].attrs["SchayeSF:MetDep_SFthresh_Slope"])
+KS_thresh_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_norm_HpCM3"])
+KS_gas_fraction = float(f["/Parameters"].attrs["SchayeSF:fg"])
+KS_thresh_max_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_max_norm_HpCM3"])
+KS_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
+EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
+
+# 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"][:]
+
+# Centre the box
+gas_pos[:,0] -= centre[0]
+gas_pos[:,1] -= centre[1]
+gas_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
+
+# Equations of state
+eos_cool_rho = np.logspace(-5, 5, 1000)
+eos_cool_T = eos_cool_rho**0. * 8000.
+eos_Jeans_rho = np.logspace(-1, 5, 1000)
+eos_Jeans_T = (eos_Jeans_rho/ 10**(-1))**(1./3.) * 8000.
+
+# 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(1e-4, 3e3)
+ylim(500., 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=1)
+plot(eos_Jeans_rho, eos_Jeans_T, 'k--', lw=1)
+scatter(gas_nH[gas_SFR > 0.], gas_T[gas_SFR > 0.], s=0.2)
+xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
+xlim(1e-4, 3e3)
+ylim(500., 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
+
+# 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.) & (gas_pos[:, 2] > -1.)
+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))
+
+qv = QuickView(
+ gas_pos,
+ mass=gas_mass,
+ r="infinity",
+ xsize=len(x_edges)-1,
+ ysize=len(y_edges)-1,
+ p=0, # Viewing angle theta
+ roll=0, # Viewing angle phi
+ plot=False,
+ logscale=False,
+ hsml=gas_hsml
+)
+
+map_mass2 = qv.get_image()
+extent_mass = qv.get_extent()
+
+gas_SFR[gas_SFR<=0] = 1e-10
+
+qv = QuickView(
+ gas_pos,
+ mass=gas_SFR,
+ r="infinity",
+ xsize=len(x_edges)-1,
+ ysize=len(y_edges)-1,
+ p=0, # Viewing angle theta
+ roll=0, # Viewing angle phi
+ plot=False,
+ logscale=False,
+ hsml=gas_hsml
+)
+
+map_SFR2 = qv.get_image()
+extent_SFR = qv.get_extent()
+
+
+# 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)
+
+# Mass map 2
+figure()
+subplot(111)
+imshow(np.log10(map_mass2),extent=extent_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_SPHVIEWER.png", dpi=200)
+
+# SF map
+figure()
+subplot(111)
+pcolormesh(x_edges, y_edges, np.log10(map_SFR), vmax = -.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)
+
+plot_map_SFR2 = np.zeros(np.shape(map_SFR2))
+plot_map_SFR2[map_SFR2>1e-6] = map_SFR2[map_SFR2>1e-6]
+# SF map 2
+figure()
+subplot(111)
+imshow(np.log10(plot_map_SFR2),extent=extent_SFR, vmax = -.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_SPHVIEWER.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. * 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, KS_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()
+
+plot_map_SFR2[plot_map_SFR2<=0] = 1e-6
+
+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, KS_gamma_effective, EAGLE_Z), fontsize=7)
+scatter(map_mass2.flatten() / 1e6, plot_map_SFR2.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_SPHVIEWER.png", dpi=200)
diff --git a/examples/IsolatedGalaxy_starformation/run.sh b/examples/IsolatedGalaxy_starformation/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3f9c3c6331034eb630503973106f47b3caa93d46
--- /dev/null
+++ b/examples/IsolatedGalaxy_starformation/run.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+if [ ! -e fid.hdf5 ]
+then
+ echo "Fetching initial conditions for the isolated galaxy example..."
+ ./getIC.sh
+fi
+
+../../swift --threads=32 --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
+
+python3 ks_plotter.py