In [1]:
import yt
import math
from yt import derived_field
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm
from yt.units import pc, kpc, second, Kelvin, gram, erg, cm

%matplotlib inline

# Define some constant parameters to be used.
mp      = 1.6726e-24  * gram # g
mu      = 1.2924
kb      = 1.3806e-16  *erg / Kelvin # erg K-1
GNewton = 6.6743e-8   * cm**3 / (gram * second**2 )# cm3 g-1 s-2
Msun    = 1.9884e33   * gram
mm      = mu*mp

# Create a derived field.
@derived_field(name="numdens", units="1/cm**3", force_override=True)
def numdens(field, data):
    #mp      = 1.6726e-24 * gram
    #mu      = 1.2924
    #mm      = mu*mp
    dens_here = data["dens"].value
    dens_here = dens_here * gram / cm**3
    return dens_here/mm

yt.add_field('numdens', function=numdens, units="1/cm**3", force_override=True)
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/yt/fields/local_fields.py:46: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
  warnings.warn("Because 'sampling_type' not specified, yt will "
In [2]:
# Load the mass vs time calculation and use that to compute the accretion rate.
In [3]:
###############################
# 1pc resolution simulations
###############################

Mdot_1pc = np.array([  1.73222727e-03,   4.07646970e-04,   2.12081818e-04,   3.51164646e-05,
   7.79530303e-05,   1.05671212e-04,   3.16688889e-04,   7.95267677e-04,
   6.71101010e-05,   2.31239899e-03,   8.25696970e-04,   1.55590909e-02,
   2.04323232e-02,   8.64732323e-05,   1.19500000e-02,   1.70525252e-03,
   1.29698485e-03,   3.61572222e-04,   1.21006061e-04,   1.08310404e-04,
   2.24430808e-04,   6.74449495e-05,   3.24201010e-04,   4.22939394e-04,
   5.49661616e-02,   7.37478788e-04,   7.33408081e-04,   3.30342929e-03,
   8.82150505e-04,   4.59412626e-05,   5.89819697e-04,   1.63540303e-04,
   3.91360606e-05,   8.81868687e-06,   7.03035354e-05,   1.06215657e-05,
   5.74717172e-05,   1.86442424e-04,   1.31775758e-04,   5.44204545e-03,
   4.57732323e-04,   2.48601010e-04,   1.25735859e-01,   2.21853535e-04,
   3.69152020e-04,   2.86080808e-05,   7.26136364e-02,   2.60378788e-04,
   6.36363636e-06,   9.00505050e-06,   2.74787879e-03,   1.37448485e-03,
   6.70378788e-02,   3.42112121e-03,   2.91432828e-04,   4.76247475e-03,
   2.97229293e-04,   4.41383838e-05,   5.54771212e-04,   5.54097980e-04,
   2.89565657e-05,   2.39793939e-04,   1.47631313e-04,   2.51873737e-04,
   9.73722222e-03,   2.17712121e-03,   5.17621212e-02,   3.05095960e-04,
   5.07994950e-03,   5.58641414e-05,   3.19590909e-04,   1.34615657e-02])

M_1pc = np.array([  3.10567020e+04,   2.78899747e+03,   2.12153333e+03,   4.01471263e+02,
   1.29298030e+03,   7.75816162e+02,   7.53082323e+02,   5.53081818e+03,
   8.83619192e+02,   3.18578030e+04,   1.74074646e+04,   1.21677020e+06,
   6.85192424e+05,   6.29728788e+02,   1.87407424e+06,   2.93722980e+05,
   6.28638889e+03,   7.63907071e+02,   5.36831818e+02,   3.99252727e+02,
   2.67309091e+03,   5.45559596e+02,   1.06977121e+03,   1.52469798e+04,
   7.05624747e+05,   8.22961616e+02,   1.57744596e+03,   3.61842323e+04,
   8.79026768e+02,   4.58862879e+02,   8.17526263e+02,   4.48073838e+02,
   5.30464646e+02,   1.64824343e+03,   8.86904040e+02,   4.78114646e+02,
   1.30498889e+03,   1.04836768e+03,   1.56430859e+03,   4.91180556e+04,
   1.07681010e+04,   1.24446919e+04,   1.36647424e+06,   7.07491414e+03,
   4.80618636e+03,   9.57207576e+02,   9.99591414e+05,   1.58063384e+04,
   1.43253283e+03,   3.13474040e+04,   5.45601010e+04,   1.26932929e+04,
   8.72909091e+05,   4.91004949e+04,   1.28518990e+03,   3.89434091e+04,
   1.91172273e+03,   1.58234495e+03,   2.40907020e+03,   4.96973030e+03,
   1.02219040e+03,   2.30209293e+03,   5.17533838e+03,   1.60667172e+04,
   3.02350404e+05,   5.73079798e+04,   9.39946970e+05,   6.32778283e+03,
   5.25216162e+04,   6.33141414e+02,   2.09773687e+03,   2.33225859e+05])

Surf_Area = np.array([7.9839733672858331e+40, 1.0145190840820591e+40, 9.5615419599611985e+39, 2.1286018007813093e+39, 
             5.9566518134767256e+39, 3.8452161562501063e+39, 5.3901690761720239e+39, 2.1354682582031833e+40, 
             4.5490280419923129e+39, 9.0431244246100153e+40, 4.4340148801760666e+40, 1.6058927295417408e+42, 
             7.5555064241585968e+41, 3.5533917158204109e+39, 9.5215448454764295e+41, 2.8166208344527795e+41, 
             3.0727396962892016e+40, 3.6563885771485387e+39, 2.9869089785157086e+39, 2.3002632363281887e+39, 
             1.0471347568359661e+40, 2.4719246718750687e+39, 6.5059684072267417e+39, 3.0607233958009442e+40, 
             7.5774790879085835e+41, 3.5705578593750982e+39, 7.5702693076173956e+39, 7.9204586361334755e+40, 
             4.0512098789063619e+39, 2.6092538203125727e+39, 2.6092538203125724e+39, 2.8324136865235163e+39, 
             2.5577553896485088e+39, 7.2956110107423869e+39, 4.6176926162110648e+39, 1.9397742216797414e+39, 
             5.2871722148438955e+39, 5.0640123486329519e+39, 7.2612787236330122e+39, 1.1760524949316846e+41, 
             3.3044826342775303e+40, 2.6487359504883914e+40, 1.5626340477837493e+42, 1.8556601182617694e+40, 
             8.8920623613283678e+39, 4.016877591796986e+39, 4.9874513483779765e+41, 2.4959572728516589e+40, 
             5.3215045019532714e+39, 6.4888022636723975e+40, 1.0533145685156732e+41, 2.2796638640625692e+40, 
             5.8536549521473043e+41, 1.3961224553027576e+41, 5.853654952148599e+39, 5.9463521273441814e+40, 
             1.0059360123047151e+40, 6.849291278320501e+39, 8.9778930791018089e+39, 1.445389287304727e+40, 
             5.6133289423829674e+39, 5.7506580908204712e+39, 7.8620937480470898e+39, 2.0736701414063066e+40, 
             2.1639640565037106e+41, 8.461192158106203e+40, 4.6537415176749751e+41, 2.1114356572266201e+40, 
             1.1973385129394981e+41, 2.986908978515708e+39, 5.836488808593911e+39, 2.1213920204880957e+41])
In [5]:
Surf_Area = Surf_Area / (3.0856e18)**2
In [6]:
#print len(M_1pc)
In [7]:
#print len(Surf_Area)
In [8]:
min_SA = np.min(Surf_Area)
max_SA = np.max(Surf_Area)

big_circle   = 100
small_circle = 10
In [9]:
print(min_SA)
print(max_SA)
203.737949922
168669.780665
In [10]:
# Make circle sizes consistent between the two plots. With and without self-gravity.
min_SA = 182
max_SA = 1.7e5
In [11]:
def calculate_circle_size(Surf_Area):
    
    Surf_Area = np.log10(Surf_Area)
    
    min_SA = np.min(Surf_Area)
    max_SA = np.max(Surf_Area)

    big_circle   = 500
    small_circle = 50
    
    m = (big_circle - small_circle) / (max_SA - min_SA)
    b = small_circle - min_SA * m
    
    circle_size = m * Surf_Area + b
    
    return circle_size
In [12]:
circle_size = calculate_circle_size(Surf_Area)
In [13]:
nism = 1.0
vism = 10.
In [14]:
M_here = np.array([1.0e2, 3.0e6])
#Mdot_G = 8.51e-8 * (nism/1.0)*(M_here/10.)**(1.15)
#Mdot_T = 4.1e-6 * (nism/1.0)*(vism/10.0)*(M_here/10.0)**(0.87)
In [19]:
Mdot_G = 6.0e-6   * (nism/1.0)*            (M_here/1.0e3)**(1.25)
Mdot_T = 1.026e-4 * (nism/1.0)*(vism/10.0)*(M_here/1.0e3)
In [20]:
csize_scale = np.logspace(np.log10(min_SA), np.log10(max_SA), num=5)
scale_sizes = calculate_circle_size(csize_scale)

csx = np.logspace(2.3, 4, num=5)
csy = np.array([3.0e-1, 3.0e-1, 3.0e-1, 3.0e-1, 3.0e-1])

size_text = []
for i in range(len(csize_scale)):
    size_text.append("%.1e" %csize_scale[i])
In [21]:
import matplotlib.patches as patches
In [34]:
cmtopc  = 3.2407557442396e-19
kmtopc  = 3.24078e-14
gtoMsun = 5.0e-34
stoMyr  = 1.0/3.15e13
stoyr  = 1.0/3.15e7


rho_inf = 1.0 * 2.35 * 1.6733e-24 # g cm-3
vinf    = 7.5                     # km s-1

Mass_ZA = np.array([1.0e2, 3.8e4, 6.0e4, 8.0e4]) # Msun

Radius  = np.array([64.0, 64.0, 58., 48.])       # pc

Macc_ZA = 2.0*rho_inf * vinf * np.pi * Radius**2  * (gtoMsun / cmtopc**3 * kmtopc / stoyr)
In [74]:
Masses = np.logspace(np.log10(5.0e2), np.log10(2.0e6))

cmtopc  = 3.2407557442396e-19
gtoMsun = 5.0e-34
stoMyr  = 1.0/3.15e13

Gnewton   = 6.67259E-8 * cmtopc**3 * gtoMsun**(-1)* stoMyr**(-2)
Sigma_low = 8.0
Sigma_high= 16.0
M0        = 5.0e4

#Macc_low  = 4.64 * Gnewton**2 * Sigma_low**3 * np.log10(Masses/M0)**3 
#Macc_high = 4.64 * Gnewton**2 * Sigma_high**3 * np.log10(Masses/M0)**3

Macc_low  = 81.48 * Gnewton**2 * Sigma_low**3 * ((np.log10(Masses) - 4.69)/2.5225225225225215e-02)**3 /1.0e6
Macc_high = 81.48 * Gnewton**2 * Sigma_high**3 * ((np.log10(Masses) - 4.69)/4.4848484848484836e-02)**3/1.0e6
In [83]:
fig  = plt.figure(figsize=(10,9))

#ax   = fig.add_subplot(111)
ax   = fig.add_axes([0.09, 0.06, 0.9, 0.93])

R    = 40.0
mdot = 4*math.pi*(R*3.08e18)**2 * (7*1.0e5) * (10 * mm.value) / 1.98e33 * 3.155e7

Fukui_Mdot = np.array([mdot, mdot])
Fukui_M    = np.array([1.0e4, 3e6])

ax.text(1.0e4, 0.052, "Fukui & Kawamura", fontsize=12)
ax.plot(Fukui_M, Fukui_Mdot, ':k', linewidth=3)

# 1 pc resolution clouds.

ax.plot(Mass_ZA, Macc_ZA, "-b", linewidth=2, label="Zamora-Aviles (2012)")

ax.plot(Masses, Macc_low, "-r", linewidth=2,  label="G11 $\Sigma_{res}=8$  M$_{\odot}$pc$^{-2}$")
ax.plot(Masses, Macc_high,"--r", linewidth=2, label="G11 $\Sigma_{res}=16$M$_{\odot}$pc$^{-2}$")

ax.scatter(M_1pc, Mdot_1pc, s=circle_size, alpha=0.8, c='#1f78b4', edgecolor="k", linewidth='2', label="0.95 pc simulation")

ax.plot(M_here, Mdot_G, "-k", linewidth=2,  label="Grav. Accretion. $\dot{M} \propto M^{\,1.25}$")
#ax.plot(M_here, Mdot_G*2, "-k", linewidth=1  )
#ax.plot(M_here, Mdot_G*4, "-k", linewidth=0.5)
ax.plot(M_here, Mdot_T, "--k", linewidth=2, label="Turb. Accretion. $\dot{M} \propto M$")

ax.scatter(csx, csy, s=scale_sizes, alpha=0.8, c='#1f78b4', edgecolor="k", linewidth='2')
ax.plot([2.0e2, 1.0e4], [2.0e-1, 2.0e-1], "-k", linewidth=2)
for i in range(len(size_text)):
    if i%2 ==0:
        ax.plot([csx[i], csx[i]], [1.9e-1, 2.1e-1], "-k", linewidth=2)
        ax.text(csx[i]*0.7, csy[i]*0.5, size_text[i], rotation=15, fontsize=13)

ax.text(2.1e4, 1.4e-1, "pc$^{2}$", fontsize=14, rotation=15)

ax.add_patch( patches.Rectangle(
        (1.25e2, 1.0e-1),
        3.3e4,
        3.0e-1,
        fill=False      # remove background
    )
)
        
ax.set_xscale("log")
ax.set_yscale("log")

ax.set_xlim(1.0e2, 3.0e6)
ax.set_ylim(1.0e-6, 5.0e-1)

ax.set_ylabel("Mass Accretion rate [M$_{\odot}$ yr$^{-1}$]", fontsize=15)
ax.set_xlabel("Cloud Mass [M$_{\odot}$]", fontsize=15)

ax.xaxis.set_tick_params(which="minor", width=1.5, length=3.0, direction="in")
ax.xaxis.set_tick_params(which="major", width=1.5, length=6.0, direction="in")

ax.yaxis.set_tick_params(which="minor", width=1.5, length=3.0, direction="in")
ax.yaxis.set_tick_params(which="major", width=1.5, length=6.0, direction="in")

ax.legend(loc=4, fontsize=15, framealpha=1.0, fancybox=True)

for tick in ax.xaxis.get_major_ticks():
    tick.label.set_fontsize(15) 
for tick in ax.yaxis.get_major_ticks():
    tick.label.set_fontsize(15) 

fig.show()

#save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/"
save_dir = "/home/jcibanezm/Dropbox/Projects/Papers/Submitted/AccretionPaper/Figures/"

#fig.savefig(save_dir + "Mass_Acc_Cloud_Pop_SG.pdf", format='pdf')
/home/jcibanezm/codes/libs/yt3.4/yt-conda/lib/python3.6/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [ ]: