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 Reconstruct_UG as RUG
import Vector_computations as vec

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

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

yt.add_field('numdens', function=numdens, units="1/cm**3", force_override=True)
In [2]:
Cloud_name = "M4e3"
resolution = 0.06

tinit  = 0.0
tfin   = 9.8

snapshot = 50

#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]:
particle_one = Rd.Read_One_GlobalAcc(Cloud_name, resolution, snapshot)
Reading the File M4e3_GlobalAccParts_snpsht50.dat
In [4]:
particle_evol = Rd.Read_GlobalAcc_Evolution(Cloud_name, resolution, tinit, tfin)
Reading the File M4e3_GlobalAccParts_snpsht00.dat
Reading the File M4e3_GlobalAccParts_snpsht00.dat
Reading the File M4e3_GlobalAccParts_snpsht01.dat
Reading the File M4e3_GlobalAccParts_snpsht02.dat
Reading the File M4e3_GlobalAccParts_snpsht03.dat
Reading the File M4e3_GlobalAccParts_snpsht04.dat
Reading the File M4e3_GlobalAccParts_snpsht05.dat
Reading the File M4e3_GlobalAccParts_snpsht06.dat
Reading the File M4e3_GlobalAccParts_snpsht07.dat
Reading the File M4e3_GlobalAccParts_snpsht08.dat
Reading the File M4e3_GlobalAccParts_snpsht09.dat
Reading the File M4e3_GlobalAccParts_snpsht10.dat
Reading the File M4e3_GlobalAccParts_snpsht11.dat
Reading the File M4e3_GlobalAccParts_snpsht12.dat
Reading the File M4e3_GlobalAccParts_snpsht13.dat
Reading the File M4e3_GlobalAccParts_snpsht14.dat
Reading the File M4e3_GlobalAccParts_snpsht15.dat
Reading the File M4e3_GlobalAccParts_snpsht16.dat
Reading the File M4e3_GlobalAccParts_snpsht17.dat
Reading the File M4e3_GlobalAccParts_snpsht18.dat
Reading the File M4e3_GlobalAccParts_snpsht19.dat
Reading the File M4e3_GlobalAccParts_snpsht20.dat
Reading the File M4e3_GlobalAccParts_snpsht21.dat
Reading the File M4e3_GlobalAccParts_snpsht22.dat
Reading the File M4e3_GlobalAccParts_snpsht23.dat
Reading the File M4e3_GlobalAccParts_snpsht24.dat
Reading the File M4e3_GlobalAccParts_snpsht25.dat
Reading the File M4e3_GlobalAccParts_snpsht26.dat
Reading the File M4e3_GlobalAccParts_snpsht27.dat
Reading the File M4e3_GlobalAccParts_snpsht28.dat
Reading the File M4e3_GlobalAccParts_snpsht29.dat
Reading the File M4e3_GlobalAccParts_snpsht30.dat
Reading the File M4e3_GlobalAccParts_snpsht31.dat
Reading the File M4e3_GlobalAccParts_snpsht32.dat
Reading the File M4e3_GlobalAccParts_snpsht33.dat
Reading the File M4e3_GlobalAccParts_snpsht34.dat
Reading the File M4e3_GlobalAccParts_snpsht35.dat
Reading the File M4e3_GlobalAccParts_snpsht36.dat
Reading the File M4e3_GlobalAccParts_snpsht37.dat
Reading the File M4e3_GlobalAccParts_snpsht38.dat
Reading the File M4e3_GlobalAccParts_snpsht39.dat
Reading the File M4e3_GlobalAccParts_snpsht40.dat
Reading the File M4e3_GlobalAccParts_snpsht41.dat
Reading the File M4e3_GlobalAccParts_snpsht42.dat
Reading the File M4e3_GlobalAccParts_snpsht43.dat
Reading the File M4e3_GlobalAccParts_snpsht44.dat
Reading the File M4e3_GlobalAccParts_snpsht45.dat
Reading the File M4e3_GlobalAccParts_snpsht46.dat
Reading the File M4e3_GlobalAccParts_snpsht47.dat
Reading the File M4e3_GlobalAccParts_snpsht48.dat
Reading the File M4e3_GlobalAccParts_snpsht49.dat
Reading the File M4e3_GlobalAccParts_snpsht50.dat
Reading the File M4e3_GlobalAccParts_snpsht51.dat
Reading the File M4e3_GlobalAccParts_snpsht52.dat
Reading the File M4e3_GlobalAccParts_snpsht53.dat
Reading the File M4e3_GlobalAccParts_snpsht54.dat
Reading the File M4e3_GlobalAccParts_snpsht55.dat
Reading the File M4e3_GlobalAccParts_snpsht56.dat
Reading the File M4e3_GlobalAccParts_snpsht57.dat
Reading the File M4e3_GlobalAccParts_snpsht58.dat
Reading the File M4e3_GlobalAccParts_snpsht59.dat
Reading the File M4e3_GlobalAccParts_snpsht60.dat
Reading the File M4e3_GlobalAccParts_snpsht61.dat
Reading the File M4e3_GlobalAccParts_snpsht62.dat
Reading the File M4e3_GlobalAccParts_snpsht63.dat
Reading the File M4e3_GlobalAccParts_snpsht64.dat
Reading the File M4e3_GlobalAccParts_snpsht65.dat
Reading the File M4e3_GlobalAccParts_snpsht66.dat
Reading the File M4e3_GlobalAccParts_snpsht67.dat
Reading the File M4e3_GlobalAccParts_snpsht68.dat
Reading the File M4e3_GlobalAccParts_snpsht69.dat
Reading the File M4e3_GlobalAccParts_snpsht70.dat
Reading the File M4e3_GlobalAccParts_snpsht71.dat
Reading the File M4e3_GlobalAccParts_snpsht72.dat
Reading the File M4e3_GlobalAccParts_snpsht73.dat
Reading the File M4e3_GlobalAccParts_snpsht74.dat
Reading the File M4e3_GlobalAccParts_snpsht75.dat
Reading the File M4e3_GlobalAccParts_snpsht76.dat
Reading the File M4e3_GlobalAccParts_snpsht77.dat
Reading the File M4e3_GlobalAccParts_snpsht78.dat
Reading the File M4e3_GlobalAccParts_snpsht79.dat
Reading the File M4e3_GlobalAccParts_snpsht80.dat
Reading the File M4e3_GlobalAccParts_snpsht81.dat
Reading the File M4e3_GlobalAccParts_snpsht82.dat
Reading the File M4e3_GlobalAccParts_snpsht83.dat
Reading the File M4e3_GlobalAccParts_snpsht84.dat
Reading the File M4e3_GlobalAccParts_snpsht85.dat
Reading the File M4e3_GlobalAccParts_snpsht86.dat
Reading the File M4e3_GlobalAccParts_snpsht87.dat
Reading the File M4e3_GlobalAccParts_snpsht88.dat
Reading the File M4e3_GlobalAccParts_snpsht89.dat
Reading the File M4e3_GlobalAccParts_snpsht90.dat
Reading the File M4e3_GlobalAccParts_snpsht91.dat
Reading the File M4e3_GlobalAccParts_snpsht92.dat
Reading the File M4e3_GlobalAccParts_snpsht93.dat
Reading the File M4e3_GlobalAccParts_snpsht94.dat
Reading the File M4e3_GlobalAccParts_snpsht95.dat
Reading the File M4e3_GlobalAccParts_snpsht96.dat
Reading the File M4e3_GlobalAccParts_snpsht97.dat
Reading the File M4e3_GlobalAccParts_snpsht98.dat
In [5]:
cloud_one = Rd.Read_One_Snapshot(Cloud_name, resolution, snapshot)
Reading the File M4e3_CloudObject_snp050.dat
In [6]:
cloud_data, cloud_props = Rd.Restore_Cloud_Evolution(Cloud_name, resolution, tinit, tfin)
Reading the File M4e3_CloudObject_snp000.dat
Reading the File M4e3_CloudObject_snp000.dat
Reading the File M4e3_CloudObject_snp001.dat
Reading the File M4e3_CloudObject_snp002.dat
Reading the File M4e3_CloudObject_snp003.dat
Reading the File M4e3_CloudObject_snp004.dat
Reading the File M4e3_CloudObject_snp005.dat
Reading the File M4e3_CloudObject_snp006.dat
Reading the File M4e3_CloudObject_snp007.dat
Reading the File M4e3_CloudObject_snp008.dat
Reading the File M4e3_CloudObject_snp009.dat
Reading the File M4e3_CloudObject_snp010.dat
Reading the File M4e3_CloudObject_snp011.dat
Reading the File M4e3_CloudObject_snp012.dat
Reading the File M4e3_CloudObject_snp013.dat
Reading the File M4e3_CloudObject_snp014.dat
Reading the File M4e3_CloudObject_snp015.dat
Reading the File M4e3_CloudObject_snp016.dat
Reading the File M4e3_CloudObject_snp017.dat
Reading the File M4e3_CloudObject_snp018.dat
Reading the File M4e3_CloudObject_snp019.dat
Reading the File M4e3_CloudObject_snp020.dat
Reading the File M4e3_CloudObject_snp021.dat
Reading the File M4e3_CloudObject_snp022.dat
Reading the File M4e3_CloudObject_snp023.dat
Reading the File M4e3_CloudObject_snp024.dat
Reading the File M4e3_CloudObject_snp025.dat
Reading the File M4e3_CloudObject_snp026.dat
Reading the File M4e3_CloudObject_snp027.dat
Reading the File M4e3_CloudObject_snp028.dat
Reading the File M4e3_CloudObject_snp029.dat
Reading the File M4e3_CloudObject_snp030.dat
Reading the File M4e3_CloudObject_snp031.dat
Reading the File M4e3_CloudObject_snp032.dat
Reading the File M4e3_CloudObject_snp033.dat
Reading the File M4e3_CloudObject_snp034.dat
Reading the File M4e3_CloudObject_snp035.dat
Reading the File M4e3_CloudObject_snp036.dat
Reading the File M4e3_CloudObject_snp037.dat
Reading the File M4e3_CloudObject_snp038.dat
Reading the File M4e3_CloudObject_snp039.dat
Reading the File M4e3_CloudObject_snp040.dat
Reading the File M4e3_CloudObject_snp041.dat
Reading the File M4e3_CloudObject_snp042.dat
Reading the File M4e3_CloudObject_snp043.dat
Reading the File M4e3_CloudObject_snp044.dat
Reading the File M4e3_CloudObject_snp045.dat
Reading the File M4e3_CloudObject_snp046.dat
Reading the File M4e3_CloudObject_snp047.dat
Reading the File M4e3_CloudObject_snp048.dat
Reading the File M4e3_CloudObject_snp049.dat
Reading the File M4e3_CloudObject_snp050.dat
Reading the File M4e3_CloudObject_snp051.dat
Reading the File M4e3_CloudObject_snp052.dat
Reading the File M4e3_CloudObject_snp053.dat
Reading the File M4e3_CloudObject_snp054.dat
Reading the File M4e3_CloudObject_snp055.dat
Reading the File M4e3_CloudObject_snp056.dat
Reading the File M4e3_CloudObject_snp057.dat
Reading the File M4e3_CloudObject_snp058.dat
Reading the File M4e3_CloudObject_snp059.dat
Reading the File M4e3_CloudObject_snp060.dat
Reading the File M4e3_CloudObject_snp061.dat
Reading the File M4e3_CloudObject_snp062.dat
Reading the File M4e3_CloudObject_snp063.dat
Reading the File M4e3_CloudObject_snp064.dat
Reading the File M4e3_CloudObject_snp065.dat
Reading the File M4e3_CloudObject_snp066.dat
Reading the File M4e3_CloudObject_snp067.dat
Reading the File M4e3_CloudObject_snp068.dat
Reading the File M4e3_CloudObject_snp069.dat
Reading the File M4e3_CloudObject_snp070.dat
Reading the File M4e3_CloudObject_snp071.dat
Reading the File M4e3_CloudObject_snp072.dat
Reading the File M4e3_CloudObject_snp073.dat
Reading the File M4e3_CloudObject_snp074.dat
Reading the File M4e3_CloudObject_snp075.dat
Reading the File M4e3_CloudObject_snp076.dat
Reading the File M4e3_CloudObject_snp077.dat
Reading the File M4e3_CloudObject_snp078.dat
Reading the File M4e3_CloudObject_snp079.dat
Reading the File M4e3_CloudObject_snp080.dat
Reading the File M4e3_CloudObject_snp081.dat
Reading the File M4e3_CloudObject_snp082.dat
Reading the File M4e3_CloudObject_snp083.dat
Reading the File M4e3_CloudObject_snp084.dat
Reading the File M4e3_CloudObject_snp085.dat
Reading the File M4e3_CloudObject_snp086.dat
Reading the File M4e3_CloudObject_snp087.dat
Reading the File M4e3_CloudObject_snp088.dat
Reading the File M4e3_CloudObject_snp089.dat
Reading the File M4e3_CloudObject_snp090.dat
Reading the File M4e3_CloudObject_snp091.dat
Reading the File M4e3_CloudObject_snp092.dat
Reading the File M4e3_CloudObject_snp093.dat
Reading the File M4e3_CloudObject_snp094.dat
Reading the File M4e3_CloudObject_snp095.dat
Reading the File M4e3_CloudObject_snp096.dat
Reading the File M4e3_CloudObject_snp097.dat
Reading the File M4e3_CloudObject_snp098.dat
In [7]:
# Generalize for all the other clouds.
if Cloud_name == "M4e3":

    px = 180
    py = -25
    pz =  10
    cloud_alpha = 0.4

    rad = 30

    px_str = str(px)
    py_str = str(py)
    pz_str = str(pz)

    pz = pz - 3.

    # Mean bulk velocity of the cloud.
    bvx_mean = -1.
    bvy_mean = -1.
    bvz_mean = -1.
    
    txt_xoffset = 15
    txt_yoffset = -25

elif Cloud_name == "M3e3":
    px = 450
    py = -380
    pz = 25
    cloud_alpha = 0.4

    rad = 45

    px_str = str(px)
    py_str = str(py)
    pz_str = str(pz)

    px = px + 8
    py = py 
    pz = pz - 8

    # Mean bulk velocity of the cloud.
    bvx_mean = 0.
    bvy_mean = 3.0
    bvz_mean = -3.0
    
    txt_xoffset = +30
    txt_yoffset = -30


elif Cloud_name == "M8e3":
    px = 60
    py = 370
    pz = 30
    cloud_alpha = 0.8

    rad = 60

    px_str = str(px)
    py_str = str(py)
    pz_str = str(pz)

    # Mean bulk velocity of the cloud.
    bvx_mean =  1.
    bvy_mean = -1.
    bvz_mean = -1.

    px = px + 5
    py = py - 10
    pz = pz - 5
    
    txt_xoffset = 30
    txt_yoffset = -60


cloud_dir = Cloud_name + "_" + "a%.2i" %(cloud_alpha*10) + "_x" + px_str + "_y" + py_str + "_z" + pz_str

if resolution < 0.1:
    resolution_str = "%.3i" %(resolution*100)
else:
    resolution_str = "%.2i" %(resolution*10)

data_dir = "/data/manda/jcibanezm/StratBox/RealProd/1pc_and_AccClouds/AccPaper_Resims/Particles/" + cloud_dir + "/%spc/" %resolution_str

# Set the address and the name of the data set.
basename = "Strat_Box_" + Cloud_name + "_%spc_hdf5_" %(resolution_str)
In [8]:
plt_file  = data_dir + basename + "plt_cnt_%.4i" %snapshot
part_file = data_dir + basename + "part_%.4i"    %snapshot

# Load data set
pf = yt.load(plt_file , particle_filename=part_file)

# Move the center of the frame.
px_now = px + bvx_mean * snapshot/10.
py_now = py + bvy_mean * snapshot/10.
pz_now = pz + bvz_mean * snapshot/10.

c = [px_now*ppc, py_now*ppc, pz_now*ppc]
#rad = 60

le = [c[0] - rad*ppc, c[1]-rad*ppc, c[2]-rad*ppc]
re = [c[0] + rad*ppc, c[1]+rad*ppc, c[2]+rad*ppc]

# Generate regions to access the particle and fluid data.
box  = pf.region(c, le, re)
sph  = pf.sphere(c, rad*pc)
yt : [INFO     ] 2017-04-21 09:17:24,898 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-04-21 09:17:24,899 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-04-21 09:17:24,900 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-04-21 09:17:24,900 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-04-21 09:17:24,901 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-04-21 09:17:24,933 Parameters: current_time              = 7.72494805221e+15
yt : [INFO     ] 2017-04-21 09:17:24,933 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-04-21 09:17:24,934 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-04-21 09:17:24,935 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-04-21 09:17:24,935 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-04-21 09:17:26,836 Loading field plugins.
yt : [INFO     ] 2017-04-21 09:17:26,838 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,839 Loaded astro (15 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,840 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,841 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,843 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,844 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,845 Loaded local (111 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,846 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,847 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-04-21 09:17:26,848 Loaded species (121 new fields)
In [9]:
UG = RUG.Reconstruct_UG(Cloud_name, resolution, snapshot)
Reading the File M4e3_CloudObject_snp050.dat
===============================================
Reconstructing a uniform grid for this cloud
Uniform Grid properties:
Lz = 28.90   Ly = 22.85   Lz = 28.19
dx = 0.06 pc
Nx = 489,    Ny = 387,    Nz = 476
xmin = 162.00      xmax = 190.90
ymin = -44.51      ymax = -21.66
zmin = -17.21      zmax = 10.98
===============================================
In [10]:
NxUG = np.shape(UG["x"])[0]
NyUG = np.shape(UG["x"])[1]
NzUG = np.shape(UG["x"])[2]

xminUG, xmaxUG = np.min(UG["x"]), np.max(UG["x"])
yminUG, ymaxUG = np.min(UG["y"]), np.max(UG["y"])
zminUG, zmaxUG = np.min(UG["z"]), np.max(UG["z"])
In [11]:
surface_x = np.zeros((NyUG, NzUG))
for i in range(NyUG):
    for j in range(NzUG):
        surface_x[i,j] = np.max(UG["numdens"][:, i, j])
        
surface_y = np.zeros((NzUG, NxUG))
for i in range(NzUG):
    for j in range(NxUG):
        surface_y[i,j] = np.max(UG["numdens"][j, :, i])
        
surface_z = np.zeros((NxUG, NyUG))
surf_temp = np.zeros((NxUG, NyUG))

for i in range(NxUG):
    for j in range(NyUG):
        surface_z[i,j] = np.max(UG["numdens"][i, j, :])
        surf_temp[i,j] = np.max(UG["temp"][i, j, :])
In [12]:
window_frac = 1.0

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

ax = fig.add_axes([0.035, 0.50, 0.43, 0.43])

los     = 'z'
extentz  = [le[1]/ppc + rad*(1.-window_frac), re[1]/ppc - rad*(1.-window_frac), le[0]/ppc + rad*(1.-window_frac), re[0]/ppc- rad*(1.-window_frac)]
prj     = pf.proj("numdens", los, center=c, data_source=box)
frb_z = prj.to_frb(window_frac*2*rad*pc, 256)["numdens"].value

cax = ax.imshow(np.rot90(np.log10(frb_z)), cmap="RdYlBu_r", extent=extentz, vmin=18, vmax=23)
#ax.invert_yaxis()

ax.contour(np.rot90(np.transpose(surface_z)),1, levels=[100], extent=[yminUG/ppc, ymaxUG/ppc, xmaxUG/ppc, xminUG/ppc], colors='k', linewidth=2)

ax.scatter(particle_one["posy"]/ppc, particle_one["posx"]/ppc, s=5, c='k')
ax.text(c[1]/ppc - 0.3*txt_yoffset*window_frac, c[0]/ppc - txt_xoffset*1.2*(window_frac), 't = %2.1f Myr' %(snapshot/10.), style='italic', bbox={'facecolor':'white', 'alpha':0.8, 'pad':10}, fontsize=15)

ax.set_xlim(extentz[0], extentz[1])
ax.set_ylim(extentz[3], extentz[2])

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) 

start, end = ax.get_xlim()
ax.xaxis.set_ticks(np.arange(start, end, 20))

#ax.xaxis.tick_top()

#plt.xticks()

ax.xaxis.set_tick_params(labeltop='on', labelbottom='off', labelsize=15)
    
ax.xaxis.set_label_position('top')

ax.set_xlabel("position y [pc]", fontsize=18)
ax.set_ylabel("position x [pc]", fontsize=18)

################
los     = 'x'
extentx  = [le[1]/ppc + rad*(1.-window_frac), re[1]/ppc - rad*(1.-window_frac), re[2]/ppc - rad*(1.-window_frac), le[2]/ppc + rad*(1.-window_frac)]
prj     = pf.proj("numdens", los, center=c, data_source=box)
frb_x = prj.to_frb(window_frac*2*rad*pc, 256)["numdens"].value

ax1 = fig.add_axes([0.035, 0.07, 0.43, 0.43])
cax = ax1.imshow(np.log10(frb_x), cmap="RdYlBu_r", extent=extentx, vmin=18, vmax=23)
ax1.invert_yaxis()
ax1.contour(np.rot90(surface_x),1, levels=[100], extent=[yminUG/ppc, ymaxUG/ppc, zmaxUG/ppc, zminUG/ppc], colors='k', linewidth=2)
ax1.scatter(particle_one["posy"]/ppc, particle_one["posz"]/ppc, s=5, c='k')

ax1.text(c[1]/ppc - 0.3*txt_yoffset*window_frac, c[2]/ppc - txt_xoffset*1.2*(window_frac), 't = %2.1f Myr' %(snapshot/10.), style='italic', bbox={'facecolor':'white', 'alpha':0.8, 'pad':10}, fontsize=15)


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

ax1.set_xlabel("position y [pc]", fontsize=18)
ax1.set_ylabel("position z [pc]", fontsize=18)

ax1.set_xlim(extentx[0], extentx[1])
ax1.set_ylim(extentx[3], extentx[2])

cbar_ax = fig.add_axes([0.42, 0.07, 0.025, 0.86])
cbar = fig.colorbar(cax, cax=cbar_ax)

cbar.set_label("Column Density [cm$^{-2}$]", fontsize=16)

cbar.set_ticks([18, 19, 20, 21, 22, 23,])
cbar.set_ticklabels(["10$^{18}$", "10$^{19}$", "10$^{20}$", "10$^{21}$", "10$^{22}$", "10$^{23}$"])

cbar.ax.tick_params(labelsize=13) 

for pp in range(10):
    parx = []
    pary = []
    parz = []

    for i in range(len(particle_evol["time"])):
        parx.append(particle_evol["posx"][i][pp])
        pary.append(particle_evol["posy"][i][pp])
        parz.append(particle_evol["posz"][i][pp])

    parx = np.array(parx)
    pary = np.array(pary)
    parz = np.array(parz)

    ax.plot(pary[0:snapshot]/ppc, parx[0:snapshot]/ppc,  ":k", linewidth=0.3)
    ax1.plot(pary[0:snapshot]/ppc, parz[0:snapshot]/ppc, ":k", linewidth=0.3)

start, end = ax1.get_xlim()
ax1.xaxis.set_ticks(np.arange(start, end, 20))

    
###########################################################################
# Plot some of the gas properties traced by the particles
###########################################################################

time = []
for i in range(len(particle_evol["time"])):
    time.append(particle_evol["time"][i])
time = np.array(time)

x00 = 0.585
y00 = 0.06

hpad = 0.02

height = 0.21
width  = 0.40

# --------------------------------------------------------------------------
# number density
# --------------------------------------------------------------------------

ax2 = fig.add_axes([x00, y00+3*height+3*hpad, width, height])

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

ax2.get_xaxis().tick_bottom()  
ax2.get_yaxis().tick_left()  

for pp in range(10):
    ndens = []

    for i in range(len(particle_evol["time"])):
        ndens.append(particle_evol["numdens"][i][pp])

    ndens = np.array(ndens)

    ax2.plot(time[snapshot:-1], np.log10(ndens[snapshot:-1]),  "-k", linewidth=1, alpha=0.5)
    ax2.plot(time[0:snapshot], np.log10(ndens[0:snapshot]),   "-k", linewidth=1)
    #ax2.plot(time, np.log10(ndens),  "-k", linewidth=1)
    
    
ax2.set_ylabel("log$_{10}$(n) [cm$^{-3}$]", fontsize=12)

ax2.set_ylim(-3.9, 6.2)

ax2.plot([0,10], [2,2], '--r', linewidth=0.8)

for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(0) 
    
# --------------------------------------------------------------------------
# velocity pperpendicular to the density gradient
# --------------------------------------------------------------------------
    
ax3 = fig.add_axes([x00, y00+2*height+2*hpad, 0.4, height])

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

ax3.get_xaxis().tick_bottom()  
ax3.get_yaxis().tick_left()  

for tick in ax3.xaxis.get_major_ticks():
    tick.label.set_fontsize(0) 


# --------------------------------------------------------------------------
# velocity parallel to the density gradient
# --------------------------------------------------------------------------

ax4 = fig.add_axes([x00, y00+height+hpad, width, height])

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

ax4.get_xaxis().tick_bottom()  
ax4.get_yaxis().tick_left()  

for tick in ax4.xaxis.get_major_ticks():
    tick.label.set_fontsize(0) 

# --------------------------------------------------------------------------
# add the velocities
# --------------------------------------------------------------------------
V_par  = []
V_perp = []

for tt in range(len(particle_evol["time"])):
    
    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)

for pp in range(10):
    vpart_par  = []
    vpart_perp = []
    
    for i in range(len(particle_evol["time"])):
        vpart_par.append(V_par[i][pp])
        vpart_perp.append(V_perp[i][pp])
    
    vpart_par  = np.array(np.abs(vpart_par))/1.0e5
    vpart_perp = np.array(np.abs(vpart_perp))/1.0e5

    #ax3.plot(time, np.log10(vpart_perp),  "-k", linewidth=1)
    #ax4.plot(time, np.log10(vpart_par),  "-k", linewidth=1)
    ax3.plot(time[0:snapshot], np.log10(vpart_perp[0:snapshot]),  "-k", linewidth=1)
    ax4.plot(time[0:snapshot], np.log10(vpart_par[0:snapshot]),  "-k", linewidth=1)
    ax3.plot(time[snapshot:-1], np.log10(vpart_perp[snapshot:-1]),  "-k", linewidth=1, alpha=0.5)
    ax4.plot(time[snapshot:-1], np.log10(vpart_par[snapshot:-1]),  "-k", linewidth=1, alpha=0.5)

    
ax3.set_ylabel("log$_{10}$(v$_{\perp}$) [km s$^{-1}$]", fontsize=12)
ax4.set_ylabel("log$_{10}$(v$_\parallel$) [km s$^{-1}$]", fontsize=12)

ax3.set_ylim(-2, 3)
ax4.set_ylim(-2, 3)


# --------------------------------------------------------------------------
# Mach numbers of the particles with respect to the gas.
# --------------------------------------------------------------------------

ax5 = fig.add_axes([x00, y00, width, height])

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

ax5.get_xaxis().tick_bottom()  
ax5.get_yaxis().tick_left()  

for tick in ax5.xaxis.get_major_ticks():
    tick.label.set_fontsize(12) 
    
# cs = sqrt( \gamma * kb * T / mm  )
for pp in range(10):
    Mach  = []
    
    for i in range(len(particle_evol["time"])):
        
        # calculate the adiabatic sound speed.
        cs = np.sqrt(5./3. * kb * particle_evol["temp"][i][pp] / mm)
        
        # calculate the relative velocity of the particle with respect to the cloud center of mass
        vrelx = np.array(particle_evol["velx"][i][pp]) - cloud_props["BV_x"][i]
        vrely = np.array(particle_evol["vely"][i][pp]) - cloud_props["BV_y"][i]
        vrelz = np.array(particle_evol["velz"][i][pp]) - cloud_props["BV_z"][i]
        
        # calculate the velocity magnitude
        vmag = np.sqrt(vrelx**2 + vrely**2 + vrelz**2)
        
        # calculate the mach number for this particle.
        Mach.append(vmag/cs)
                    
    #ax5.plot(time, Mach,  "-k", linewidth=1)
    ax5.plot(time[0:snapshot], Mach[0:snapshot],  "-k", linewidth=1)
    ax5.plot(time[snapshot:-1], Mach[snapshot:-1],  "-k", linewidth=1, alpha=0.5)

    
ax5.set_ylabel("Mach number", fontsize=12)
ax5.plot([0, 10],[1,1], '--r', linewidth=0.8)
ax5.set_yscale("log")
ax5.set_ylim(0.05, 20)

ax5.set_xlabel("time [Myr]", fontsize=14)

fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/General"
save_animation = "/data/gamera/jcibanezm/StratBox/AccretionPaper/KineticEnergy/Animations/GeneralAccretion/"

#fig.savefig("%s/%s_GeneralAcc_%spc.pdf" %(save_dir, Cloud_name, resolution_str), dpi=100)
#fig.savefig("%s/%s/%s_GeneralAcc_%spc_%.2i.png" %(save_animation, Cloud_name, Cloud_name, resolution_str, snapshot))

print("Saving figure: %s/%s_GeneralAcc_%spc.pdf" %(save_dir, Cloud_name, resolution_str))
yt : [INFO     ] 2017-04-21 09:17:59,885 Projection completed
yt : [INFO     ] 2017-04-21 09:17:59,889 Making a fixed resolution buffer of (numdens) 256 by 256
yt : [INFO     ] 2017-04-21 09:18:01,979 Projection completed
yt : [INFO     ] 2017-04-21 09:18:02,003 Making a fixed resolution buffer of (numdens) 256 by 256
Saving figure: /home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/General/M4e3_GeneralAcc_006pc.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 [13]:
fig = plt.figure(figsize=(14,6))

# --------------------------------------------------------------------------
# velocity parallel to the density gradient
# --------------------------------------------------------------------------

ax4 = fig.add_subplot(111)

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

ax4.get_xaxis().tick_bottom()  
ax4.get_yaxis().tick_left()  

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

# --------------------------------------------------------------------------
# add the velocities
# --------------------------------------------------------------------------
V_par  = []
V_perp = []

for tt in range(len(particle_evol["time"])):
#for tt in range(1):
    
    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)

for pp in range(10):
    
    vpart_par  = []
    vpart_perp = []
    
    for i in range(len(particle_evol["time"])):
        vpart_par.append(V_par[i][pp])
        vpart_perp.append(V_perp[i][pp])
    
    vpart_par  = np.array(np.abs(vpart_par))/1.0e5
    vpart_perp = np.array(vpart_perp)/1.0e5

    #ax3.plot(time, np.log10(vpart_perp),  "-k", linewidth=1)
    #ax4.plot(time, np.log10(vpart_par),  "-k", linewidth=1)
    #ax3.plot(time[0:snapshot], np.log10(vpart_perp[0:snapshot]),  "-k", linewidth=1)
    ax4.plot(time[0:snapshot], np.log10(vpart_par[0:snapshot]),  "-k", linewidth=1)
    #ax3.plot(time[snapshot:-1], np.log10(vpart_perp[snapshot:-1]),  "-k", linewidth=1, alpha=0.5)
    ax4.plot(time[snapshot-1:-1], np.log10(vpart_par[snapshot-1:-1]),  "-k", linewidth=1, alpha=0.5)

    
ax4.set_ylabel("log$_{10}$(v$_\parallel$) [km s$^{-1}$]", fontsize=12)

ax4.set_ylim(-2, 3)

fig.show()
In [14]:
fig = plt.figure(figsize=(14,12))

# --------------------------------------------------------------------------
# velocity parallel to the density gradient
# --------------------------------------------------------------------------

ax4 = fig.add_subplot(211)

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

ax4.get_xaxis().tick_bottom()  
ax4.get_yaxis().tick_left()  

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

# --------------------------------------------------------------------------
# add the velocities
# --------------------------------------------------------------------------
V_par  = []
V_perp = []

for tt in range(len(particle_evol["time"])):
#for tt in range(1):
    
    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["accel_x"][tt], particle_evol["accel_y"][tt], particle_evol["accel_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)

for pp in range(10):
    
    vpart_par  = []
    vpart_perp = []
    
    for i in range(len(particle_evol["time"])):
        vpart_par.append(V_par[i][pp])
        vpart_perp.append(V_perp[i][pp])
    
    vpart_par  = np.array(np.abs(vpart_par))/1.0e5
    vpart_perp = np.array(vpart_perp)/1.0e5

    #ax3.plot(time, np.log10(vpart_perp),  "-k", linewidth=1)
    #ax4.plot(time, np.log10(vpart_par),  "-k", linewidth=1)
    #ax3.plot(time[0:snapshot], np.log10(vpart_perp[0:snapshot]),  "-k", linewidth=1)
    ax4.plot(time[0:snapshot], np.log10(vpart_par[0:snapshot]),  "-k", linewidth=1)
    #ax3.plot(time[snapshot:-1], np.log10(vpart_perp[snapshot:-1]),  "-k", linewidth=1, alpha=0.5)
    ax4.plot(time[snapshot-1:-1], np.log10(vpart_par[snapshot-1:-1]),  "-k", linewidth=1, alpha=0.5)

    
ax4.set_ylabel("log$_{10}$(v$_\parallel$) [km s$^{-1}$]", fontsize=12)

ax4.set_ylim(-2, 3)


ax4 = fig.add_subplot(212)

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

ax4.get_xaxis().tick_bottom()  
ax4.get_yaxis().tick_left()  

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

for pp in range(10):
    
    vpart_par  = []
    vpart_perp = []
    
    for i in range(len(particle_evol["time"])):
        vpart_par.append(V_par[i][pp])
        vpart_perp.append(V_perp[i][pp])
    
    vpart_par  = np.array(np.abs(vpart_par))/1.0e5
    vpart_perp = np.array(vpart_perp)/1.0e5

    ax4.plot(time[0:snapshot], np.log10(vpart_perp[0:snapshot]),  "-k", linewidth=1)
    ax4.plot(time[snapshot-1:-1], np.log10(vpart_par[snapshot-1:-1]),  "-k", linewidth=1, alpha=0.5)
    
ax4.set_ylabel("log$_{10}$(v$_\perp$) [km s$^{-1}$]", fontsize=12)
ax4.set_ylim(-2, 3)

    
    
fig.show()