Added Fig 3

Fig3_zone_of_avoidance.py

+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
+
+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
+
+
+eta_res = 1.2348
+XH      = 0.74
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig3_zone_of_avoidance.png"
+
+############################################################################
+# set resolution properties to plot
+# epsarr ... gravitational softening scale in pc
+# mBarr  ... gas particle mass in Msol
+# h_min_ratio_arr ... minimum smoothing length,
+#                     defined as h_min = 1.55 * h_min_ratio * eps
+############################################################################
+epsarr  = np.asarray([31.6, 100., 316., 700., 1000.])
+mBarr   = np.asarray([1.81e5, 1.81e6, 1.81e6, 1.81e7])
+h_min_ratio_arr = np.asarray([0.01, 0.0316, 0.1, 0.316])
+
+lwfat = 4
+lwthin = 2
+
+############################################################################
+# equations from paper
+############################################################################
+# mB in units of solar mass
+# hmin in units of pc
+def zone_min_density(mB, hmin):
+    return 58. * (mB / 1.e6) * np.power(hmin / 100., -3)
+
+# hmin in units of pc
+# eps in units of pc
+# nH in units of cm-3
+def zone_max_temperature(hmin, eps, nH):
+    return 27. * np.power(hmin / 100., 5) * np.power(eps / 700., -3) * (nH / 100.)   
+
+############################################################################
+# set plot limits and define 2D arrays for contour plots
+############################################################################
+idefault_eps  = 3
+idefault_mB   = 2
+idefault_hmin = 2
+
+ncolors = np.maximum(len(epsarr), len(mBarr))
+ncolors = np.maximum(ncolors, len(h_min_ratio_arr))
+c       = pl.cm.copper(np.linspace(0,0.9,ncolors))
+ls= ['dotted', 'dashed', 'dotted', 'dashed', 'dotted', 'dashed', 'dotted', 'dashed']
+
+lognH_min    = -6.
+lognH_max    = 8.
+dlognH       = 0.01
+
+logT_min      = 1.
+logT_max      = 6.
+dlogT         = 0.1 
+
+xticks = np.arange(-6., 10., 2)
+
+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]]
+
+def logTeff_EAGLE(lognH):
+    return 1./3. * (lognH +1.) + np.log10(8000.)
+
+############################################################################
+# make plot
+############################################################################
+def make_plot(outputname):
+
+    fig, ax = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=True, squeeze=True, figsize = (3*1.3,7.4))
+    fig.subplots_adjust(top = 0.915, bottom = 0.07 , left = 0.15, right = 0.95)
+    fig.suptitle('Runaway collapse zone')
+    
+    for axloc in ax:
+        axloc.set_xlim(lognH_min, lognH_max)
+        axloc.set_xlim(logT_min, logT_max)
+
+    ################### Plot 0 ############################
+    ax[0].set_title('Vary baryon particle mass')
+    #######################################################
+    ax[0].set_ylabel('log T [K]')
+    ax[0].xaxis.set_ticks(xticks)
+    ax[0].set_ylim(logT_min, logT_max)
+    ieps  = idefault_eps
+    ihmin = idefault_hmin
+    h = []
+    l = []
+    for imB in range(len(mBarr)):
+        eps  = epsarr[ieps]
+        h_min = 1.55 * eps * h_min_ratio_arr[ihmin]
+        mB   = mBarr[imB] 
+        
+
+        if imB == idefault_mB:
+            currentcolor = "black"
+            currentlinestyle = "solid"
+            lw = lwfat
+        else:  
+            currentcolor = c[imB]
+            currentlinestyle = ls[imB]
+            lw = lwthin
+
+
+        ax[0].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+                   np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+                   lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '10$^{%.1f}$'%((np.log10(mBarr[imB]))))
+
+        ax[0].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+                     lw = lw, color = currentcolor, linestyle = currentlinestyle)
+      
+ 
+    # add EAGLE EOS:       
+    ax[0].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+    # with rotated annotation:
+    x = [ 0., 4. ]
+    y = [logTeff_EAGLE(0.), logTeff_EAGLE(4.)]
+    dy = y[1] - y[0]
+    dx = x[1] - x[0]
+    angle = np.rad2deg(np.arctan2(dy, dx))
+
+    ax[0].text(x[0], y[0]+0.1, "EAGLE T$_{\mathrm{eff}}$", ha = 'left', va = 'bottom',\
+               transform_rotates_text=True, rotation=angle, rotation_mode='anchor')    
+
+    ax[0].text(0.05,0.05, '$\epsilon$ = %i pc'%(round(eps))+'\n'+\
+               'h$_{\mathrm{min}}$ = %.1f pc'%(h_min),\
+                va = 'bottom', ha = 'left', transform=ax[0].transAxes)    
+    ax[0].text(0.97, 0.05, "runaway\n collapse zone", va = 'bottom', ha = 'right', transform=ax[0].transAxes)
+    ax[0].legend(loc = 'upper left', title = 'm$_{\mathrm{B}}$ [M$_{\odot}$]')
+
+    ################### Plot 1 ############################
+    ax[1].set_title('Vary softening scale')
+    #######################################################
+    ax[1].set_ylabel('log T [K]')
+    ax[1].xaxis.set_ticks(xticks)
+    imB   = idefault_mB
+    ihmin = idefault_hmin
+    h = []
+    l = []
+    for ieps in range(len(epsarr)):
+        eps  = epsarr[ieps]
+        h_min = 1.55 * epsarr[idefault_eps] * h_min_ratio_arr[idefault_hmin]
+        mB   = mBarr[imB] 
+        
+        if ieps == idefault_eps:
+            currentcolor = "black"
+            currentlinestyle = "solid"
+            lw = lwfat
+        else:  
+            currentcolor = c[ieps]
+            currentlinestyle = ls[ieps]    
+            lw = lwthin
+     
+        ax[1].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+                   np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+                   lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '%i'%(round(epsarr[ieps])))
+
+        ax[1].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+                     lw = lw, color = currentcolor, linestyle = currentlinestyle)
+
+    # add EAGLE EOS:       
+    ax[1].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+
+    ax[1].text(0.05,0.05, 'm$_{\mathrm{B}}$ = 10$^{%.1f}$ M$_{\odot}$'%(np.log10(mB))+'\n'+\
+               'h$_{\mathrm{min}}$ = %.1f pc'%(h_min),\
+                va = 'bottom', ha = 'left', transform=ax[1].transAxes)        
+    ax[1].legend(loc = 'upper left', title = '$\epsilon$ [pc]')
+
+
+    ################### Plot 2 ############################
+    ax[2].set_title('Vary minimum smoothing length')
+    #######################################################
+    ax[2].set_ylabel('log T [K]')
+    ax[2].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+    ax[2].xaxis.set_ticks(xticks)
+    imB  = idefault_mB
+    ieps = idefault_eps
+    h = []
+    l = []
+    for ihmin in range(len(h_min_ratio_arr)):
+        eps  = epsarr[ieps]
+        h_min = 1.55 * eps * h_min_ratio_arr[ihmin]
+        mB   = mBarr[imB] 
+        
+        if ihmin == idefault_hmin:
+            currentcolor = "black"
+            currentlinestyle = "solid"
+            lw = lwfat
+        else:  
+            currentcolor = c[ihmin]
+            currentlinestyle = ls[ihmin]    
+            lw = lwthin
+        
+        ax[2].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+                   np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+                   lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '%.1f'%(h_min))
+
+        ax[2].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+                     lw = lw, color = currentcolor, linestyle = currentlinestyle)
+
+    h, l = ax[2].get_legend_handles_labels()
+
+    # add EAGLE EOS:       
+    ax[2].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+
+    ax[2].text(0.05,0.05, 'm$_{\mathrm{B}}$ = 10$^{%.1f}$ M$_{\odot}$'%(np.log10(mB))+'\n' +\
+               '$\epsilon$ = %i pc'%(round(eps)),\
+                va = 'bottom', ha = 'left', transform=ax[2].transAxes)    
+    ax[2].legend(h[::-1],l[::-1], loc = 'upper left', title = 'h$_{\mathrm{min}}$ [pc]')
+
+    fig.savefig(outputname, dpi = 150)
+    plt.close()
+    print ("Figure saved: "+outputname)
+
+if __name__ == "__main__":
+    make_plot(outputname)