Added Fig 5

Fig5_zones_at_lsmooth.py

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