diff --git a/Fig5_zones_at_lsmooth.py b/Fig5_zones_at_lsmooth.py
new file mode 100644
index 0000000000000000000000000000000000000000..3bc7e3da09f13a73cd70313eafea89314f399bd4
--- /dev/null
+++ b/Fig5_zones_at_lsmooth.py
@@ -0,0 +1,247 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from unyt import msun, pc, cm, K
+from unyt import proton_mass_cgs as mH
+from unyt import gravitational_constant_cgs as G
+from unyt import boltzmann_constant_cgs as kB
+import matplotlib.patches as mpatches
+from matplotlib.collections import PatchCollection
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'small',
+ 'axes.labelsize': 'large',
+ 'axes.titlesize':'medium',
+ 'figure.titlesize':'large',
+ 'xtick.labelsize':'medium',
+ 'ytick.labelsize':'medium'}
+pylab.rcParams.update(params)
+from matplotlib.lines import Line2D
+
+
+import sys
+sys.path.insert(0, "helpers/")
+from equations import Jeans_length_W
+from equations import h_smooth
+
+eta_res = 1.2348
+XH = 0.74
+gamma_hydro = 5./3.
+
+kernel_support_over_smoothing_length = 1.936492
+softening_length_over_eps = 1.5
+gamma_kernel = kernel_support_over_smoothing_length
+
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig5_zones_at_lsmooth.png"
+
+
+lognH_min = -4.
+lognH_max = 8.
+dlognH = 0.01
+
+logT_min = 1.
+logT_max = 7.
+dlogT = 0.01
+
+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]]
+
+############################################################################
+# 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.
+
+############################################################################
+# Add annotation
+############################################################################
+def add_annotation(ax):
+ dy = 1.
+ dx = 3.
+ angle = np.rad2deg(np.arctan2(dy, dx))
+ # annotations
+ ax.text(-3.5, 2.5, "gravitationally stable in nature at l$_{\mathrm{smooth}}$ ($\lambda_{\mathrm{J,N}} > l_{\mathrm{smooth}}$)", ha='left', va='bottom',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+ ax.text(-3.5, 2.3, "gravitationally unstable in nature at l$_{\mathrm{smooth}}$ ($\lambda_{\mathrm{J,N}} < l_{\mathrm{smooth}}$)", ha='left', va='top',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+
+ dy = -2.
+ dx = 3.
+ angle = np.rad2deg(np.arctan2(dy, dx))
+ ax.text(3.8, 3.5, "gravitationally stable in sim\nat l$_{\mathrm{smooth}}$ ($\lambda_{\mathrm{J,s}} > l_{\mathrm{smooth}}$)", ha='left', va='top',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+ ax.text(2.0, 3.7, "gravitationally unstable in sim\nat l$_{\mathrm{smooth}}$ ($\lambda_{\mathrm{J,s}} < l_{\mathrm{smooth}}$)", ha='left', va='top',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+
+ return
+
+############################################################################
+# lambda_J,s = l_smooth (or lambda_J,N = l_smooth for lsoft = 0)
+############################################################################
+def add_lambdaJs_equal_lsmooth(ax, mB, lsoft, lsoftdefault, fill, ls = 'solid', lc = 'cadetblue', lw = 3., lf = 1.):
+
+ if np.isscalar(mB):
+ kernel_size = kernel_support_over_smoothing_length * h_smooth(mB, nH_2Darr, 0.0)
+ else:
+ print ("add_lambdaJs_equal_lsmooth not implemented for multiple values of mB")
+ return
+
+ if fill:
+ if np.isscalar(lsoft):
+ eps = lsoft/softening_length_over_eps
+ else:
+ eps = lsoft[lsoftdefault]/softening_length_over_eps
+
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps) < kernel_size) ] = lf
+ all_zones[all_zones == 0] = np.nan
+ im = ax.imshow(all_zones, extent = [lognH_2Darr.min(), lognH_2Darr.max(), logT_2Darr.min(), logT_2Darr.max()], \
+ origin = 'lower', alpha = 0.3, vmin = 0, vmax = 3, aspect = 'auto', cmap = 'viridis')
+
+ else:
+ # add line for default softening length
+ if np.isscalar(lsoft):
+ eps = lsoft/softening_length_over_eps
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[:] = Jeans_length_W(nH_2Darr, T_2Darr, eps) / kernel_size
+ CS = ax.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(1.)], colors = lc, linewidths = lw, linestyles = ls)
+
+ # add lines for other softening values
+ if not np.isscalar(lsoft):
+ for l in lsoft:
+ eps = l/softening_length_over_eps
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[:] = Jeans_length_W(nH_2Darr, T_2Darr, eps) / kernel_size
+ CS = ax.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(1.)], colors = lc, linewidths = lw, linestyles = ls)
+ ax.clabel(CS, CS.levels, inline=True, fmt="%i pc"%(int(l)), fontsize=10, colors = 'black')
+
+ return
+
+############################################################################
+def add_ISM_patches(ax):
+ patchcolor = 'black'
+ 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 = patchcolor, edgecolor = patchcolor, alpha = 0.3, \
+ zorder = 100)
+ 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 = patchcolor, edgecolor = patchcolor, alpha = 0.3, \
+ zorder = 100)
+ 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 = patchcolor, edgecolor = patchcolor, alpha = 0.3, \
+ zorder = 100)
+ ax.add_patch(WNM)
+ ax.add_patch(CNM)
+ ax.add_patch(MCs)
+ return
+
+############################################################################
+def setup_axis(axarr):
+ axarr[0].set_ylabel('log T [K]')
+ for ax in axarr:
+ ax.set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ ax.yaxis.set_ticks(yticks)
+ ax.xaxis.set_ticks(xticks)
+ ax.set_xlim(lognH_min, lognH_max)
+ ax.set_ylim(logT_min, logT_max)
+ return
+
+############################################################################
+mB = [1.e5/8., 1.e5, 8.e5] # msun
+mB_default = 1
+lsoft = [10., 20., 50., 100., 200., 500.] # pc
+lsoft_default = 3
+hmin = 0.2 / kernel_support_over_smoothing_length
+
+ls = 'dashed'
+lc = '#00CC66'
+lw = 3.
+
+fig, ax = plt.subplots(nrows = 1, ncols = 3, width_ratios = [0.9, 0.5, 0.5], figsize = (9.0, 3.7))
+fig.suptitle("Constant l$_{\mathrm{soft}}$: Gravitational stability at the hydro spatial resolution limit l$_{\mathrm{smooth}}$", y = 0.92)
+fig.subplots_adjust(bottom = 0.18, top = 0.98, left = 0.07, right = 0.98)
+
+# setup each axis
+setup_axis(ax)
+ax[0].set_title("m$_{\mathrm{B}}$ = 10$^5$ M$_{\odot}$, l$_{\mathrm{soft}}$ = 100 pc")
+ax[1].set_title("m$_{\mathrm{B}}$ = 10$^5$ M$_{\odot}$")
+ax[2].set_title("l$_{\mathrm{soft}}$ = 100 pc")
+
+######### Plot nr 1. ########
+
+add_lambdaJs_equal_lsmooth(ax[0], mB[mB_default], 0., 0, True, ls, lc, lw, lf = 2)
+add_lambdaJs_equal_lsmooth(ax[0], mB[mB_default], lsoft[lsoft_default], lsoft_default, True)
+add_lambdaJs_equal_lsmooth(ax[0], mB[mB_default], lsoft[lsoft_default], lsoft_default, False)
+add_lambdaJs_equal_lsmooth(ax[0], mB[mB_default], 0., 0, False, ls, lc, lw)
+add_annotation(ax[0])
+
+######### Plot nr 2. ########
+
+add_lambdaJs_equal_lsmooth(ax[1], mB[mB_default], 0., 0, True, ls, lc, lw, lf = 2)
+add_lambdaJs_equal_lsmooth(ax[1], mB[mB_default], lsoft[lsoft_default], lsoft_default, True)
+add_lambdaJs_equal_lsmooth(ax[1], mB[mB_default], lsoft, lsoft_default, False, 'solid', 'black', 1.)
+add_lambdaJs_equal_lsmooth(ax[1], mB[mB_default], 0., 0, False, ls, lc, lw)
+
+add_ISM_patches(ax[1])
+ax[1].text(0.02, 0.98, "solid line labels: l$_{\mathrm{soft}}$", transform = ax[1].transAxes, ha = 'left', va = 'top')
+
+######### Plot nr 3. ########
+add_lambdaJs_equal_lsmooth(ax[2], mB[mB_default], 0., 0, True, ls, lc, lw, lf = 2)
+add_lambdaJs_equal_lsmooth(ax[2], mB[mB_default], lsoft[lsoft_default], lsoft_default, True)
+
+# individuall add some more lines
+add_lambdaJs_equal_lsmooth(ax[2], mB[0], lsoft[lsoft_default], lsoft_default, False, 'dotted', 'black', 1.5)
+add_lambdaJs_equal_lsmooth(ax[2], mB[mB_default], lsoft[lsoft_default], lsoft_default, False, 'dashed', 'black', 1.5)
+add_lambdaJs_equal_lsmooth(ax[2], mB[2], lsoft[lsoft_default], lsoft_default, False, 'dashdot', 'black', 1.5)
+add_lambdaJs_equal_lsmooth(ax[2], mB[0], 0., 0, False, 'dotted', '#00CC66', 3.)
+add_lambdaJs_equal_lsmooth(ax[2], mB[mB_default], 0., 0, False, ls, lc, lw)
+add_lambdaJs_equal_lsmooth(ax[2], mB[2], 0., 0, False, 'dashdot', '#00CC66', 3.)
+
+add_ISM_patches(ax[2])
+
+#manually add legend
+ls = ['dotted', 'dashed', 'dashdot']
+lines = [Line2D([0], [0], color = 'black', linewidth = 2, linestyle = s) for s in ls]
+labels = [r"$10^5\,/\,8$", r"$10^5$", r"$10^5 \times 8$"]
+ax[2].legend(lines[::-1], labels[::-1], loc = 'upper left', title = 'm$_{\mathrm{B}}$ [M$_{\odot}$]', handlelength=3)
+
+plt.tight_layout()
+fig.savefig(outputname, dpi = 250)
+plt.close()
+print ("Figure saved: "+outputname)
+
+
+
+