Add a cleaned up example before removing the other two examples

examples/IsolatedGalaxy_starformation/README

+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.

examples/IsolatedGalaxy_starformation/getIC.sh

+#!/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

examples/IsolatedGalaxy_starformation/isolated_galaxy.yml

+# 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

examples/IsolatedGalaxy_starformation/ks_plotter.py

+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)

examples/IsolatedGalaxy_starformation/run.sh

+#!/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