diff --git a/Fig4_runaway_collapse_zone_examples.py b/Fig4_runaway_collapse_zone_examples.py
new file mode 100644
index 0000000000000000000000000000000000000000..aaf000f48a28a9dd42c97a11a04babe9781ebf61
--- /dev/null
+++ b/Fig4_runaway_collapse_zone_examples.py
@@ -0,0 +1,271 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 2 10:04:11 2022
+
+@author: sylviaploeckinger
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from unyt import msun, pc
+from unyt import proton_mass_cgs as mH
+import matplotlib.pylab as pl
+from astropy.cosmology import Planck18
+import matplotlib.patches as mpatches
+from matplotlib.collections import PatchCollection
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'large',
+ 'axes.labelsize': 'x-large',
+ 'axes.titlesize':'x-large',
+ 'figure.titlesize':'x-large',
+ 'xtick.labelsize':'large',
+ 'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+import sys
+sys.path.insert(0, "helpers/")
+from equations import Jeans_length_0
+from equations import Jeans_length_W
+from equations import density_critical_hmin
+from equations import Jeans_mass_0
+from equations import Jeans_mass_P
+from equations import Jeans_mass_W
+
+from simulation_plotparameter import color_zone_C
+
+eta_res = 1.2348
+XH = 0.74
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig4_runaway_collapse_zone_examples.png"
+
+############################################################################
+# Read data for the Springel & Hernquist 2003 effective equation of state
+# and convert the overdensity to a physical density for z=0
+############################################################################
+SH03_logrho, SH03_logT = np.loadtxt('data/Springel_Hernquist_temperature_data.dat', unpack = True, delimiter=',')
+log_avgdens_z0 = np.log10((Planck18.critical_density0 * Planck18.Ob0 * XH / mH).value)
+
+############################################################################
+# Set some resolution parameters for simulation projects from the literature
+############################################################################
+APOSTLEL1_eps = 134. # pc
+APOSTLEL1_mB = 1.e4 # Msun
+APOSTLEL1_hmin = 1.55 * 0.1 * APOSTLEL1_eps
+
+APOSTLEL2_eps = 307. # pc
+APOSTLEL2_mB = 1.2e5 # Msun
+APOSTLEL2_hmin = 1.55 * 0.1 * APOSTLEL2_eps
+
+APOSTLEL3_eps = 711. # pc
+APOSTLEL3_mB = 1.5e6 # Msun
+APOSTLEL3_hmin = 1.55 * 0.1 * APOSTLEL3_eps
+
+EAGLE_eps = 700. # pc
+EAGLE_mB = 1.81e6 # Msun
+EAGLE_hmin = 1.55 * 0.1 * EAGLE_eps
+
+EAGLE_highres_eps = 350.
+EAGLE_highres_mB = 2.26e5
+EAGLE_highres_hmin = 1.55 * 0.1 * EAGLE_highres_eps
+
+#####################################
+#Pillepich et al. (2019) Table 1
+#####################################
+TNG50_eps = 72.
+TNG50_mB = 8.5e4
+TNG50_hmin_avgcellsizeSF = 140. * 4. # force resolution needs ~ 4 resolution elements, similar to number of particles in kernel [~60**(1/3)]
+TNG50_hmin_mincellsize = 6.5 * 4.
+
+TNG100_eps = 190.
+TNG100_mB = 1.4e6
+TNG100_hmin_avgcellsizeSF = 355. * 4.
+
+TNG300_eps = 370.
+TNG300_mB = 5.9e7
+TNG300_hmin_avgcellsizeSF = 715. * 4.
+
+############################################################################
+# Set some general properties of the gas phases for reference
+############################################################################
+WNM_nHmin = 0.1
+WNM_nHmax = 0.5
+WNM_Tmin = 6000.
+WNM_Tmax = 8000.
+
+CNM_nHmin = 20.
+CNM_nHmax = 50.
+CNM_Tmin = 50.
+CNM_Tmax = 100.
+
+MC_nHmin = 1.e2
+MC_nHmax = 1.e6
+MC_Tmin = 15.
+MC_Tmax = 20.
+
+############################################################################
+# set plot limits and define 2D arrays for contour plots
+############################################################################
+
+color_zone = color_zone_C
+color_res = '#00CCCC'
+
+lognH_min = -6.
+lognH_max = 8.
+dlognH = 0.1
+
+logT_min = 1.
+logT_max = 8.
+dlogT = 0.1
+
+xticks = np.arange(-6., 10., 2)
+yticks = np.arange(1., 9., 1)
+
+lognH_arr = np.arange(lognH_min, lognH_max + dlognH, dlognH)
+logT_arr = np.arange(logT_min , logT_max + dlogT, dlogT)
+logT_arr_large = np.arange(logT_min , logT_max + 4. + dlogT, dlogT)
+
+lognH_2Darr = np.tile(lognH_arr, (len(logT_arr), 1))
+logT_2Darr = (np.tile(logT_arr, (len(lognH_arr), 1))).T
+
+nH_2Darr = np.power(10., lognH_2Darr)
+T_2Darr = np.power(10., logT_2Darr)
+
+eplot = [lognH_arr[0], lognH_arr[-1], logT_arr[0], logT_arr[-1]]
+
+MJmin = np.log10(1.e1)
+MJmax = np.log10(1.e10)
+
+manual_locations = [ (4.5,1.5), \
+ (3.0,1.5), \
+ (1.0,1.5), \
+ (-1.,1.5), \
+ (-2.0,1.5), \
+ (-4.5,1.5), \
+ (-2.0,3.), \
+ (-4.5,3.), \
+ (-2.0,4.5), \
+ (-4.5,4.) ]
+
+
+############################################################################
+# make subplot:
+# axloc ... matplotlib axes object
+# eps ... gravitational softening scale in pc
+# h_min ... minimum smoothing length in pc
+# mB ... gas particle mass resolution
+# title ... title for plot
+############################################################################
+def make_subplot(axloc, eps, h_min, mB, title):
+ axloc.set_title(title)
+ axloc.set_ylabel('log T [K]')
+ axloc.set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ axloc.yaxis.set_ticks(yticks)
+ axloc.xaxis.set_ticks(xticks)
+ axloc.set_ylim(logT_min, logT_max)
+ axloc.set_xlim(lognH_min, lognH_max)
+
+ logMJ_W_2D = np.log10(Jeans_mass_W(nH_2Darr, T_2Darr, 3.*eps))
+ CF_W = axloc.contour(lognH_2Darr, logT_2Darr, logMJ_W_2D, levels = np.arange(MJmin, MJmax + 1, 1), linestyles = 'solid', colors = 'black', alpha=0.3)
+ axloc.clabel(CF_W, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations[5:])
+
+ zone_of_avoidance = np.zeros_like(lognH_2Darr)
+ zone_of_avoidance[:,:] = 1.
+ zone_of_avoidance[(Jeans_length_W(nH_2Darr, T_2Darr, eps) / Jeans_length_0(nH_2Darr, T_2Darr) > 1.)] = 0.
+
+ # gravity resolved
+ #axloc.contourf(lognH_2Darr, logT_2Darr, Jeans_length_0(nH_2Darr, T_2Darr) / (1.5 * eps), levels = np.asarray([1.01, 1.e10]), colors = color_res, alpha = 0.3)
+ #axloc.contour(lognH_2Darr, logT_2Darr, Jeans_length_0(nH_2Darr, T_2Darr) / (1.5 * eps), levels = np.asarray([1.01]), colors = color_res)
+
+ if 'EAGLE' in title or 'APOSTLE' in title:
+ axloc.plot(lognH_arr[lognH_arr>=-1.], 1./3. * (lognH_arr[lognH_arr>=-1.] +1.) + np.log10(8000.), color = 'grey', lw = 5, alpha = 1.0, label = 'T$_{\mathrm{eff}}$')
+ elif 'TNG' in title:
+ axloc.plot(SH03_logrho + log_avgdens_z0, SH03_logT, color = 'grey', lw = 5, alpha = 1.0, label = 'SH2003')
+
+ zone_of_avoidance = np.zeros_like(lognH_2Darr)
+ zone_of_avoidance[:,:] = 1.
+ zone_of_avoidance[(Jeans_length_W(nH_2Darr, T_2Darr, eps) / h_min < 1.) & (lognH_2Darr > np.log10(density_critical_hmin(mB, h_min)))] = 0.
+ CS = axloc.contour(lognH_2Darr, logT_2Darr, zone_of_avoidance, levels = np.asarray([0.99]), colors = color_zone, linestyles = 'dashed')
+ if len(CS.allsegs[0]) > 0:
+ axloc.text(0.99, 0.15, 'runaway\ncollapse zone', va = 'bottom', ha = 'right', transform=axloc.transAxes)
+ CS = axloc.contourf(lognH_2Darr, logT_2Darr, zone_of_avoidance, levels = np.asarray([0., 0.99]), colors = color_zone, alpha = 0.8)
+
+ if 'EAGLE' in title or 'APOSTLE' in title or 'TNG' in title:
+ axloc.legend(loc = 'upper left', framealpha = 1.0)
+
+
+ # add the patches for the gas phases
+ WNM = mpatches.FancyBboxPatch((np.log10(WNM_nHmin),np.log10(WNM_Tmin)), \
+ (np.log10(WNM_nHmax) - np.log10(WNM_nHmin)), \
+ (np.log10(WNM_Tmax) - np.log10(WNM_Tmin)), \
+ boxstyle=mpatches.BoxStyle("Round", pad=0.2), facecolor = 'black', edgecolor = 'none', alpha = 0.6)
+ CNM = mpatches.FancyBboxPatch((np.log10(CNM_nHmin),np.log10(CNM_Tmin)), \
+ (np.log10(CNM_nHmax) - np.log10(CNM_nHmin)), \
+ (np.log10(CNM_Tmax) - np.log10(CNM_Tmin)), \
+ boxstyle=mpatches.BoxStyle("Round", pad=0.2), facecolor = 'black', edgecolor = 'none', alpha = 0.6)
+ MCs = mpatches.FancyBboxPatch((np.log10(MC_nHmin),np.log10(MC_Tmin)), \
+ (np.log10(MC_nHmax) - np.log10(MC_nHmin)), \
+ (np.log10(MC_Tmax) - np.log10(MC_Tmin)), \
+ boxstyle=mpatches.BoxStyle("Round", pad=0.2), facecolor = 'black', edgecolor = 'none', alpha = 0.6)
+ axloc.add_patch(WNM)
+ axloc.add_patch(CNM)
+ axloc.add_patch(MCs)
+ #axloc.annotate('WNM', (np.log10(WNM_nHmax), np.log10(WNM_Tmax)), xytext=(np.log10(WNM_nHmax)+1., np.log10(WNM_Tmax)-0.5),\
+ # arrowprops=dict(facecolor='black', width=.5, headwidth = 4., headlength = 4.), fontsize = 14)
+ #axloc.annotate('CNM', (np.log10(CNM_nHmin), np.log10(CNM_Tmax)), xytext=(np.log10(CNM_nHmin)-2., np.log10(CNM_Tmax)+0.5),\
+ # arrowprops=dict(facecolor='black', width=.5, headwidth = 4., headlength = 4.), fontsize = 14)
+ #axloc.annotate('MCs', (np.log10(MC_nHmin), np.log10(MC_Tmin)), xytext=(np.log10(MC_nHmin)-3., np.log10(MC_Tmax)),\
+ # arrowprops=dict(facecolor='black', width=.5, headwidth = 4., headlength = 4.), fontsize = 14)
+
+
+############################################################################
+# make plot
+############################################################################
+def make_plot(outputname):
+
+ fig, ax = plt.subplots(nrows=4, ncols=1, sharex=True, sharey=True, squeeze=True, figsize = (3.5*1.3,8*1.2))
+ fig.subplots_adjust(top = 0.92, bottom = 0.07 , left = 0.13,right = 0.95)
+ fig.suptitle('Runaway collapse zone: examples')
+
+ #############################################################
+ # EAGLE
+ #############################################################
+ make_subplot(ax[0], EAGLE_eps, EAGLE_hmin, EAGLE_mB, 'EAGLE (100 Mpc)$^{3}$')
+ ax[0].set_xlabel('')
+
+ ### #########################################################
+ # EAGLE highres
+ #############################################################
+ make_subplot(ax[1], EAGLE_highres_eps, EAGLE_highres_hmin, EAGLE_highres_mB, 'EAGLE high-res (25 Mpc)$^{3}$')
+ ax[1].set_xlabel('')
+
+ #############################################################
+ # Apostle
+ #############################################################
+ make_subplot(ax[2], APOSTLEL1_eps, APOSTLEL1_hmin, APOSTLEL1_mB, r'APOSTLE L1')
+ ax[2].set_xlabel('')
+
+ ### #########################################################
+ # suggested parameters
+ #############################################################
+ #make_subplot(ax[3], 500., 1.55 * 0.01 * 500., EAGLE_highres_mB, r'Alternative: $\epsilon$ = 500 pc, h$_{\mathrm{min}}$ = 7.75 pc')
+ make_subplot(ax[3], 200., 2., EAGLE_highres_mB, r'Alternative: $\epsilon$ = 200 pc, h$_{\mathrm{min}}$ = 2 pc')
+
+ #############################################################
+ #other options:
+ #############################################################
+
+ #make_subplot(ax[0], TNG100_eps, TNG100_hmin_avgcellsizeSF, TNG100_mB, 'TNG100')
+ #make_subplot(ax[1], APOSTLEL1_eps, APOSTLEL1_hmin, APOSTLEL1_mB, r'APOSTLE L1')
+ #make_subplot(ax[2], APOSTLEL2_eps, APOSTLEL2_hmin, APOSTLEL2_mB, r'APOSTLE L2')
+ #make_subplot(ax[2], TNG50_eps, TNG50_hmin_mincellsize, TNG50_mB, 'TNG50')
+
+
+ fig.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+make_plot(outputname)
diff --git a/data/Springel_Hernquist_temperature_data.dat b/data/Springel_Hernquist_temperature_data.dat
new file mode 100644
index 0000000000000000000000000000000000000000..8f7a66f5a49c336471f64a7703675d6ee7afe31c
--- /dev/null
+++ b/data/Springel_Hernquist_temperature_data.dat
@@ -0,0 +1,35 @@
+#############################################
+# Springel & Hernquist (2003) Fig. 1, top panel
+# digitized, the final value is extrapolated
+# using neff -> 0.76 for the highest densities
+# (middle panel)
+# column 1: log10 baryonic overdensity
+# column 2: log10 effective temperature K
+#############################################
+5.89, 4.00
+5.99, 4.14
+6.10, 4.28
+6.27, 4.49
+6.46, 4.68
+6.67, 4.85
+6.85, 4.98
+7.11, 5.13
+7.36, 5.26
+7.61, 5.38
+7.82, 5.48
+8.09, 5.58
+8.32, 5.67
+8.53, 5.75
+8.81, 5.84
+9.09, 5.91
+9.40, 5.97
+9.67, 5.99
+9.98, 6.00
+10.23, 5.98
+10.49, 5.95
+10.76, 5.90
+11.10, 5.84
+11.35, 5.77
+11.65, 5.71
+11.93, 5.64
+16.00, 4.66
diff --git a/helpers/simulation_plotparameter.py b/helpers/simulation_plotparameter.py
new file mode 100644
index 0000000000000000000000000000000000000000..110a4a59444df0e3c586ebcc713bc5079523cf07
--- /dev/null
+++ b/helpers/simulation_plotparameter.py
@@ -0,0 +1,55 @@
+
+import numpy as np
+
+lognH_min = -4.
+lognH_max = 8.
+dlognH_ticks = 2.
+
+logT_min = 0.8
+logT_max = 4.5
+
+nbins = 128
+
+lognH_arr = np.linspace(lognH_min, lognH_max, nbins)
+logT_arr = np.linspace(logT_min , logT_max, nbins)
+
+lognH_2Darr = np.tile(lognH_arr, (len(logT_arr), 1))
+logT_2Darr = (np.tile(logT_arr, (len(lognH_arr), 1))).T
+
+nH_2Darr = np.power(10., lognH_2Darr)
+T_2Darr = np.power(10., logT_2Darr)
+
+massfrac_min = 0.005
+massfrac_max = 2.
+
+colors = []
+colors.append('#CC0066')
+colors.append('grey')
+colors.append('#0066CC')
+
+color_zone = '#CC0066'
+import matplotlib
+norm = matplotlib.colors.Normalize(vmin=0.0, vmax=8.0)
+cmap = matplotlib.cm.get_cmap('coolwarm')
+color_zone_A = [cmap(norm(4.))]
+color_zone_B = [cmap(norm(3.))]
+color_zone_C = [cmap(norm(7.))]
+color_zone_D = [cmap(norm(6.))]
+color_zone_E = [cmap(norm(1.5))]
+color_zone_F = [cmap(norm(0.2))]
+
+mask_central_kpc = 3.
+
+MaxAgeYoungStarsMyr = 1000.
+r_img_kpc = 25.
+slice_thickness_in_kpc = 10.
+npix = int(4096/4)
+vmin = 5.02
+vmax = 8.48
+vmax_sigma_line = 12.
+cmap_stars = 'bone'
+sigma_min = np.power(10., vmin)
+radialbin_kpc = 0.5
+pixsize_kpc = r_img_kpc / float(npix)
+
+eta_res = 1.2348