diff --git a/Fig6_zones_at_lsmooth_adap.py b/Fig6_zones_at_lsmooth_adap.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a9f6bc5943a1b86e28a8e6ce5ab867a3a6a47f4
--- /dev/null
+++ b/Fig6_zones_at_lsmooth_adap.py
@@ -0,0 +1,233 @@
+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': 'medium',
+ 'axes.titlesize':'small',
+ 'figure.titlesize':'medium',
+ '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 Jeans_length_0
+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 = "Fig6_zones_at_lsmooth_adap.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(0.2, 3.7, r"$\lambda_{\mathrm{J,s}}\approx\lambda_{\mathrm{J,N}} > \mathrm{max}(l_{\mathrm{smooth}}, l_{\mathrm{smooth,min}})$", ha='left', va='bottom',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+ ax.text(0.2, 3.0, r"$\lambda_{\mathrm{J,s}}\approx\lambda_{\mathrm{J,N}} < \mathrm{max}(l_{\mathrm{smooth}}, l_{\mathrm{smooth,min}})$", ha='left', va='top',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+ return
+
+############################################################################
+# Add annotation - lsoftmin =/ lsmoothmin
+############################################################################
+def add_annotation_2(ax):
+ dy = 1.
+ dx = 3.
+ angle = np.rad2deg(np.arctan2(dy, dx))
+ # annotations
+ ax.text(4.0, 5.0, r"$\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(4.0, 4.7, r"$\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(2.5, 3.5, r"$\lambda_{\mathrm{J,s}} > l_{\mathrm{smooth}}$", ha='left', va='bottom',
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor', color = 'black', fontsize = 7)
+ ax.text(2.0, 3.4, r"$\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
+
+
+############################################################################
+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
+
+############################################################################
+# lambda_J,s = l_smooth (or lambda_J,N = l_smooth for Newtonian = True)
+############################################################################
+def add_lambda_Js_equal_lsmooth_adap(ax, mB, lsoftmin, lsmoothmin, fill, ls = 'solid', lc = 'cadetblue', lw = 3., lf = 1., Newtonian = False):
+
+ hmin = lsmoothmin / kernel_support_over_smoothing_length
+ epsmin = lsoftmin / softening_length_over_eps
+
+ kernel_size = kernel_support_over_smoothing_length * h_smooth(mB, nH_2Darr, hmin)
+ if Newtonian:
+ eps = 0.
+ else:
+ eps = np.maximum(kernel_size/softening_length_over_eps, epsmin)
+
+ if fill:
+ 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:
+ 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)
+ return
+
+
+
+ls = 'dashed'
+lc = '#00CC66'
+lw = 3.
+
+fig, ax = plt.subplots(nrows = 1, ncols = 2, width_ratios = [1.0, 1.0], figsize = (7.0, 3.1))
+fig.suptitle("Adaptive l$_{\mathrm{soft}}$: Gravitational stability at the length scale max(l$_{\mathrm{smooth}}$, l$_{\mathrm{smooth,min}}$)", y = 0.92)
+fig.subplots_adjust(bottom = 0.18, top = 0.98, left = 0.07, right = 0.98)
+
+# setup each axis
+setup_axis(ax)
+
+######### Plot nr 1. ########
+
+lsoftmin = 2.25
+lsmoothmin = 2.25
+mB = 6.4e4
+ax[0].set_title(r"m$_{\mathrm{B}}$ = %.1f$\times$10$^4$ M$_{\odot}$, l$_{\mathrm{soft,min}}$ = l$_{\mathrm{smooth,min}}$ = %.1f pc"%(mB/1.e4, lsoftmin))
+add_lambda_Js_equal_lsmooth_adap(ax[0], mB, lsoftmin, lsmoothmin, True, ls = ls, lc = lc, lw = lw, lf = 2., Newtonian = True)
+add_lambda_Js_equal_lsmooth_adap(ax[0], mB, lsoftmin, lsmoothmin, True, ls = 'solid', lc = 'cadetblue', lw = 3.)
+add_lambda_Js_equal_lsmooth_adap(ax[0], mB, lsoftmin, lsmoothmin, False, ls = 'solid', lc = 'cadetblue', lw = 3.)
+add_lambda_Js_equal_lsmooth_adap(ax[0], mB, lsoftmin, lsmoothmin, False, ls = ls, lc = lc, lw = lw, lf = 2., Newtonian = True)
+
+
+add_annotation(ax[0])
+add_ISM_patches(ax[0])
+
+######### Plot nr 2. ########
+lsoftmin = 108.
+lsmoothmin = 0.
+mB = 8.5e4
+ax[1].set_title(r"m$_{\mathrm{B}}$ = %.1f$\times$10$^4$ M$_{\odot}$, l$_{\mathrm{soft,min}}$ = %.1f pc, l$_{\mathrm{smooth,min}}$ = %.1f pc"%(mB/1.e4, lsoftmin, lsmoothmin))
+
+add_lambda_Js_equal_lsmooth_adap(ax[1], mB, lsoftmin, lsmoothmin, True, ls = ls, lc = lc, lw = lw, lf = 2., Newtonian = True)
+add_lambda_Js_equal_lsmooth_adap(ax[1], mB, lsoftmin, lsmoothmin, True, ls = 'solid', lc = 'cadetblue', lw = 3.)
+add_lambda_Js_equal_lsmooth_adap(ax[1], mB, lsoftmin, lsmoothmin, False, ls = 'solid', lc = 'cadetblue', lw = 3.)
+add_lambda_Js_equal_lsmooth_adap(ax[1], mB, lsoftmin, lsmoothmin, False, ls = ls, lc = lc, lw = lw, lf = 2., Newtonian = True)
+
+add_annotation_2(ax[1])
+add_ISM_patches(ax[1])
+
+
+
+plt.tight_layout()
+fig.savefig(outputname, dpi = 250)
+plt.close()
+print ("Figure saved: "+outputname)
+
+
+
+