In [1]:
import yt
import os
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

import Read_binary as Rd
import Vector_computations as vec
import Locate_Nearby_SN as NSN

import copy
from matplotlib.colors import LogNorm

%matplotlib inline

#import Read_particles as Rd

ppc = 3.0856e18

mp      = 1.6726e-24# g
mu      = 1.2924
kb      = 1.3806e-16   # erg K-1
GNewton = 6.6743e-8   # cm3 g-1 s-2
Msun    = 1.9884e33   # gram
mm      = mu*mp
In [2]:
Cloud_name = "M3e3"


if Cloud_name == "M4e3":
    resolution = 0.06
#    resolution = 0.5
    tinit  = 0.0
    tfin   = 9.7

elif Cloud_name == "M8e3":
    resolution = 0.1
    #resolution = 0.5
    tinit  = 0.0
    tfin   = 9.8
elif Cloud_name == "M3e3":
    resolution = 0.06
    tinit  = 0.0
    tfin   = 9.8


if resolution < 0.1:
    resolution_str = "%.3i" %(resolution*100)
else:
    resolution_str = "%.2i" %(resolution*10)
In [3]:
pp_test = Rd.Read_One_Particle_Snapshot(Cloud_name, resolution, 40)
Reading the File M3e3_Particles_snp040.dat
In [4]:
particle_one = Rd.Read_One_InstAcc(Cloud_name, resolution, 40)
Reading the File M3e3_InstAccParts_snpsht40.dat
In [5]:
particle_evol = Rd.Read_InstAcc_Evolution(Cloud_name, resolution, tinit, tfin)
Reading the File M3e3_InstAccParts_snpsht00.dat
Reading the File M3e3_InstAccParts_snpsht00.dat
Reading the File M3e3_InstAccParts_snpsht01.dat
Reading the File M3e3_InstAccParts_snpsht02.dat
Reading the File M3e3_InstAccParts_snpsht03.dat
Reading the File M3e3_InstAccParts_snpsht04.dat
Reading the File M3e3_InstAccParts_snpsht05.dat
Reading the File M3e3_InstAccParts_snpsht06.dat
Reading the File M3e3_InstAccParts_snpsht07.dat
Reading the File M3e3_InstAccParts_snpsht08.dat
Reading the File M3e3_InstAccParts_snpsht09.dat
Reading the File M3e3_InstAccParts_snpsht10.dat
Reading the File M3e3_InstAccParts_snpsht11.dat
Reading the File M3e3_InstAccParts_snpsht12.dat
Reading the File M3e3_InstAccParts_snpsht13.dat
Reading the File M3e3_InstAccParts_snpsht14.dat
Reading the File M3e3_InstAccParts_snpsht15.dat
Reading the File M3e3_InstAccParts_snpsht16.dat
Reading the File M3e3_InstAccParts_snpsht17.dat
Reading the File M3e3_InstAccParts_snpsht18.dat
Reading the File M3e3_InstAccParts_snpsht19.dat
Reading the File M3e3_InstAccParts_snpsht20.dat
Reading the File M3e3_InstAccParts_snpsht21.dat
Reading the File M3e3_InstAccParts_snpsht22.dat
Reading the File M3e3_InstAccParts_snpsht23.dat
Reading the File M3e3_InstAccParts_snpsht24.dat
Reading the File M3e3_InstAccParts_snpsht25.dat
Reading the File M3e3_InstAccParts_snpsht26.dat
Reading the File M3e3_InstAccParts_snpsht27.dat
Reading the File M3e3_InstAccParts_snpsht28.dat
Reading the File M3e3_InstAccParts_snpsht29.dat
Reading the File M3e3_InstAccParts_snpsht30.dat
Reading the File M3e3_InstAccParts_snpsht31.dat
Reading the File M3e3_InstAccParts_snpsht32.dat
Reading the File M3e3_InstAccParts_snpsht33.dat
Reading the File M3e3_InstAccParts_snpsht34.dat
Reading the File M3e3_InstAccParts_snpsht35.dat
Reading the File M3e3_InstAccParts_snpsht36.dat
Reading the File M3e3_InstAccParts_snpsht37.dat
Reading the File M3e3_InstAccParts_snpsht38.dat
Reading the File M3e3_InstAccParts_snpsht39.dat
Reading the File M3e3_InstAccParts_snpsht40.dat
Reading the File M3e3_InstAccParts_snpsht41.dat
Reading the File M3e3_InstAccParts_snpsht42.dat
Reading the File M3e3_InstAccParts_snpsht43.dat
Reading the File M3e3_InstAccParts_snpsht44.dat
Reading the File M3e3_InstAccParts_snpsht45.dat
Reading the File M3e3_InstAccParts_snpsht46.dat
Reading the File M3e3_InstAccParts_snpsht47.dat
Reading the File M3e3_InstAccParts_snpsht48.dat
Reading the File M3e3_InstAccParts_snpsht49.dat
Reading the File M3e3_InstAccParts_snpsht50.dat
Reading the File M3e3_InstAccParts_snpsht51.dat
Reading the File M3e3_InstAccParts_snpsht52.dat
Reading the File M3e3_InstAccParts_snpsht53.dat
Reading the File M3e3_InstAccParts_snpsht54.dat
Reading the File M3e3_InstAccParts_snpsht55.dat
Reading the File M3e3_InstAccParts_snpsht56.dat
Reading the File M3e3_InstAccParts_snpsht57.dat
Reading the File M3e3_InstAccParts_snpsht58.dat
Reading the File M3e3_InstAccParts_snpsht59.dat
Reading the File M3e3_InstAccParts_snpsht60.dat
Reading the File M3e3_InstAccParts_snpsht61.dat
Reading the File M3e3_InstAccParts_snpsht62.dat
Reading the File M3e3_InstAccParts_snpsht63.dat
Reading the File M3e3_InstAccParts_snpsht64.dat
Reading the File M3e3_InstAccParts_snpsht65.dat
Reading the File M3e3_InstAccParts_snpsht66.dat
Reading the File M3e3_InstAccParts_snpsht67.dat
Reading the File M3e3_InstAccParts_snpsht68.dat
Reading the File M3e3_InstAccParts_snpsht69.dat
Reading the File M3e3_InstAccParts_snpsht70.dat
Reading the File M3e3_InstAccParts_snpsht71.dat
Reading the File M3e3_InstAccParts_snpsht72.dat
Reading the File M3e3_InstAccParts_snpsht73.dat
Reading the File M3e3_InstAccParts_snpsht74.dat
Reading the File M3e3_InstAccParts_snpsht75.dat
Reading the File M3e3_InstAccParts_snpsht76.dat
Reading the File M3e3_InstAccParts_snpsht77.dat
Reading the File M3e3_InstAccParts_snpsht78.dat
Reading the File M3e3_InstAccParts_snpsht79.dat
Reading the File M3e3_InstAccParts_snpsht80.dat
Reading the File M3e3_InstAccParts_snpsht81.dat
Reading the File M3e3_InstAccParts_snpsht82.dat
Reading the File M3e3_InstAccParts_snpsht83.dat
Reading the File M3e3_InstAccParts_snpsht84.dat
Reading the File M3e3_InstAccParts_snpsht85.dat
Reading the File M3e3_InstAccParts_snpsht86.dat
Reading the File M3e3_InstAccParts_snpsht87.dat
Reading the File M3e3_InstAccParts_snpsht88.dat
Reading the File M3e3_InstAccParts_snpsht89.dat
Reading the File M3e3_InstAccParts_snpsht90.dat
Reading the File M3e3_InstAccParts_snpsht91.dat
Reading the File M3e3_InstAccParts_snpsht92.dat
Reading the File M3e3_InstAccParts_snpsht93.dat
Reading the File M3e3_InstAccParts_snpsht94.dat
Reading the File M3e3_InstAccParts_snpsht95.dat
Reading the File M3e3_InstAccParts_snpsht96.dat
Reading the File M3e3_InstAccParts_snpsht97.dat
Reading the File M3e3_InstAccParts_snpsht98.dat
In [6]:
cloud_data, cloud_props = Rd.Restore_Cloud_Evolution(Cloud_name, resolution, tinit, tfin)
Reading the File M3e3_CloudObject_snp000.dat
Reading the File M3e3_CloudObject_snp000.dat
Reading the File M3e3_CloudObject_snp001.dat
Reading the File M3e3_CloudObject_snp002.dat
Reading the File M3e3_CloudObject_snp003.dat
Reading the File M3e3_CloudObject_snp004.dat
Reading the File M3e3_CloudObject_snp005.dat
Reading the File M3e3_CloudObject_snp006.dat
Reading the File M3e3_CloudObject_snp007.dat
Reading the File M3e3_CloudObject_snp008.dat
Reading the File M3e3_CloudObject_snp009.dat
Reading the File M3e3_CloudObject_snp010.dat
Reading the File M3e3_CloudObject_snp011.dat
Reading the File M3e3_CloudObject_snp012.dat
Reading the File M3e3_CloudObject_snp013.dat
Reading the File M3e3_CloudObject_snp014.dat
Reading the File M3e3_CloudObject_snp015.dat
Reading the File M3e3_CloudObject_snp016.dat
Reading the File M3e3_CloudObject_snp017.dat
Reading the File M3e3_CloudObject_snp018.dat
Reading the File M3e3_CloudObject_snp019.dat
Reading the File M3e3_CloudObject_snp020.dat
Reading the File M3e3_CloudObject_snp021.dat
Reading the File M3e3_CloudObject_snp022.dat
Reading the File M3e3_CloudObject_snp023.dat
Reading the File M3e3_CloudObject_snp024.dat
Reading the File M3e3_CloudObject_snp025.dat
Reading the File M3e3_CloudObject_snp026.dat
Reading the File M3e3_CloudObject_snp027.dat
Reading the File M3e3_CloudObject_snp028.dat
Reading the File M3e3_CloudObject_snp029.dat
Reading the File M3e3_CloudObject_snp030.dat
Reading the File M3e3_CloudObject_snp031.dat
Reading the File M3e3_CloudObject_snp032.dat
Reading the File M3e3_CloudObject_snp033.dat
Reading the File M3e3_CloudObject_snp034.dat
Reading the File M3e3_CloudObject_snp035.dat
Reading the File M3e3_CloudObject_snp036.dat
Reading the File M3e3_CloudObject_snp037.dat
Reading the File M3e3_CloudObject_snp038.dat
Reading the File M3e3_CloudObject_snp039.dat
Reading the File M3e3_CloudObject_snp040.dat
Reading the File M3e3_CloudObject_snp041.dat
Reading the File M3e3_CloudObject_snp042.dat
Reading the File M3e3_CloudObject_snp043.dat
Reading the File M3e3_CloudObject_snp044.dat
Reading the File M3e3_CloudObject_snp045.dat
Reading the File M3e3_CloudObject_snp046.dat
Reading the File M3e3_CloudObject_snp047.dat
Reading the File M3e3_CloudObject_snp048.dat
Reading the File M3e3_CloudObject_snp049.dat
Reading the File M3e3_CloudObject_snp050.dat
Reading the File M3e3_CloudObject_snp051.dat
Reading the File M3e3_CloudObject_snp052.dat
Reading the File M3e3_CloudObject_snp053.dat
Reading the File M3e3_CloudObject_snp054.dat
Reading the File M3e3_CloudObject_snp055.dat
Reading the File M3e3_CloudObject_snp056.dat
Reading the File M3e3_CloudObject_snp057.dat
Reading the File M3e3_CloudObject_snp058.dat
Reading the File M3e3_CloudObject_snp059.dat
Reading the File M3e3_CloudObject_snp060.dat
Reading the File M3e3_CloudObject_snp061.dat
Reading the File M3e3_CloudObject_snp062.dat
Reading the File M3e3_CloudObject_snp063.dat
Reading the File M3e3_CloudObject_snp064.dat
Reading the File M3e3_CloudObject_snp065.dat
Reading the File M3e3_CloudObject_snp066.dat
Reading the File M3e3_CloudObject_snp067.dat
Reading the File M3e3_CloudObject_snp068.dat
Reading the File M3e3_CloudObject_snp069.dat
Reading the File M3e3_CloudObject_snp070.dat
Reading the File M3e3_CloudObject_snp071.dat
Reading the File M3e3_CloudObject_snp072.dat
Reading the File M3e3_CloudObject_snp073.dat
Reading the File M3e3_CloudObject_snp074.dat
Reading the File M3e3_CloudObject_snp075.dat
Reading the File M3e3_CloudObject_snp076.dat
Reading the File M3e3_CloudObject_snp077.dat
Reading the File M3e3_CloudObject_snp078.dat
Reading the File M3e3_CloudObject_snp079.dat
Reading the File M3e3_CloudObject_snp080.dat
Reading the File M3e3_CloudObject_snp081.dat
Reading the File M3e3_CloudObject_snp082.dat
Reading the File M3e3_CloudObject_snp083.dat
Reading the File M3e3_CloudObject_snp084.dat
Reading the File M3e3_CloudObject_snp085.dat
Reading the File M3e3_CloudObject_snp086.dat
Reading the File M3e3_CloudObject_snp087.dat
Reading the File M3e3_CloudObject_snp088.dat
Reading the File M3e3_CloudObject_snp089.dat
Reading the File M3e3_CloudObject_snp090.dat
Reading the File M3e3_CloudObject_snp091.dat
Reading the File M3e3_CloudObject_snp092.dat
Reading the File M3e3_CloudObject_snp093.dat
Reading the File M3e3_CloudObject_snp094.dat
Reading the File M3e3_CloudObject_snp095.dat
Reading the File M3e3_CloudObject_snp096.dat
Reading the File M3e3_CloudObject_snp097.dat
Reading the File M3e3_CloudObject_snp098.dat
In [7]:
Nearby_SN = NSN.Locate_Nearby_SN(Cloud_name, min_dist=100)
NSN.Print_Nearby_SN(Nearby_SN)
('Total Number of SN:', 4627)
SN at d = 60.33 	, time = 1.17
dx = 30.00, dy = 51.06, dz = 11.48
SN at d = 79.85 	, time = 1.81
dx = 14.84, dy = 77.65, dz = -11.30
SN at d = 56.73 	, time = 2.62
dx = 7.24, dy = 51.48, dz = -22.70
SN at d = 72.57 	, time = 4.26
dx = 60.39, dy = 13.08, dz = 38.07
t = 1.17 	, d = 60.33 	, px = -942.23, py = -51.06, pz = 11.48
t = 1.81 	, d = 79.85 	, px = 14.84, py = -77.64, pz = -11.30
t = 2.62 	, d = 56.73 	, px = 7.24, py = 51.48, pz = -22.70
t = 4.26 	, d = 72.57 	, px = -911.85, py = -13.08, pz = 38.07

In [8]:
# Plot the number of accreted particles as a function of time
In [9]:
nparts_acc = np.zeros_like(particle_evol["time"])

for i in range(len(nparts_acc)):
    nparts_acc[i] = len(particle_evol["numdens"][i])

nparts_acc = np.array(nparts_acc)
In [10]:
fig = plt.figure(figsize=(12,10))

ax = fig.add_subplot(111)

ax.plot(particle_evol["time"], nparts_acc, "-k", linewidth=2)

ax.set_xlabel("time [Myr]", fontsize=15)
ax.set_ylabel("Number of acc particles", fontsize=15)
ax.set_title("Instantaneous accretion of particles. Cloud %s" %Cloud_name, fontsize=15)

for tick in ax.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax.yaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
    
fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEvolution_props/"
filename = Cloud_name + "_%s" %resolution_str + "_AccParticles"

fig.savefig(save_dir+filename+".pdf", format='pdf')

print("Saving Figure: %s%s.pdf" %(save_dir, filename))
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEvolution_props/M3e3_006_AccParticles.pdf
/home/jcibanezm/codes/ytAcc2/yt-x86_64/lib/python2.7/site-packages/matplotlib-1.4.3-py2.7-linux-x86_64.egg/matplotlib/figure.py:387: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [11]:
dist = np.histogram(np.log10(particle_evol["numdens"][0]), bins=100, normed=False)
    
nparts = dist[0] + 0.1
dens   = dist[1][0:-1]
In [12]:
#for i in range(len(particle_evol["numdens"][0])):
    #if particle_evol["numdens"][0][i] > 100:
        #print "This particle, tag = no tag, has a density of %f" %(particle_evol["numdens"][0][i])    
In [13]:
index_cut = 0
for i in range(len(dens)):
    if 10**dens[i] > 120.:
        print 10**dens[i], i
        index_cut = i
        break
        
if index_cut == 0:
    index_cut = len(nparts)
120.614628796 75

In [14]:
# Find the mean density of the incoming gas.
In [16]:
########################################################################
fig = plt.figure(figsize=(12,5))

nmin = -2 
nmax = 2.1

plot_range = [0.07, 0.1, 0.98, 0.87]

ax0  = fig.add_axes(plot_range, frameon=True)
#ax0  = fig.add_subplot(111)

#cmap = "gist_stern_r"
cmap = "gnuplot2"

#ax0.fill_between([0, particle_evol["time"][-1]], -3, 2.2, color='k', alpha=0.5)
ax0.set_axis_bgcolor('k')
for i in range(len(particle_evol["time"])-1):
    
    # This is the time bin where the data will be plotted.
    dt = particle_evol["time"][i+1] - particle_evol["time"][i]
    t0 = particle_evol["time"][i]   - dt/2.0
    t1 = particle_evol["time"][i+1] - dt/2.0 + 0.05
    
    dist = np.histogram(np.log10(particle_evol["numdens"][i]), range=(nmin, nmax), bins=50,  normed=False)
   
    nparts = dist[0] + 0.1
    dens   = dist[1][0:-1]

    index_cut = 0
    for i in range(len(dens)):
        if 10**dens[i] > 100.:
            #print 10**dens[i], i
            index_cut = i
            break

    if index_cut == 0:
        index_cut = len(nparts)

    nparts = dist[0][0:index_cut] + 0.1
    dens   = dist[1][0:index_cut]

    # Make a matrix with two entries of the histogram of velocities.
    matrix = []

    for i in range(len(dens)):
        matrix.append(0)
        matrix[i] = []
        for j in range(2):
            matrix[i].append(nparts[i])

    matrix = np.array(matrix)

    #extent = [t0, t1, np.max(dens), np.min(dens)]

    extent = [t0, t1, nmin, 2]
    cax = ax0.imshow(np.log10(matrix), cmap=cmap, extent=extent, vmin=-1, vmax=3.477, origin="lower")

    
#ax0.invert_yaxis()

ax0.set_aspect('auto')


ax0.spines["top"].set_visible(False)  
ax0.spines["right"].set_visible(True)  

ax0.get_xaxis().tick_bottom()  
ax0.get_yaxis().tick_left()  

ax0.set_ylim(nmin, 2.2)
ax0.set_xlim(0, t1)

ax0.set_ylabel("number density [cm$^{-3}$]", fontsize=15, labelpad=0)
ax0.set_xlabel("time [Myr]", fontsize=15)

yt     = [-2, -1, 0, 1, 2]
labels = ["10$^{-2}$", "10$^{-1}$", "10$^{0}$", "10$^{1}$", "10$^{2}$"]

# You can specify a rotation for the tick labels in degrees or with keywords.
plt.yticks(yt, labels)


for tick in ax0.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax0.yaxis.get_major_ticks():
    tick.label.set_fontsize(16) 

    
# Add the SN lines
yy = ax0.get_ylim()
for sn in range(len(Nearby_SN["time"])):
    # Give the SN line a width with respect to the distance to the cloud.
    lwidth = -(1./25.)*Nearby_SN["distance"][sn] + 4
    lalpha = -(1/100.)*Nearby_SN["distance"][sn] + 2.0
    if lalpha > 1: lalpha = 1
    print lwidth, lalpha
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth)
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
    
    if (Nearby_SN["distance"][sn] < 95):
        ax0.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--w', linewidth=lwidth, alpha=lalpha)
        ytxt = (yy[1] - yy[0])*0.95 + yy[0]
        #ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)
        

cbar = fig.colorbar(cax, aspect=12, shrink=1.0, pad=0.06)
cbar.set_label("Number of particles", fontsize=15)
cbar.set_ticks([-1, 0, 0.477, 1, 1.477, 2, 2.477, 3, 3.477])
cbar.set_ticklabels([0, 1, 3, 10, 30, 100, 300, 1000, 3000])
cbar.ax.tick_params(labelsize=14)

ax1 = ax0.twinx()

ax1.spines["top"].set_visible(False)  
ax1.spines["left"].set_visible(False)  
ax1.spines["right"].set_visible(True)  

ax1.plot(particle_evol["time"], np.log10(nparts_acc), "-w", linewidth=1.5)
#ax1.plot(particle_evol["time"], nparts_acc, "-w", linewidth=1.5)

#ax1.set_xlim(0, t1)
ax1.set_xlim(0, 10.1)


for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
    

#ax1.yaxis.set_label_position("right")
#ax1.set_ylabel("log$_{10}$(Accreted particles)", fontsize=15)
ax1.set_ylabel("log$_{10}$ (Num. particles) ", fontsize=14)

plt.tick_params(axis="y", labelcolor="k", labelsize=14)

start, end = ax0.get_xlim()
ax0.xaxis.set_ticks(np.arange(start, end, 1))

#fig.tight_layout()

Save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/"

fig.savefig("%s%s_InstAcc_ndens.pdf"%(Save_dir, Cloud_name), dpi=100)

print("Saving Figure: %s%s_InstAcc_ndens.pdf" %(Save_dir, Cloud_name))
1.58696553655 1
0.805846316861 1
1.73087893414 1
1.09705093261 1
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/M3e3_InstAcc_ndens.pdf
$$ v_{inf} = \left( \frac{G M}{2 R} \right)^{1/2} $$$$ v_{inf} = \left( \frac{2G M}{5 R} \right)^{1/2} $$
In [17]:
v_inf2 = np.sqrt(GNewton / 2.) * np.sqrt(np.array(cloud_props["mass"]) / np.array(cloud_props["radius"]))
vin    = np.sqrt(2.*GNewton / 5.) * np.sqrt(np.array(cloud_props["mass"]) / np.array(cloud_props["radius"]))

Compute the relative velocity with respect to the density gradient

In [18]:
# Make a function that takes two vectors vec_1 and vec_2 and returns the projected
# magnitude of vector 1, parallel and perpendicular to vector 2. 
# For example, vec1=velocity, vec2=density gradient.
# Returns, velocity parallel to density gradient and velocity perpendicular to density gradient.
In [19]:
V_par  = []
V_perp = []

for tt in range(len(particle_evol["time"])):
    
    #print ("Running time %.2f" %(particle_evol["time"][tt]))
    
    vrelx = np.array(particle_evol["velx"][tt]) - cloud_props["BV_x"][tt]
    vrely = np.array(particle_evol["vely"][tt]) - cloud_props["BV_y"][tt]
    vrelz = np.array(particle_evol["velz"][tt]) - cloud_props["BV_z"][tt]
    
    vec1 = np.array([vrelx, vrely, vrelz])
    vec2 = np.array([particle_evol["dens_grad_x"][tt], particle_evol["dens_grad_y"][tt], particle_evol["dens_grad_z"][tt]])

    vparallel      = vec.dot_product(vec1, vec2)
    vperpendicular = vec.cross_product(vec1, vec2)[1]
    
    V_par.append(vparallel)
    V_perp.append(vperpendicular)
    
V_par  = np.array(V_par)
V_perp = np.array(V_perp)
In [20]:
mean_vperp = []
mean_vpar  = []

# Find the Mean, parallel and perpendicular velocity.
for i in range(len(cloud_props["time"])):
    mean_vperp.append(np.mean(np.abs(V_perp[i])/1.0e5))
    mean_vpar.append(np.mean(np.abs(V_par[i])/1.0e5))

mean_vperp = np.array(mean_vperp)
mean_vpar  = np.array(mean_vpar)
In [21]:
########################################################################
fig = plt.figure(figsize=(12,4))


ax0  = fig.add_axes([0.08, 0.13, 0.98, 0.845], frameon=True)
ax0.set_axis_bgcolor('k')

cmap = "gnuplot2"

for i in range(len(particle_evol["time"])-1):

    # This is the time bin where the data will be plotted.
    dt = particle_evol["time"][i+1] - particle_evol["time"][i]
    t0 = particle_evol["time"][i]   - dt/2.0
    t1 = particle_evol["time"][i+1] - dt/2.5 
    
    dist = np.histogram(np.log10(np.abs(V_par[i])/1.0e5), range=(-4,2.2), bins=100, normed=False)
   
    nparts = dist[0] + 0.1
    dens   = dist[1][0:-1]

    index_cut = 0
    for i in range(len(dens)):
        if 10**dens[i] > 100.:
            #print 10**dens[i], i
            index_cut = i
            break

    if index_cut == 0:
        index_cut = len(nparts)

    nparts = dist[0][0:index_cut] + 0.1
    dens   = dist[1][0:index_cut]

    # Make a matrix with two entries of the histogram of velocities.
    matrix = []

    for i in range(len(dens)):
        matrix.append(0)
        matrix[i] = []
        for j in range(2):
            matrix[i].append(nparts[i])

    matrix = np.array(matrix)

    extent = [t0, t1, np.max(dens), np.min(dens)]

    cax = ax0.imshow(np.log10(matrix), cmap=cmap, extent=extent, vmin=0, vmax=3)
  

#ax0.plot(cloud_props["time"], np.log10(v_inf2/1.0e5), "-k", linewidth=1)
ax0.plot(cloud_props["time"], np.log10(vin/1.0e5), "-k", linewidth=2)

ax0.plot(cloud_props["time"], np.log10(mean_vpar), "--w", linewidth=2)
ax0.plot(cloud_props["time"], np.log10(mean_vperp), "ow-", linewidth=0.5, markersize=4)


yt     = [-2, -1, 0, 1, 2]
labels = ["10$^{-2}$", "10$^{-1}$", "10$^{0}$", "10$^{1}$", "10$^{2}$"]

# You can specify a rotation for the tick labels in degrees or with keywords.
plt.yticks(yt, labels)

ax0.set_aspect('auto')

ax0.spines["top"].set_visible(False)  
ax0.spines["right"].set_visible(False)  

ax0.get_xaxis().tick_bottom()  
ax0.get_yaxis().tick_left()  

#ax0.set_ylim(-3, 2.1)
ax0.set_ylim(-2, 2.2)
ax0.set_xlim(0, 10.1)

start, end = ax0.get_xlim()
ax0.xaxis.set_ticks(np.arange(start, end, 1))

#ax0.set_ylabel("log$_{10}$ (v$_{rel}$) [km s$^{-1}$]", fontsize=15)
ax0.set_ylabel("$v_{rel \, \parallel}$ [km s$^{-1}$]", fontsize=15)
ax0.set_xlabel("time [Myr]", fontsize=15)

for tick in ax0.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax0.yaxis.get_major_ticks():
    tick.label.set_fontsize(16) 
    
# Add the SN lines
yy = ax0.get_ylim()
for sn in range(len(Nearby_SN["time"])):
    # Give the SN line a width with respect to the distance to the cloud.
    lwidth = -(1./25.)*Nearby_SN["distance"][sn] + 4
    lalpha = -(1/100.)*Nearby_SN["distance"][sn] + 2.0
    if lalpha > 1: lalpha = 1
    print lwidth, lalpha
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth)
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
    
    if (Nearby_SN["distance"][sn] < 95):
        ax0.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--w', linewidth=lwidth, alpha=lalpha)
        ytxt = (yy[1] - yy[0])*0.95 + yy[0]
        #ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)
        

cbar = fig.colorbar(cax, aspect=12, shrink=1, pad=0.06)
cbar.set_label("Number of particles", fontsize=15)
cbar.set_ticks([-1, 0, 0.477, 1, 1.477, 2, 2.477, 3, 3.477])
cbar.set_ticklabels([0, 1, 3, 10, 30, 100, 300, 1000, 3000])
cbar.ax.tick_params(labelsize=14)

Save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/"

print("Saving Figure: %s%s_InstAcc_vrel_Par.pdf"%(Save_dir, Cloud_name))

fig.savefig(Save_dir + Cloud_name+"_InstAcc_vrel_Par.pdf", dpi=100)
1.58696553655 1
0.805846316861 1
1.73087893414 1
1.09705093261 1
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/M3e3_InstAcc_vrel_Par.pdf
-c:17: RuntimeWarning: divide by zero encountered in log10
In [22]:
########################################################################
fig = plt.figure(figsize=(12,4))

ax0  = fig.add_axes([0.08, 0.13, 0.98, 0.845], frameon=True)
ax0.set_axis_bgcolor('k')

cmap = "gnuplot2"

for i in range(len(particle_evol["time"])-1):
    
    # This is the time bin where the data will be plotted.
    dt = particle_evol["time"][i+1] - particle_evol["time"][i]
    t0 = particle_evol["time"][i]   - dt/2.0
    t1 = particle_evol["time"][i+1] - dt/2.5 
    
    dist = np.histogram(np.log10(V_perp[i]/1.0e5), range=(-3,2.5),  bins=100, normed=False)
   
    nparts = dist[0] + 0.1
    dens   = dist[1][0:-1]

    index_cut = 0
    for i in range(len(dens)):
        if 10**dens[i] > 100.:
            #print 10**dens[i], i
            index_cut = i
            break

    if index_cut == 0:
        index_cut = len(nparts)

    nparts = dist[0][0:index_cut] + 0.1
    dens   = dist[1][0:index_cut]

    # Make a matrix with two entries of the histogram of velocities.
    matrix = []

    for i in range(len(dens)):
        matrix.append(0)
        matrix[i] = []
        for j in range(2):
            matrix[i].append(nparts[i])

    matrix = np.array(matrix)

    extent = [t0, t1, np.max(dens), np.min(dens)]

    cax = ax0.imshow(np.log10(matrix), cmap=cmap, extent=extent, vmin=0, vmax=3)
  
#ax0.invert_yaxis()

#ax0.plot(global_cloud_props["time"], np.log10(v_inf/1.0e5), "--k", linewidth=1)

ax0.plot(cloud_props["time"], np.log10(vin/1.0e5), "-k", linewidth=2)

#ax0.plot(cloud_props["time"], np.log10(mean_vpar), "--", color="#92c5de", linewidth=2)
ax0.plot(cloud_props["time"], np.log10(mean_vpar), "--w", linewidth=2)
ax0.plot(cloud_props["time"], np.log10(mean_vperp), "ow-", linewidth=0.5, markersize=4)

yt     = [-2, -1, 0, 1, 2]
labels = ["10$^{-2}$", "10$^{-1}$", "10$^{0}$", "10$^{1}$", "10$^{2}$"]
# You can specify a rotation for the tick labels in degrees or with keywords.
plt.yticks(yt, labels)


ax0.set_aspect('auto')

ax0.spines["top"].set_visible(False)  
ax0.spines["right"].set_visible(False)  

ax0.get_xaxis().tick_bottom()  
ax0.get_yaxis().tick_left()  

#ax0.set_ylim(-3, 2.5)
ax0.set_ylim(-2, 2.2)
#ax0.set_xlim(0, t1)
ax0.set_xlim(0, 10.1)

start, end = ax0.get_xlim()
ax0.xaxis.set_ticks(np.arange(start, end, 1))

#ax0.set_ylabel("log$_{10}$ (v$_{rel}$) [km s$^{-1}$]", fontsize=15)
ax0.set_ylabel("$v_{rel \, \perp}$  [km s$^{-1}$]", fontsize=15, labelpad=0)
ax0.set_xlabel("time [Myr]", fontsize=15)

for tick in ax0.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax0.yaxis.get_major_ticks():
    tick.label.set_fontsize(16) 

# Add the SN lines
yy = ax0.get_ylim()
for sn in range(len(Nearby_SN["time"])):
    # Give the SN line a width with respect to the distance to the cloud.
    lwidth = -(1./25.)*Nearby_SN["distance"][sn] + 4
    lalpha = -(1/100.)*Nearby_SN["distance"][sn] + 2.0
    if lalpha > 1: lalpha = 1
    print lwidth, lalpha
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth)
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
    
    if (Nearby_SN["distance"][sn] < 95):
        ax0.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--w', linewidth=lwidth, alpha=lalpha)
        ytxt = (yy[1] - yy[0])*0.95 + yy[0]
        #ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)
        
    
cbar = fig.colorbar(cax, aspect=12, shrink=1.0, pad=0.06)
cbar.set_label("Number of particles", fontsize=18)
cbar.set_ticks([-1, 0, 0.477, 1, 1.477, 2, 2.477, 3, 3.477])
cbar.set_ticklabels([0, 1, 3, 10, 30, 100, 300, 1000, 3000])
cbar.ax.tick_params(labelsize=14)

Save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/"

fig.savefig(Save_dir + Cloud_name+"_InstAcc_vrel_Perp.pdf", dpi=100)

print("Saving Figure: %s%s_InstAcc_vrel_Perp.pdf" %(Save_dir, Cloud_name))
1.58696553655 1
0.805846316861 1
1.73087893414 1
1.09705093261 1
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/M3e3_InstAcc_vrel_Perp.pdf
-c:16: RuntimeWarning: divide by zero encountered in log10
In [22]:
 

In [23]:
# Compute the Mass Accretion rate

Mass_accreted = np.zeros_like(cloud_props["mass"])

for i in range(len(Mass_accreted)-1):
    Mass_accreted[i] = (cloud_props["mass"][i+1] - cloud_props["mass"][i]) / Msun 
                       

# Make the last input equal to the second to last input.
# Just to make the array equal size to the other arrays.
# Gravitational collapse is already unresolved by now.
Mass_accreted[-1] = Mass_accreted[-2]
In [24]:
# Now divide the accreted mass equally over all particles.
#E_kin_parts = np.zeros_like()
Ekin_in = np.zeros_like(cloud_props["time"])

for i in range(len(cloud_props["time"])):
    M_per_part = Mass_accreted[i] / nparts_acc[i]
    
    v2 = ((np.array(particle_evol["velx"][i]) - cloud_props["BV_x"][i])**2 + 
          (np.array(particle_evol["vely"][i]) - cloud_props["BV_y"][i])**2 + 
          (np.array(particle_evol["velz"][i]) - cloud_props["BV_z"][i])**2 ) / (1.0e5)**2

    e_kin_per_part = 0.5 * M_per_part * v2
    
    Ekin_in[i] = np.sum(e_kin_per_part)
    
Ekin_in = np.array(Ekin_in)

[Mper_part] = Mass

[v2] = Velocity^2 = Length^2 / Time^2

Ekin_in = Mass Length^2 / Time^2 = Energy.

If I divide this by the interval, it's the instantaneous accretion.

In [25]:
# Now, save the cloud data to a binary file.

import cPickle as cP

save_Ein= "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/InstantaneousAccretion"

Ein_dict = {"info":"This is the influx of kinetic energy as a function of time for cloud %s"%(Cloud_name)}

Ekin_in_cgs = Ekin_in * Msun * (1.0e5)**2

Ein_dict["Ein"]  = np.array(Ekin_in_cgs)
Ein_dict["time"] = np.array(particle_evol["time"])

f = open("%s/%s_Ein.dat" %(save_Ein, Cloud_name), 'wb')
cP.dump(Ein_dict, f, protocol=cP.HIGHEST_PROTOCOL)
f.close()

print("Saving the information of the influ of kinetic energy from the accretion")
Saving the information of the influ of kinetic energy from the accretion
In [26]:
# Calculate the expected Ein from my equation.
Const = 2*math.pi / 2.0**(2./3.) * (3/(4*math.pi))**(1./6.) * (6.6743e-8)**(3./2.)

avg_dens = np.array(cloud_props["mass"]) / np.array(cloud_props["volume"])

Ein_pred = Const * (10*mm) * np.array(cloud_props["mass"])**(5./3.) / avg_dens**(1./6.)
In [27]:
Ein_AstroUnit = Ein_pred / (Msun * 1.0e5**2)
Ein_01Myr = Ein_AstroUnit * 3.1556e12
In [28]:
fig = plt.figure(figsize=(12,10))

ax = fig.add_subplot(111)

ax.plot(particle_evol["time"], Ekin_in, "-k", linewidth=2)

ax.set_xlabel("time [Myr]", fontsize=15)
ax.set_ylabel("Accreted Kinetic Energy [M$_{\odot}$ km$^{2}$ s$^{-2}$]", fontsize=15)
ax.set_title("$\dot{E}_{in}$ . Cloud %s" %Cloud_name, fontsize=15)

for tick in ax.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax.yaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
    
max_here = np.max(Ekin_in)

ax.plot(cloud_props["time"], Ein_01Myr, "--r", linewidth=1)

ax.set_ylim(1.0e0, max_here*1.1)

ax.set_yscale("log")
    
ax.set_ylim(1.0e0, 1.0e6)
    
ax.set_xlim(0, 10)
    
fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEvolution_props/"
filename = Cloud_name + "_%s" %resolution_str + "_Ekin_in"

fig.savefig(save_dir+filename+".pdf", format='pdf')

print("Saving Figure: %s%s.pdf" %(save_dir, filename))
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEvolution_props/M3e3_006_Ekin_in.pdf

Compute the local Mach number for each accreted particle:

In [29]:
print particle_evol.keys()
['info', 'magz', 'magy', 'magx', 'accelz', 'accelx', 'temp', 'accely', 'velx', 'posy', 'posx', 'tag', 'posz', 'velz', 'dx', 'time', 'numdens', 'dens_grad_y', 'dens_grad_x', 'vely', 'dens_grad_z']
In [30]:
cs = np.sqrt(kb*particle_evol["temp"][0]*5.0/3.0 / mm)
In [31]:
Mach  = []

for tt in range(len(particle_evol["time"])):
    
    #print ("Running time %.2f" %(particle_evol["time"][tt]))
    
    vrelx = np.array(particle_evol["velx"][tt]) - cloud_props["BV_x"][tt]
    vrely = np.array(particle_evol["vely"][tt]) - cloud_props["BV_y"][tt]
    vrelz = np.array(particle_evol["velz"][tt]) - cloud_props["BV_z"][tt]
    
    vel_mag = np.sqrt(vrelx**2 + vrely**2 + vrelz**2)

    # Get the sounds speed

    cs = np.sqrt(kb*particle_evol["temp"][tt]*5.0/3.0 / mm)
    
    mach_here = vel_mag / cs
    
    Mach.append(mach_here)
    
Mach = np.array(Mach)
In [33]:
#Mach[90] = Mach[91]
In [34]:
########################################################################
fig = plt.figure(figsize=(12,4))


ax0  = fig.add_axes([0.08, 0.13, 0.98, 0.845], frameon=True)
ax0.set_axis_bgcolor('k')

cmap = "gnuplot2"

for i in range(len(particle_evol["time"])-1):

    # This is the time bin where the data will be plotted.
    dt = particle_evol["time"][i+1] - particle_evol["time"][i]
    t0 = particle_evol["time"][i]   - dt/2.0
    t1 = particle_evol["time"][i+1] - dt/2.5 
    
    dist = np.histogram(np.log10(np.abs(Mach[i])), range=(-1.5,1.5), bins=100, normed=False)
   
    nparts = dist[0] + 0.1
    dens   = dist[1][0:-1]

    index_cut = 0
    for i in range(len(dens)):
        if 10**dens[i] > 100.:
            #print 10**dens[i], i
            index_cut = i
            break

    if index_cut == 0:
        index_cut = len(nparts)

    nparts = dist[0][0:index_cut] + 0.1
    dens   = dist[1][0:index_cut]

    # Make a matrix with two entries of the histogram of velocities.
    matrix = []

    for i in range(len(dens)):
        matrix.append(0)
        matrix[i] = []
        for j in range(2):
            matrix[i].append(nparts[i])

    matrix = np.array(matrix)

    extent = [t0, t1, np.max(dens), np.min(dens)]

    cax = ax0.imshow(np.log10(matrix), cmap=cmap, extent=extent, vmin=0, vmax=3)
  

ax0.plot(cloud_props["time"], np.ones_like(cloud_props["time"])/10, "-w", linewidth=2)
#ax0.plot(cloud_props["time"], np.log10(vin/1.0e5), "-k", linewidth=2)

#ax0.plot(cloud_props["time"], np.log10(mean_vpar), "--w", linewidth=2)

yt     = [-2, -1, 0, 1, 2]
labels = ["10$^{-2}$", "10$^{-1}$", "10$^{0}$", "10$^{1}$", "10$^{2}$"]

# You can specify a rotation for the tick labels in degrees or with keywords.
plt.yticks(yt, labels)

ax0.set_aspect('auto')

ax0.spines["top"].set_visible(False)  
ax0.spines["right"].set_visible(False)  

ax0.get_xaxis().tick_bottom()  
ax0.get_yaxis().tick_left()  

#ax0.set_ylim(-3, 2.1)
ax0.set_ylim(-2, 2.2)
ax0.set_xlim(0, 10.1)

start, end = ax0.get_xlim()
ax0.xaxis.set_ticks(np.arange(start, end, 1))

#ax0.set_ylabel("log$_{10}$ (v$_{rel}$) [km s$^{-1}$]", fontsize=15)
ax0.set_ylabel("Mach Number", fontsize=15)
ax0.set_xlabel("time [Myr]", fontsize=15)

for tick in ax0.xaxis.get_major_ticks():
    tick.label.set_fontsize(14) 
for tick in ax0.yaxis.get_major_ticks():
    tick.label.set_fontsize(16) 
    
# Add the SN lines
yy = ax0.get_ylim()
for sn in range(len(Nearby_SN["time"])):
    # Give the SN line a width with respect to the distance to the cloud.
    lwidth = -(1./25.)*Nearby_SN["distance"][sn] + 4
    lalpha = -(1/100.)*Nearby_SN["distance"][sn] + 2.0
    if lalpha > 1: lalpha = 1
    print lwidth, lalpha
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth)
    #ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
    
    if (Nearby_SN["distance"][sn] < 95):
        ax0.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--w', linewidth=lwidth, alpha=lalpha)
        ytxt = (yy[1] - yy[0])*0.95 + yy[0]
        #ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)
        

cbar = fig.colorbar(cax, aspect=12, shrink=1, pad=0.06)
cbar.set_label("Number of particles", fontsize=15)
cbar.set_ticks([-1, 0, 0.477, 1, 1.477, 2, 2.477, 3, 3.477])
cbar.set_ticklabels([0, 1, 3, 10, 30, 100, 300, 1000, 3000])
cbar.ax.tick_params(labelsize=14)

Save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/"

print("Saving Figure: %s%s_InstAcc_Mach.pdf"%(Save_dir, Cloud_name))

fig.savefig(Save_dir + Cloud_name+"_InstAcc_Mach.pdf", dpi=100)
1.58696553655 1
0.805846316861 1
1.73087893414 1
1.09705093261 1
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/InstAccretion/M3e3_InstAcc_Mach.pdf
In [ ]: