Skip to content
Snippets Groups Projects
Commit b6d0dc0e authored by Sylvia Plöckinger's avatar Sylvia Plöckinger
Browse files

Added Fig 1

parent 8e37f7c9
No related branches found
No related tags found
No related merge requests found
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib.pylab as pylab
params = {'legend.fontsize': 'large',
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
'xtick.labelsize':'large',
'ytick.labelsize':'large'}
pylab.rcParams.update(params)
import sys
sys.path.insert(0, "helpers/")
from equations import Jeans_length_0
from equations import Jeans_mass_0
from equations import Jeans_mass_P
from equations import Jeans_mass_W
############################################################################
# Set Plummer softening scale in pc
############################################################################
eps = 100. # in pc
H = 3. * eps
############################################################################
# Set output filename
############################################################################
outputname = 'Fig1_Jeans_mass_softening_%4.4ipc.png'%(round(eps))
############################################################################
# set plot limits and define 2D arrays for contour plots
############################################################################
lognH_min = -6.
lognH_max = 6.
dlognH = 0.1
logT_min = 1.
logT_max = 6.
dlogT = 0.1
LJmin = np.log10(1.)
LJmax = np.log10(10. * 1000.)
MJmin = np.log10(1.e1)
MJmax = np.log10(1.e10)
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]]
# manual locations for contour labels
manual_locations = [ (4.5,1.5), \
(3.0,1.5), \
(1.0,1.5), \
(-1.,1.5), \
(-2.0,1.5), \
(-4.5,1.5), \
(-2.0,3.), \
(-4.5,3.), \
(-2.0,4.5), \
(-4.5,4.) ]
ls = ['solid', 'dotted']
lines = [Line2D([0], [0], color='black', linewidth=2, linestyle=l) for l in ls]
labels = ['Plummer', 'Wendland C2']
linesN = [Line2D([0], [0], color='black', linewidth=2, linestyle='solid')]
labelsN = ['Newton']
linesL = [Line2D([0], [0], color='#CC0066', linewidth=2, linestyle='dashed')]
labelsL = ['$\lambda_{\mathrm{J,N}} \equiv \epsilon = %i$ pc'%(round(eps))]
############################################################################
# get Jeans masses and Newtonian Jeans lengths for full phase-space
############################################################################
logMJ_0_2D = np.log10(Jeans_mass_0(nH_2Darr, T_2Darr))
logMJ_P_2D = np.log10(Jeans_mass_P(nH_2Darr, T_2Darr, eps))
logMJ_W_2D = np.log10(Jeans_mass_W(nH_2Darr, T_2Darr, H))
logLJ_0_2D = np.log10(Jeans_length_0(nH_2Darr, T_2Darr))
############################################################################
# make plot
############################################################################
def make_plot(outputname):
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, squeeze=True, figsize = (5*1.1,4*1.1))
fig.subplots_adjust(top = 0.9, bottom = 0.15)
fig.suptitle('Jeans mass (log M$_{\mathrm{J}}$ [M$_{\odot}$])', fontsize = 14)
ax[0].set_ylabel('log T [K]')
CF_0 = ax[0].contour(lognH_2Darr, logT_2Darr, logMJ_0_2D, levels = np.arange(MJmin, MJmax + 1, 1), colors = 'black')
ax[0].clabel(CF_0, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations)
ax[0].legend(linesN,labelsN, loc = 'upper left', framealpha = 1.0)
ax[1].set_ylabel('log T [K]')
ax[1].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
CF_P = ax[1].contour(lognH_2Darr, logT_2Darr, logMJ_P_2D, levels = np.arange(MJmin, MJmax + 1, 1), colors = 'black')
CF_W = ax[1].contour(lognH_2Darr, logT_2Darr, logMJ_W_2D, levels = np.arange(MJmin, MJmax + 1, 1), linestyles = 'dotted', colors = 'black')
ax[1].clabel(CF_P, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations[5:], colors = 'white')
ax[1].clabel(CF_W, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations[5:])
legend1 = ax[1].legend(lines,labels, loc = 'upper left', framealpha = 1.0)
CF_extra = ax[1].contour(lognH_2Darr, logT_2Darr, logLJ_0_2D-np.log10(eps), levels = np.asarray([0.]), colors ='#CC0066', linestyles = 'dashed')
legend2 = ax[1].legend(linesL,labelsL, loc = 'lower right', framealpha = 1.0)
ax[1].add_artist(legend1)
ax[1].add_artist(legend2)
fig.savefig(outputname, dpi = 150)
plt.close()
print ("Figure saved: "+outputname)
make_plot(outputname)
import numpy as np
from unyt import msun, pc
from unyt import proton_mass_cgs as mH
XH = 0.74
eta_res = 1.2348
############################################################################
# nH ... gas density in cm-3
# T ... gas temperature in K
# returns the Newtonian Jeans length in pc
############################################################################
def Jeans_length_0(nH,T):
return 5.4 * np.power(T, 1./2.) * np.power(nH, -1./2.)
############################################################################
# nH ... gas density in cm-3
# T ... gas temperature in K
# returns the Newtonian Jeans mass in Msol
############################################################################
def Jeans_mass_0(nH, T):
return 21.2 * np.power(T, 3./2.) * np.power(nH, -1./2.)
############################################################################
# nH ... gas density in cm-3
# T ... gas temperature in K
# returns the Plummer softened Jeans mass in Msol
############################################################################
def Jeans_mass_P(nH, T, eps):
return Jeans_mass_0(nH, T) * np.power(1. + 1.42 * np.power(eps/Jeans_length_0(nH,T), 3./2.), 6./5.)
############################################################################
# nH ... gas density in cm-3
# T ... gas temperature in K
# returns the softened Jeans length (for Wendland C2) in pc
############################################################################
def Jeans_length_W(nH,T,eps):
H = 3. * eps
return Jeans_length_0(nH,T) * np.power(1. + 0.27 * np.power(H / Jeans_length_0(nH,T), 2.), 0.3)
############################################################################
# nH ... gas density in cm-3
# T ... gas temperature in K
# returns the Wendland C2 softened Jeans mass in Msol
############################################################################
def Jeans_mass_W(nH, T, H):
return Jeans_mass_0(nH, T) * np.power(1. + 0.27 * np.power(H/Jeans_length_0(nH,T), 2.), 9./10.)
############################################################################
# mB ... particle mass resolution in Msol
# hmin ... minimum smoothing length in pc
# returns the density above which the smoothing length is limited by hmin
############################################################################
def density_critical_hmin(mB, hmin):
mBsol = mB * msun
hminpc = hmin * pc
nH = XH / mH * np.power(eta_res,3) * mBsol / np.power(hminpc, 3.)
nH.convert_to_units('cm**-3')
return nH
############################################################################
# mB ... particle mass resolution in Msol
# nH ... gas density in cm-3
# hmin ... minimum smoothing length in pc
# returns the smoothing length in pc
############################################################################
def h_smooth(mB, nH, hmin):
h_pc = 82. * np.power(mB / 1.e6, 1./3.) * np.power(nH / 100., -1./3.)
h_pc = np.maximum(h_pc, hmin)
return h_pc
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment