diff --git a/Fig7_zones_below_lsmooth.py b/Fig7_zones_below_lsmooth.py
new file mode 100644
index 0000000000000000000000000000000000000000..51ea17dc3917e1b29115503cd30d67f64ec5cf6e
--- /dev/null
+++ b/Fig7_zones_below_lsmooth.py
@@ -0,0 +1,218 @@
+
+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': 'small',
+ 'axes.labelsize': 'large',
+ 'axes.titlesize':'large',
+ 'figure.titlesize':'large',
+ 'xtick.labelsize':'medium',
+ 'ytick.labelsize':'medium'}
+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 equations import h_smooth
+
+eta_res = 1.2348
+XH = 0.74
+
+kernel_support_over_smoothing_length = 1.936492
+softening_length_over_eps = 1.5
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig7_zones_below_lsmooth.png"
+
+
+lognH_min = -4.
+lognH_max = 8.
+dlognH = 0.01
+
+logT_min = 1.
+logT_max = 6.
+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.
+
+############################################################################
+# make subplot:
+# axloc ... matplotlib axes object
+# eps ... gravitational softening scale in pc
+# h_min ... minimum smoothing length in pc
+# mB ... gas particle mass resolution
+############################################################################
+def make_subplot(axloc, h_min, mB, epsarr):
+
+ 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)
+
+ kernel_size = kernel_support_over_smoothing_length * h_smooth(mB, nH_2Darr, h_min)
+ kernel_size_min = kernel_support_over_smoothing_length * h_min
+
+ # Gravitational instabilities in nature:
+ #all_zones = np.zeros_like(lognH_2Darr)
+ #all_zones[ (Jeans_length_0(nH_2Darr, T_2Darr) < kernel_size) ] = 2.
+ #all_zones[all_zones == 0] = np.nan
+ #im = axloc.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')
+
+ # Gravitational instabilities in simulations
+ if len(epsarr) == 1:
+ eps = epsarr[0]
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) < kernel_size) ] = 1.
+ all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) < kernel_size) & \
+ (kernel_size > eps) ] = 0.2
+ all_zones[all_zones == 0] = np.nan
+ im = axloc.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')
+
+ axloc.text(0.02, 0.02, "l$_{\mathrm{soft}}$ = %i pc"%(eps), transform = axloc.transAxes, ha = 'left', va = 'bottom')
+
+
+ # Add contours for various softening lengths
+ for eps in epsarr:
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[:] = Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) / kernel_size
+ if len(epsarr) > 1:
+ CS = axloc.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(1.)], colors = 'black')
+ axloc.clabel(CS, CS.levels, inline=True, fmt="%i pc"%(int(eps)), fontsize=10)
+ else:
+ CS = axloc.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(1.)], colors = 'cadetblue', linewidths = 3.)
+
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) < kernel_size) & \
+ (kernel_size > eps) ] = 0.2
+ all_zones[all_zones == 0] = np.nan
+ im = axloc.imshow(all_zones, extent = [lognH_2Darr.min(), lognH_2Darr.max(), logT_2Darr.min(), logT_2Darr.max()], \
+ origin = 'lower', alpha = 0.1, vmin = 0, vmax = 3, aspect = 'auto', cmap = 'viridis')
+ if len(epsarr) > 1:
+ all_zones = np.zeros_like(lognH_2Darr)
+ all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) < kernel_size) & \
+ (kernel_size > eps) ] = 0.2
+ CS = axloc.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(0.19)], colors = 'purple', linewidths = 2., linestyles = 'dotted')
+
+ #eps = 0.
+ #all_zones = np.zeros_like(lognH_2Darr)
+ #all_zones[:] = Jeans_length_W(nH_2Darr, T_2Darr, eps/softening_length_over_eps) / kernel_size
+ #CS = axloc.contour(lognH_2Darr, logT_2Darr, all_zones, levels = [(1.)], colors = '#00CC66', linestyles = 'dashed', linewidths = 3.)
+
+ # add the hmin line
+ #axloc.axvline(np.log10(density_critical_hmin(mB, h_min)), linestyle = 'dotted', color = 'black', ymin = 0., ymax = 0.2)
+
+ # add the patches for the gas phases
+ 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)
+ if len(epsarr) > 1:
+ axloc.add_patch(WNM)
+ axloc.add_patch(CNM)
+ axloc.add_patch(MCs)
+
+ for eps in epsarr:
+ nH = np.power(kernel_support_over_smoothing_length * eta_res / (eps * pc), 3) * XH * mB * msun/ mH
+ nH.convert_to_cgs()
+
+ if len(epsarr) == 1:
+ eps = epsarr[0]
+ nH = np.power(kernel_support_over_smoothing_length * eta_res / (eps * pc), 3) * XH * mB * msun/ mH
+ nH.convert_to_cgs()
+ axloc.axvline(np.log10(nH), color = 'purple', linestyle = 'dotted', linewidth = 2.)
+
+ axloc.text(np.log10(nH)-0.3, 5.8, "$l_{\mathrm{soft}} < l_{\mathrm{smooth}}$", ha = 'right', va = 'top')
+ axloc.text(np.log10(nH)+0.3, 5.8, "$l_{\mathrm{soft}} > l_{\mathrm{smooth}}$", ha = 'left', va = 'top')
+ return im
+
+def get_title(eps, hmin):
+ if eps < 0.:
+ title = 'l$_{\mathrm{soft,min}}$ = l$_{\mathrm{smooth,min}}$ = %i pc'%(round(hmin * kernel_support_over_smoothing_length))
+ else:
+ title = 'l$_{\mathrm{soft}}$ = %i pc, l$_{\mathrm{smooth,min}}$ = %i pc'%(round(eps * softening_length_over_eps), round(hmin * kernel_support_over_smoothing_length))
+
+ return title
+
+def make_one_plot(outputname):
+ mB = 1.e5
+ fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (3.5*1.3*1.8,3.5))
+ fig.suptitle("Gravitational stability below the hydro spatial resolution limit l$_{\mathrm{smooth}}$")
+ fig.subplots_adjust(bottom = 0.18, top = 0.8, left = 0.07, right = 0.98)
+
+ title = 'Constant softening, m$_{\mathrm{B}}$ = 10$^5$ M$_{\odot}$, l$_{\mathrm{smooth,min}}$ = 0pc'
+ plt.figtext(0.5,0.84,title, va="bottom", ha="center")
+
+ hmin = 0.2 / kernel_support_over_smoothing_length
+
+ epsarr = [100.]
+ im = make_subplot(ax[0], hmin, mB, epsarr)
+
+ epsarr = [10., 20., 50., 100., 200., 500.]
+ im = make_subplot(ax[1], hmin, mB, epsarr)
+
+ fig.savefig(outputname, dpi = 250)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+
+make_one_plot(outputname)