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 cPickle as cP

import Read_binary as Rd
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 * gram
mu      = 1.2924
kb      = 1.3806e-16   # erg K-1
GNewton = 6.6743e-8   * cm**3 * gram**(-1) * second**(-2)
Msun    = 1.9884e33   * gram
mm      = mu*mp
In [2]:
Cloud_name = "M3"


if Cloud_name == "M4":
    resolution = 0.06
    tunres       = 6.0
elif Cloud_name == "M8":
    resolution = 0.1
    tunres       = 6.3
elif Cloud_name == "M3":
    resolution = 0.06
    tunres       = 7.2
    
tinit  = 0.0
tfin   = 9.0


#Cloud_name = "M8e3"
#resolution = 0.5

#tinit  = 0.0
#tfin   = 7.2

if resolution < 0.1:
    resolution_str = "%.3i" %(resolution*100)
else:
    resolution_str = "%.2i" %(resolution*10)
In [3]:
c0, cp0 = Rd.Restore_One_Snapshot(Cloud_name+"e3", resolution, tinit, tfin)
Reading the File M3e3_CloudObject_snp000.dat
In [4]:
print "Cloud %s, global free fall time = %.2f Myr" %(Cloud_name+"e3", cp0["tff"]/3.1556e13)
Cloud M3e3, global free fall time = 3.48 Myr
In [5]:
def Read_Cloud_Gpot_Evolution(Cloud_name, resolution, tinit, tfin):

    global_data = {'info':'This dictionary contains the information of the cloud data as and its time evolution.'}

    global_data["time"] = []

    snapshot_data0 = Rd.Read_Gpot_cloud(Cloud_name, resolution, 0.)
    for key_name in snapshot_data0.keys():
        global_data[key_name] = []

    # Number of snapshots.
    nsnaps = int((tfin - tinit)*10)

    for i in range(nsnaps+1):
        tnow = tinit + i/10.
        snapshot_data = Rd.Read_Gpot_cloud(Cloud_name, resolution, i)
        global_data["time"].append(tnow)
        for key_name in snapshot_data0.keys():
            global_data[key_name].append(snapshot_data['%s'%key_name])

    return global_data
In [6]:
cloud_data, cloud_props = Rd.Restore_Cloud_Evolution(Cloud_name+"e3", resolution, tinit, tfin)
cloud_gpot              = Read_Cloud_Gpot_Evolution(Cloud_name+"e3", resolution, tinit, tfin)
#cloud_gpot = Rd.Read_Cloud_Gpot_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_Gpot_snpsht00.dat
Reading the File M3e3_Gpot_snpsht00.dat
Reading the File M3e3_Gpot_snpsht01.dat
Reading the File M3e3_Gpot_snpsht02.dat
Reading the File M3e3_Gpot_snpsht03.dat
Reading the File M3e3_Gpot_snpsht04.dat
Reading the File M3e3_Gpot_snpsht05.dat
Reading the File M3e3_Gpot_snpsht06.dat
Reading the File M3e3_Gpot_snpsht07.dat
Reading the File M3e3_Gpot_snpsht08.dat
Reading the File M3e3_Gpot_snpsht09.dat
Reading the File M3e3_Gpot_snpsht10.dat
Reading the File M3e3_Gpot_snpsht11.dat
Reading the File M3e3_Gpot_snpsht12.dat
Reading the File M3e3_Gpot_snpsht13.dat
Reading the File M3e3_Gpot_snpsht14.dat
Reading the File M3e3_Gpot_snpsht15.dat
Reading the File M3e3_Gpot_snpsht16.dat
Reading the File M3e3_Gpot_snpsht17.dat
Reading the File M3e3_Gpot_snpsht18.dat
Reading the File M3e3_Gpot_snpsht19.dat
Reading the File M3e3_Gpot_snpsht20.dat
Reading the File M3e3_Gpot_snpsht21.dat
Reading the File M3e3_Gpot_snpsht22.dat
Reading the File M3e3_Gpot_snpsht23.dat
Reading the File M3e3_Gpot_snpsht24.dat
Reading the File M3e3_Gpot_snpsht25.dat
Reading the File M3e3_Gpot_snpsht26.dat
Reading the File M3e3_Gpot_snpsht27.dat
Reading the File M3e3_Gpot_snpsht28.dat
Reading the File M3e3_Gpot_snpsht29.dat
Reading the File M3e3_Gpot_snpsht30.dat
Reading the File M3e3_Gpot_snpsht31.dat
Reading the File M3e3_Gpot_snpsht32.dat
Reading the File M3e3_Gpot_snpsht33.dat
Reading the File M3e3_Gpot_snpsht34.dat
Reading the File M3e3_Gpot_snpsht35.dat
Reading the File M3e3_Gpot_snpsht36.dat
Reading the File M3e3_Gpot_snpsht37.dat
Reading the File M3e3_Gpot_snpsht38.dat
Reading the File M3e3_Gpot_snpsht39.dat
Reading the File M3e3_Gpot_snpsht40.dat
Reading the File M3e3_Gpot_snpsht41.dat
Reading the File M3e3_Gpot_snpsht42.dat
Reading the File M3e3_Gpot_snpsht43.dat
Reading the File M3e3_Gpot_snpsht44.dat
Reading the File M3e3_Gpot_snpsht45.dat
Reading the File M3e3_Gpot_snpsht46.dat
Reading the File M3e3_Gpot_snpsht47.dat
Reading the File M3e3_Gpot_snpsht48.dat
Reading the File M3e3_Gpot_snpsht49.dat
Reading the File M3e3_Gpot_snpsht50.dat
Reading the File M3e3_Gpot_snpsht51.dat
Reading the File M3e3_Gpot_snpsht52.dat
Reading the File M3e3_Gpot_snpsht53.dat
Reading the File M3e3_Gpot_snpsht54.dat
Reading the File M3e3_Gpot_snpsht55.dat
Reading the File M3e3_Gpot_snpsht56.dat
Reading the File M3e3_Gpot_snpsht57.dat
Reading the File M3e3_Gpot_snpsht58.dat
Reading the File M3e3_Gpot_snpsht59.dat
Reading the File M3e3_Gpot_snpsht60.dat
Reading the File M3e3_Gpot_snpsht61.dat
Reading the File M3e3_Gpot_snpsht62.dat
Reading the File M3e3_Gpot_snpsht63.dat
Reading the File M3e3_Gpot_snpsht64.dat
Reading the File M3e3_Gpot_snpsht65.dat
Reading the File M3e3_Gpot_snpsht66.dat
Reading the File M3e3_Gpot_snpsht67.dat
Reading the File M3e3_Gpot_snpsht68.dat
Reading the File M3e3_Gpot_snpsht69.dat
Reading the File M3e3_Gpot_snpsht70.dat
Reading the File M3e3_Gpot_snpsht71.dat
Reading the File M3e3_Gpot_snpsht72.dat
Reading the File M3e3_Gpot_snpsht73.dat
Reading the File M3e3_Gpot_snpsht74.dat
Reading the File M3e3_Gpot_snpsht75.dat
Reading the File M3e3_Gpot_snpsht76.dat
Reading the File M3e3_Gpot_snpsht77.dat
Reading the File M3e3_Gpot_snpsht78.dat
Reading the File M3e3_Gpot_snpsht79.dat
Reading the File M3e3_Gpot_snpsht80.dat
Reading the File M3e3_Gpot_snpsht81.dat
Reading the File M3e3_Gpot_snpsht82.dat
Reading the File M3e3_Gpot_snpsht83.dat
Reading the File M3e3_Gpot_snpsht84.dat
Reading the File M3e3_Gpot_snpsht85.dat
Reading the File M3e3_Gpot_snpsht86.dat
Reading the File M3e3_Gpot_snpsht87.dat
Reading the File M3e3_Gpot_snpsht88.dat
Reading the File M3e3_Gpot_snpsht89.dat
Reading the File M3e3_Gpot_snpsht90.dat
In [7]:
tt = 0

e_grav = []
e_kin  = []
e_mag  = []

for tt in range(len(cloud_props["time"])):
    # Compute the normalization of the gravitational potential.
    #Egrav       = np.abs(np.sum(cloud_data["cell_mass"][tt]*cloud_data["gpot"][tt]))
    #expec_Egrav =  3./5.*GNewton.value*cloud_props["mass"][tt]**2 / cloud_props["radius"][tt]
    #phi_norm    =  (Egrav - expec_Egrav) / cloud_props["mass"][tt]

    # Calculate the gravitaitonal potential energy density
    #epot = np.sum(cloud_data["numdens"][tt]*mm*(np.abs(cloud_data["gpot"][tt]) - phi_norm))

    # Compute the normalization of the gravitational potential.
    Egrav       = np.abs(np.sum(cloud_data["cell_mass"][tt]*cloud_gpot["phi"][tt]))

    # Calculate the gravitaitonal potential energy density
    epot = np.sum(cloud_data["numdens"][tt]*mm*(np.abs(cloud_gpot["phi"][tt])))

    
    # Calculate the relative velocity to the center of mass
    vrelx = np.array(cloud_data["velx"][tt] - cloud_props["BV_x"][tt])
    vrely = np.array(cloud_data["vely"][tt] - cloud_props["BV_y"][tt])
    vrelz = np.array(cloud_data["velz"][tt] - cloud_props["BV_z"][tt])

    # velocity magnitude
    vmag = np.sqrt(vrelx**2 + vrely**2 + vrelz**2)

    # Kinetic energy density
    ekin = np.sum(0.5*cloud_data["numdens"][tt]*mm*vmag*vmag)

    # Calculate the magnitude of the magnetic field strength
    Bmag = np.sqrt(cloud_data["magx"][tt]*cloud_data["magx"][tt] + 
                   cloud_data["magy"][tt]*cloud_data["magy"][tt] + 
                   cloud_data["magz"][tt]*cloud_data["magz"][tt]) 

    # magnetic energy density
    emag = np.sum(Bmag*Bmag/(8*math.pi))
    
    e_grav.append(epot)
    e_kin.append(ekin)
    e_mag.append(emag)
    
e_grav = np.array(e_grav)
e_kin = np.array(e_kin)
e_mag = np.array(e_mag)
In [8]:
print cloud_data.keys()
['info', 'magz', 'magy', 'magx', 'dx', 'temp', 'cell_mass', 'dy', 'gpot', 'dz', 'velz', 'vely', 'time', 'y', 'x', 'cell_volume', 'numdens', 'velx', 'z']
In [9]:
tt = 0

E_grav = []
E_kin  = []
E_mag  = []
E_therm= []

for tt in range(len(cloud_props["time"])):    
    # Compute the normalization of the gravitational potential.
    #Egrav       = np.abs(np.sum(cloud_data["cell_mass"][tt]*cloud_data["gpot"][tt]))
    #expec_Egrav =  3./5.*GNewton.value*cloud_props["mass"][tt]**2 / cloud_props["radius"][tt]
    #phi_norm    =  (Egrav - expec_Egrav) / cloud_props["mass"][tt]
    
    # Calculate the gravitaitonal potential energy density
    #Epot = np.sum(cloud_data["cell_mass"][tt]*(np.abs(cloud_data["gpot"][tt]) - phi_norm))

    
    # Calculate the gravitaitonal potential energy density.
    Epot = np.sum(cloud_data["cell_mass"][tt]*(np.abs(cloud_gpot["phi"][tt])))

    # Calculate the relative velocity to the center of mass
    vrelx = np.array(cloud_data["velx"][tt] - cloud_props["BV_x"][tt])
    vrely = np.array(cloud_data["vely"][tt] - cloud_props["BV_y"][tt])
    vrelz = np.array(cloud_data["velz"][tt] - cloud_props["BV_z"][tt])

    # velocity magnitude
    vmag = np.sqrt(vrelx**2 + vrely**2 + vrelz**2)

    # Kinetic energy density
    Ekin = np.sum(0.5*cloud_data["cell_mass"][tt]*vmag*vmag)

    # Calculate the magnitude of the magnetic field strength
    Bmag = np.sqrt(cloud_data["magx"][tt]*cloud_data["magx"][tt] + 
                   cloud_data["magy"][tt]*cloud_data["magy"][tt] + 
                   cloud_data["magz"][tt]*cloud_data["magz"][tt]) 

    # magnetic energy density
    Emag = np.sum(Bmag*Bmag*cloud_data["cell_volume"][tt]/(8*math.pi))
    
    # Pressure is an energy density.
    # Pressure times volume is energy.
    # P = nkt
    # Etherm = nkt*Vol
    Etherm = np.sum(cloud_data["numdens"][tt]*kb*cloud_data["temp"][tt]*cloud_data["cell_volume"][tt])
    
    E_grav.append(Epot)
    E_kin.append(Ekin)
    E_mag.append(Emag)
    E_therm.append(Etherm)
    
E_grav = np.array(E_grav)
E_kin = np.array(E_kin)
E_mag = np.array(E_mag)
E_therm = np.array(E_therm)

print("----Done----")
----Done----
In [10]:
Nearby_SN = NSN.Locate_Nearby_SN(Cloud_name+"e3", 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

Plot the cloud energetics

In [15]:
# Energy Density

fig = plt.figure(figsize=(10,8))

ax = fig.add_subplot(111)

if Cloud_name == "M4":
    yy = [1.0e-10, 2.0e-4]
elif Cloud_name == "M8":
    yy = [1.0e-10, 2.0e-4]
elif Cloud_name == "M3":
    yy = [1.0e-8, 1.0e-3]

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["time"][sn] < tunres:
        if (Nearby_SN["distance"][sn] < 95):
            ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
            ytxt = (yy[1] - yy[0])*0.8 + yy[0]
            ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)


ax.plot(cloud_props["time"], e_grav, "-b", linewidth=2, label="Gravitational energy")
ax.plot(cloud_props["time"], e_kin,  "-r", linewidth=2, label="Kinetic energy")
ax.plot(cloud_props["time"], e_mag, "--g", linewidth=2, label="Magnetic energy")


ax.legend(loc=4, fontsize=18)

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) 
    

ax.set_xlabel("time [Myr]", fontsize=18)
ax.set_ylabel("energy density [erg cm$^{-3}$]", fontsize=18)

t0 = 0
t1 = int(cloud_props["time"][-1]+0.5)

#ax.set_xlim(t0, t1)
ax.set_xlim(t0, tunres)

#ax.set_xticks(np.arange(t1+1))
ax.set_xticks(np.arange(tunres+1))

ax.set_ylim(yy[0], yy[1])

ax.set_yscale("log")

ax.set_title("Cloud %s"%(Cloud_name), fontsize=18)

fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEnergetics"
#fig.savefig("%s/%s_Energetics_%spc.pdf" %(save_dir, Cloud_name, resolution_str))
1.58696553655 1
0.805846316861 1
1.73087893414 1
1.09705093261 1
In [16]:
save_Ein= "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/InstantaneousAccretion"

f = open("%s/%s_Ein.dat" %(save_Ein, Cloud_name+"e3"), 'r')
Ein = cP.load(f)
f.close()
In [17]:
print Ein.keys()
['info', 'Ein', 'time']
In [18]:
if Cloud_name == "M3":
    yy = [1.0e45, 5.0e48]
elif Cloud_name == "M4":
    yy = [1.0e45, 1.0e49]
else:
    yy = [1.0e46, 1.0e50]
In [25]:
fig = plt.figure(figsize=(10,8))

ax = fig.add_subplot(111)

#yy = [1.0e44, 2.0e49]

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
    
    if Nearby_SN["time"][sn] < tunres:
        if (Nearby_SN["distance"][sn] < 95):
            ax.plot([Nearby_SN["time"][sn], Nearby_SN["time"][sn]], yy, '--b', linewidth=lwidth, alpha=lalpha)
            ytxt = (yy[1] - yy[0])*0.75 + yy[0]
            ax.text(Nearby_SN["time"][sn]+0.05, ytxt, "%.2f pc" %(Nearby_SN["distance"][sn]), rotation=90, fontsize=12)


ax.plot(cloud_props["time"], E_grav, ":ob", linewidth=0, markersize=4, label="Gravitational", )
ax.plot(cloud_props["time"], E_kin,  "--r", linewidth=2, label="Kinetic")
ax.plot(cloud_props["time"], E_mag,  "-d", color="#31a354", linewidth=1, markersize=6, label="Magnetic")
ax.plot(cloud_props["time"], E_therm, "-y", linewidth=2.5, label="Thermal")

#cumulative = np.cumsum(Ein["Ein"])
# plot the cumulative function
#ax.plot(Ein["time"], Ein["Ein"], ":", color="k", linewidth=2, label="Inflow")
#ax.plot(Ein["time"], cumulative, ":", color="k", linewidth=2, label="Inflow")


#ax.legend(loc=1, fontsize=18)

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) 
    

ax.set_xlabel("time [Myr]", fontsize=18)
ax.set_ylabel("Energy [erg]", fontsize=18)

t0 = 0
t1 = int(tunres+0.5)
#t1 = int(cloud_props["time"][-1]+0.5)

ax.set_xlim(t0, tunres)
ax.set_xticks(np.arange(t1+1))


ax.set_yscale("log")

ax.set_ylim(yy[0], yy[1])

ax.set_title("Cloud %s"%(Cloud_name), fontsize=18)

ax.legend(loc=0, fontsize=16, scatterpoints=4, numpoints=4)

# We change the fontsize of minor ticks label 
ax.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)

fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEnergetics"
#fig.savefig("%s/%s_TotalEnergy_%spc.pdf" %(save_dir, Cloud_name+"e3", resolution_str), format="pdf", dpi=100)

print("Saving figure %s/%s_TotalEnergy_%spc.pdf" %(save_dir, Cloud_name+"e3", resolution_str))
Saving figure /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEnergetics/M3e3_TotalEnergy_006pc.pdf

Energy rate of change.

$$ \frac{d E}{d t}$$
In [136]:
def dE_dt(Energy, time):
    """
    Given the energy as a function of time, return the rate of change.
    Returns the answer in units of erg per s.
    """
    
    Edot     = np.zeros_like(Energy)
    time_now = np.zeros_like(time)
    
    dt = (time[1] - time[0]) * 3.1556e13
    
    for i in range(len(Energy)):
        if i == 0:
            Edot[i]     = (Energy[i+1]-Energy[i]) / dt
            time_now[i] = time[i]
        elif i != 0 and i!= len(Energy)-1:
            Edot[i]     = (Energy[i+1]-Energy[i-1]) / (2*dt)
            time_now[i] = time[i]
        elif i == len(Energy):
            Edot[i]     = (Energy[i]-Energy[i-1]) / dt
            time_now[i] = time[i]
    
    return np.array(Edot), np.array(time_now)
In [137]:
def dE_dt_5point(Energy, time):
    """
    Given the energy as a function of time, return the rate of change.
    Returns the answer in units of erg per s.
    """
    
    Edot     = np.zeros_like(Energy)
    time_now = np.zeros_like(time)
    
    dt = (time[1] - time[0]) * 3.1556e13
    
    for i in range(len(Energy)):
        if i == 0:
            Edot[i]     = (Energy[i+1]-Energy[i]) / dt
            time_now[i] = time[i]
        elif i == 1 or i == 2:
            Edot[i]     = (Energy[i+1]-Energy[i-1]) / (2*dt)
            time_now[i] = time[i]
        elif i >= 3 and i<= len(Energy)-3:
            Edot[i]     = (-Energy[i+2] + 8*Energy[i+2] -8*Energy[i-1] + Energy[i-2]) / (12*dt)
            time_now[i] = time[i]
        elif i == len(Energy) - 2:
            Edot[i]     = (Energy[i+1]-Energy[i-1]) / (2*dt)
            time_now[i] = time[i]
        elif i == len(Energy) - 1:
            Edot[i]     = (Energy[i]-Energy[i-1]) / dt
            time_now[i] = time[i]
    
    return np.array(Edot), np.array(time_now)
In [138]:
Edot_grav, new_time = dE_dt(E_grav, cloud_props["time"])

Edot_kin, new_time = dE_dt(E_kin, cloud_props["time"])

Edot_mag, new_time = dE_dt(E_mag, cloud_props["time"])
In [139]:
Edot_grav, new_time = dE_dt_5point(E_grav, cloud_props["time"])

Edot_kin, new_time = dE_dt_5point(E_kin, cloud_props["time"])

Edot_mag, new_time = dE_dt_5point(E_mag, cloud_props["time"])
In [140]:
Edot_in = np.array(Ein["Ein"]) / (0.1*3.1556e13)
In [141]:
# Predicted decay rate of the kinetic energy:

Edec_const = GNewton**(3./2.) / (2* 5**(3./2.)) * (10 * Msun)**(5./2.) / (1 * pc)**(5./2.)

alpha = 2.3
exponent = (5 * (alpha-1)) / (2*alpha)

Edot_decay = Edec_const* (np.array(cloud_props["mass"])/(10*1.9884e33))**(exponent)
Edec_new = 1./2. * np.array(cloud_props["mass"]) * np.array(cloud_props["vel_disp_total"])**3 / np.array(cloud_props["radius"])
In [142]:
Edot_kin_pos = np.zeros_like(Edot_kin)
Edot_kin_neg = np.zeros_like(Edot_kin)

for i in range(len(Edot_kin)):
    Ehere = Edot_kin[i]
    if Ehere > 0:
        Edot_kin_pos[i] = Ehere
        Edot_kin_neg[i] = 0
    else:
        Edot_kin_neg[i] = np.abs(Ehere)
        Edot_kin_pos[i] = 0
In [143]:
fig = plt.figure(figsize=(10,8))

ax = fig.add_subplot(111)


#ax.plot(cloud_props["time"], Edot_grav, "-b", linewidth=3, label="Grav. pot. energy rate of change")
#ax.plot(cloud_props["time"], Edot_kin, "-r", linewidth=3, label="Kinetic energy")

ax.plot(cloud_props["time"], Edot_kin_pos, "-r", linewidth=3, label="$\dot{E}_{kin}$ Positive")
ax.plot(cloud_props["time"], Edot_kin_neg, ":r", linewidth=3, label="$\dot{E}_{kin}$ Negative")


ax.plot(Ein["time"], Edot_in, ":k", linewidth=3, label="Accreted energy")
ax.plot(cloud_props["time"], Edec_new, "--k", linewidth=2, label="Energy decay rate")

#ax.plot(cloud_props["time"], Edec_new, "--o", linewidth=2, label="Energy Decay")


#ax.legend(loc=1, fontsize=18)

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) 
    

ax.set_xlabel("time [Myr]", fontsize=18)
ax.set_ylabel("$\dot{E}$ [erg s$^{-1}$]", fontsize=18)

t0 = 0
t1 = int(cloud_props["time"][-1]+0.5)

ax.set_xlim(t0, tunres)
#ax.set_xticks(np.arange(t1+1))

#ax.set_ylim()

ax.legend(loc=1, fontsize=18)

ax.set_yscale("log")
#ax.set_ylim(yy[0], yy[1])
ax.set_ylim(1.0e31, 1.0e37)

ax.set_title("Cloud %s"%(Cloud_name), fontsize=18)

# We change the fontsize of minor ticks label 
ax.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)

fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEnergetics"
fig.savefig("%s/%se3_Energy_rates_%spc.pdf" %(save_dir, Cloud_name, resolution_str), format="pdf", dpi=100)
print("Saving Figure: %s/%se3_Energy_rates_%spc.pdf" %(save_dir, Cloud_name, resolution_str))
Saving Figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/CloudEnergetics/M4e3_Energy_rates_006pc.pdf
In [143]: