Added Fig 6

Fig6_zones_at_lsmooth_adap.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': '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)
+
+
+
+