Virial parameter sub-section of the accreting clouds paper

  • Load the cloud vicinity data
  • Load the Cloud gravitational potential data
  • compute the energy densities at each cell

  • pick up a slice (e.g. middle of the cloud) and plot this slice.

  • Plot three different evolutionary times and plot one next to each other. (e.g. 2,4,6, or 2, 5, 8).
In [35]:
import yt
import math
from yt import derived_field
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm
from yt.units import pc, kpc, second, Kelvin, gram, erg, cm, km

import Read_binary as Rd
import Reconstruct_UG as RUG
import Vector_computations as vec

import copy

import matplotlib
from matplotlib.ticker import MaxNLocator
from matplotlib import cm as cmap
from mpl_toolkits.mplot3d import Axes3D

from scipy.interpolate import griddata

#import plotly.plotly as py
#from plotly.graph_objs import *
#from plotly.tools import FigureFactory as FF

from skimage import measure

%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 [36]:
Cloud_name = "M3e3"

# Distance to the center of mass of the cloud
dcm = 0

npts3D = 100j
In [37]:
# Generalize for all the other clouds.
if Cloud_name == "M4e3":

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

    rad = 20

    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
    
    resolution = 0.06
    resolution3D = 0.5
    
    cloud_plt_name = "M4"

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

    rad = 25

    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

    resolution = 0.06
    resolution3D = 0.06
    
    npts3D = 60j
    
    cloud_plt_name = "M3"

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

    rad = 40

    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

    resolution = 0.1
    resolution3D = 0.5
    
    cloud_plt_name = "M8"

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 [38]:
# Restore the cloud properties and data
cloud_t2, cloud_props_t2    = Rd.Restore_One_Snapshot(Cloud_name, resolution, 20)
cloud_t4, cloud_props_t4    = Rd.Restore_One_Snapshot(Cloud_name, resolution, 40)
cloud_t6, cloud_props_t6    = Rd.Restore_One_Snapshot(Cloud_name, resolution, 60)

# Read the Gravitational potential of the self-gravitating cloud only.
cloud_gpot_t2 = Rd.Read_Gpot_cloud(Cloud_name, resolution, 20)
cloud_gpot_t4 = Rd.Read_Gpot_cloud(Cloud_name, resolution, 40)
cloud_gpot_t6 = Rd.Read_Gpot_cloud(Cloud_name, resolution, 60)
Reading the File M3e3_CloudObject_snp020.dat
Reading the File M3e3_CloudObject_snp040.dat
Reading the File M3e3_CloudObject_snp060.dat
Reading the File M3e3_Gpot_snpsht20.dat
Reading the File M3e3_Gpot_snpsht40.dat
Reading the File M3e3_Gpot_snpsht60.dat
In [39]:
# Generate a Uniform Grid for the Snapshot at time t=4 Myr. For the 3D plot.

UG = RUG.Reconstruct_UG(Cloud_name, resolution3D, 40)

# Compute some normalization of the cloud extent.
minxUG, minyUG, minzUG       = np.min(UG["x"])/ppc, np.min(UG["y"])/ppc, np.min(UG["z"])/ppc
maxxUG, maxyUG, maxzUG       = np.max(UG["x"])/ppc, np.max(UG["y"])/ppc, np.max(UG["z"])/ppc
xrangeUG, yrangeUG, zrangeUG = maxxUG - minxUG, maxyUG - minyUG, maxzUG - minzUG
xnump, ynump, znump          = np.shape(UG["x"])[0], np.shape(UG["x"])[1], np.shape(UG["x"])[2]
xscale, yscale, zscale       = xrangeUG/xnump, yrangeUG/ynump, zrangeUG/znump


#Nx = np.shape(UG["numdens"])[0]
#Ny = np.shape(UG["numdens"])[1]
#Nz = np.shape(UG["numdens"])[2]
    
#pts = np.array([[i,j, k] for i in range(Nx) for j in range(Ny) for k in range(Nz)] )
#grid_x, grid_y, grid_z = np.mgrid[0:Nx:npts3D, 0:Ny:npts3D, 0:Nz:npts3D]
#dens = np.reshape(UG["numdens"], (np.shape(pts)[0]))

#grid_new = griddata(pts, dens, (grid_x, grid_y, grid_z), method='linear')
#verts, faces = measure.marching_cubes(grid_new, 100)
    
# Generate the isosurface contour at n=100 using marching cubes.
verts, faces = measure.marching_cubes(UG["numdens"], 100)
    
xmin = 0
ymin = 0
zmin = 0

verts[:,0] = (verts[:,0] * xscale) + minxUG
verts[:,1] = (verts[:,1] * yscale) + minyUG
verts[:,2] = (verts[:,2] * zscale) + minzUG
Reading the File M3e3_CloudObject_snp040.dat
===============================================
Reconstructing a uniform grid for this cloud
Uniform Grid properties:
Lz = 37.98   Ly = 49.02   Lz = 27.00
dx = 0.06 pc
Nx = 641,    Ny = 828,    Nz = 457
xmin = 433.30      xmax = 471.28
ymin = -390.70      ymax = -341.68
zmin = -3.26      zmax = 23.74
===============================================
In [40]:
##############################################################################################################################
#                             t = 2 Myr
##############################################################################################################################

plt_file  = data_dir + basename + "plt_cnt_%.4i" %20
part_file = data_dir + basename + "part_%.4i"    %20

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

# Move the center of the frame.
px_now_t2 = px + bvx_mean * 2.
py_now_t2 = py + bvy_mean * 2.
pz_now_t2 = pz + bvz_mean * 2.

c_t2 = [px_now_t2*ppc, py_now_t2*ppc, pz_now_t2*ppc]
#rad = 60

le_t2 = [c_t2[0] - rad*ppc, c_t2[1]-rad*ppc, c_t2[2]-rad*ppc]
re_t2 = [c_t2[0] + rad*ppc, c_t2[1]+rad*ppc, c_t2[2]+rad*ppc]

# Generate regions to access the particle and fluid data.
box_t2  = pf_t2.region(c_t2, le_t2, re_t2)
sph_t2  = pf_t2.sphere(c_t2, rad*pc)


##############################################################################################################################
#                             t = 4 Myr
##############################################################################################################################

plt_file  = data_dir + basename + "plt_cnt_%.4i" %40
part_file = data_dir + basename + "part_%.4i"    %40

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

# Move the center of the frame.
px_now_t4 = px + bvx_mean * 4.
py_now_t4 = py + bvy_mean * 4.
pz_now_t4 = pz + bvz_mean * 4.

c_t4 = [px_now_t4*ppc, py_now_t4*ppc, pz_now_t4*ppc]
#rad = 60

le_t4 = [c_t4[0] - rad*ppc, c_t4[1]-rad*ppc, c_t4[2]-rad*ppc]
re_t4 = [c_t4[0] + rad*ppc, c_t4[1]+rad*ppc, c_t4[2]+rad*ppc]

# Generate regions to access the particle and fluid data.
box_t4  = pf_t4.region(c_t4, le_t4, re_t4)
sph_t4  = pf_t4.sphere(c_t4, rad*pc)

##############################################################################################################################
#                             t = 6 Myr
##############################################################################################################################

plt_file  = data_dir + basename + "plt_cnt_%.4i" %60
part_file = data_dir + basename + "part_%.4i"    %60

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

# Move the center of the frame.
px_now_t6 = px + bvx_mean * 6.
py_now_t6 = py + bvy_mean * 6.
pz_now_t6 = pz + bvz_mean * 6.

c_t6 = [px_now_t6*ppc, py_now_t6*ppc, pz_now_t6*ppc]
#rad = 60

le_t6 = [c_t6[0] - rad*ppc, c_t6[1]-rad*ppc, c_t6[2]-rad*ppc]
re_t6 = [c_t6[0] + rad*ppc, c_t6[1]+rad*ppc, c_t6[2]+rad*ppc]

# Generate regions to access the particle and fluid data.
box_t6  = pf_t6.region(c_t6, le_t6, re_t6)
sph_t6  = pf_t6.sphere(c_t6, rad*pc)
yt : [INFO     ] 2017-03-30 05:58:29,565 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:29,566 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:29,567 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:29,567 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:29,567 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:29,583 Parameters: current_time              = 7.61450049514e+15
yt : [INFO     ] 2017-03-30 05:58:29,583 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-03-30 05:58:29,584 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:29,585 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:29,586 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-03-30 05:58:30,387 Loading field plugins.
yt : [INFO     ] 2017-03-30 05:58:30,391 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,392 Loaded astro (15 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,393 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,395 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,397 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,399 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,400 Loaded local (111 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,401 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,402 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:30,403 Loaded species (121 new fields)
yt : [INFO     ] 2017-03-30 05:58:32,519 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:32,520 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:32,520 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:32,520 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:32,521 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:32,535 Parameters: current_time              = 7.6776210289e+15
yt : [INFO     ] 2017-03-30 05:58:32,535 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-03-30 05:58:32,536 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:32,537 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:32,538 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-03-30 05:58:33,486 Loading field plugins.
yt : [INFO     ] 2017-03-30 05:58:33,488 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,489 Loaded astro (15 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,490 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,491 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,493 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,494 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,495 Loaded local (111 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,496 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,497 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:33,498 Loaded species (121 new fields)
yt : [INFO     ] 2017-03-30 05:58:35,610 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:35,611 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:35,611 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:35,612 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:35,612 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-03-30 05:58:35,627 Parameters: current_time              = 7.74072560113e+15
yt : [INFO     ] 2017-03-30 05:58:35,627 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-03-30 05:58:35,628 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:35,629 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-03-30 05:58:35,629 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-03-30 05:58:37,756 Loading field plugins.
yt : [INFO     ] 2017-03-30 05:58:37,757 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,758 Loaded astro (15 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,759 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,761 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,762 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,764 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,765 Loaded local (111 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,766 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,766 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-03-30 05:58:37,767 Loaded species (121 new fields)
In [41]:
#phi_diff0_t2 = abs(cloud_t2["gpot"][0]) - abs(cloud_gpot_t2["phi"][0])
#phi_diff0_t4 = abs(cloud_t4["gpot"][0]) - abs(cloud_gpot_t4["phi"][0])
#phi_diff0_t6 = abs(cloud_t6["gpot"][0]) - abs(cloud_gpot_t6["phi"][0])

phi_diff0_t2 = abs(np.mean(cloud_t2["gpot"] - cloud_gpot_t2["phi"]))
phi_diff0_t4 = abs(np.mean(cloud_t4["gpot"] - cloud_gpot_t4["phi"]))
phi_diff0_t6 = abs(np.mean(cloud_t6["gpot"] - cloud_gpot_t6["phi"]))

In [42]:
# Center the plot on the center of mass at t=4 Myr in the xy plane, and on each clouds center of mass in z.
#CM_t2 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t2["CM_z"]]
#CM_t4 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t4["CM_z"]]
#CM_t6 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t6["CM_z"]]
In [43]:
# center of the Slice, accounting for the distance to the center of mass:
CM_t2 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t2["CM_z"] + dcm*ppc]
CM_t4 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t4["CM_z"] + dcm*ppc]
CM_t6 = [cloud_props_t4["CM_x"], cloud_props_t4["CM_y"], cloud_props_t6["CM_z"] + dcm*ppc] 
In [44]:
num_vel_vectors = 16

##############################################################################################################################
#                             t = 2 Myr
##############################################################################################################################

extent_t2  = [CM_t2[0]/ppc-rad, CM_t2[0]/ppc+rad, CM_t2[1]/ppc-rad, CM_t2[1]/ppc+rad]
slc_size_x = ( extent_t2[1] - extent_t2[0] ) * pc
slc_size_y = ( extent_t2[1] - extent_t2[0] ) * pc

slc_t2      = yt.SlicePlot(pf_t2, 'z', fields=['gpot'], center=CM_t2)
slc_frb_t2  = slc_t2.data_source.to_frb(slc_size_x, 1024, periodic=True)

slc_ndens_t2      = yt.SlicePlot(pf_t2, 'z', fields=['numdens'], center=CM_t2)
slc_frb_ndens_t2  = slc_ndens_t2.data_source.to_frb(slc_size_x, 1024, periodic=True)

velx2 = (np.array(slc_frb_t2["velx"].in_cgs()) - cloud_props_t2["BV_x"] )**2
vely2 = (np.array(slc_frb_t2["vely"].in_cgs()) - cloud_props_t2["BV_y"] )**2
velz2 = (np.array(slc_frb_t2["velz"].in_cgs()) - cloud_props_t2["BV_z"] )**2

vel_mag_t2 = np.sqrt(velx2 + vely2 + velz2)

kin_edens_t2 = (0.5 * np.array(slc_frb_t2["dens"].in_cgs().value * vel_mag_t2**2))

gpot_norm_t2  = np.array(slc_frb_t2['gpot'].in_cgs().value) + phi_diff0_t2
gpot_edens_t2 = gpot_norm_t2 * np.array(slc_frb_t2['dens'].in_cgs().value)
#gpot_edens_abs = abs(gpot_norm * np.array(slc_frb['dens'].in_cgs().value)) + 1.0e-16

thermal_energy_t2 = (slc_frb_t2["thermal_energy"]*slc_frb_t2["dens"]).in_units("erg/cm**3").value
mag_ener_dens_t2  =  slc_frb_t2["magnetic_energy"].in_units("erg/cm**3").value

total_energy_t2 = kin_edens_t2 + gpot_edens_t2 + thermal_energy_t2 + mag_ener_dens_t2

# --> --> --> --> --> --> -->  Velocity vectors <-- <-- <-- <-- <-- <-- <-- <-- <-- <--

# Generate the velocity arrays.
# --------------------------------------------------------------------------------------

slc_vx1 = yt.SlicePlot(pf_t2, 'z', 'velx', center=CM_t2, origin="native")
slc_vy1 = yt.SlicePlot(pf_t2, 'z', 'vely', center=CM_t2, origin="native")

slc_frb_vx1_t2 = slc_vx1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)
slc_frb_vy1_t2 = slc_vy1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)



##############################################################################################################################
#                             t = 4 Myr
##############################################################################################################################

extent_t4 = [CM_t4[0]/ppc-rad, CM_t4[0]/ppc+rad, CM_t4[1]/ppc-rad, CM_t4[1]/ppc+rad]
slc_size_x = ( extent_t4[1] - extent_t4[0] ) * pc
slc_size_y = ( extent_t4[1] - extent_t4[0] ) * pc

slc_t4      = yt.SlicePlot(pf_t4, 'z', fields=['gpot'], center=CM_t4)
slc_frb_t4  = slc_t4.data_source.to_frb(slc_size_x, 1024, periodic=True)

slc_ndens_t4      = yt.SlicePlot(pf_t4, 'z', fields=['numdens'], center=CM_t4)
slc_frb_ndens_t4  = slc_ndens_t4.data_source.to_frb(slc_size_x, 1024, periodic=True)

velx2 = (np.array(slc_frb_t4["velx"].in_cgs()) - cloud_props_t4["BV_x"] )**2
vely2 = (np.array(slc_frb_t4["vely"].in_cgs()) - cloud_props_t4["BV_y"] )**2
velz2 = (np.array(slc_frb_t4["velz"].in_cgs()) - cloud_props_t4["BV_z"] )**2

vel_mag_t4 = np.sqrt(velx2 + vely2 + velz2)

kin_edens_t4 = (0.5 * np.array(slc_frb_t4["dens"].in_cgs().value * vel_mag_t4**2))

gpot_norm_t4  = np.array(slc_frb_t4['gpot'].in_cgs().value) + phi_diff0_t4
gpot_edens_t4 = gpot_norm_t4 * np.array(slc_frb_t4['dens'].in_cgs().value)
#gpot_edens_abs = abs(gpot_norm * np.array(slc_frb['dens'].in_cgs().value)) + 1.0e-16

thermal_energy_t4 = (slc_frb_t4["thermal_energy"]*slc_frb_t4["dens"]).in_units("erg/cm**3").value
mag_ener_dens_t4  =  slc_frb_t4["magnetic_energy"].in_units("erg/cm**3").value

total_energy_t4 = kin_edens_t4 + gpot_edens_t4 + thermal_energy_t4 + mag_ener_dens_t4


# --> --> --> --> --> --> -->  Velocity vectors <-- <-- <-- <-- <-- <-- <-- <-- <-- <--

# Generate the velocity arrays.
# --------------------------------------------------------------------------------------

slc_vx1 = yt.SlicePlot(pf_t4, 'z', 'velx', center=CM_t4, origin="native")
slc_vy1 = yt.SlicePlot(pf_t4, 'z', 'vely', center=CM_t4, origin="native")

slc_frb_vx1_t4 = slc_vx1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)
slc_frb_vy1_t4 = slc_vy1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)


##############################################################################################################################
#                             t = 6 Myr
##############################################################################################################################

extent_t6 = [CM_t6[0]/ppc-rad, CM_t6[0]/ppc+rad, CM_t6[1]/ppc-rad, CM_t6[1]/ppc+rad]
slc_size_x = ( extent_t6[1] - extent_t6[0] ) * pc
slc_size_y = ( extent_t6[1] - extent_t6[0] ) * pc

slc_t6      = yt.SlicePlot(pf_t6, 'z', fields=['gpot'], center=CM_t6)
slc_frb_t6  = slc_t6.data_source.to_frb(slc_size_x, 1024, periodic=True)

slc_ndens_t6      = yt.SlicePlot(pf_t6, 'z', fields=['numdens'], center=CM_t6)
slc_frb_ndens_t6  = slc_ndens_t6.data_source.to_frb(slc_size_x, 1024, periodic=True)

velx2 = (np.array(slc_frb_t6["velx"].in_cgs()) - cloud_props_t6["BV_x"] )**2
vely2 = (np.array(slc_frb_t6["vely"].in_cgs()) - cloud_props_t6["BV_y"] )**2
velz2 = (np.array(slc_frb_t6["velz"].in_cgs()) - cloud_props_t6["BV_z"] )**2

vel_mag_t6 = np.sqrt(velx2 + vely2 + velz2)

kin_edens_t6 = (0.5 * np.array(slc_frb_t6["dens"].in_cgs().value * vel_mag_t6**2))

gpot_norm_t6  = np.array(slc_frb_t6['gpot'].in_cgs().value) + phi_diff0_t6
gpot_edens_t6 = gpot_norm_t6 * np.array(slc_frb_t6['dens'].in_cgs().value)
#gpot_edens_abs = abs(gpot_norm * np.array(slc_frb['dens'].in_cgs().value)) + 1.0e-16

thermal_energy_t6 = (slc_frb_t6["thermal_energy"]*slc_frb_t6["dens"]).in_units("erg/cm**3").value
mag_ener_dens_t6  =  slc_frb_t6["magnetic_energy"].in_units("erg/cm**3").value

total_energy_t6 = kin_edens_t6 + gpot_edens_t6 + thermal_energy_t6 + mag_ener_dens_t6

# --> --> --> --> --> --> -->  Velocity vectors <-- <-- <-- <-- <-- <-- <-- <-- <-- <--

# Generate the velocity arrays.
# --------------------------------------------------------------------------------------

slc_vx1 = yt.SlicePlot(pf_t6, 'z', 'velx', center=CM_t6, origin="native")
slc_vy1 = yt.SlicePlot(pf_t6, 'z', 'vely', center=CM_t6, origin="native")

slc_frb_vx1_t6 = slc_vx1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)
slc_frb_vy1_t6 = slc_vy1.data_source.to_frb(slc_size_x, num_vel_vectors, height=slc_size_y, periodic=True)
yt : [INFO     ] 2017-03-30 05:58:40,634 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:40,634 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:40,638 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:40,638 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:40,640 Making a fixed resolution buffer of (('flash', 'gpot')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:40,815 Plot image for field ('flash', 'gpot') has no positive values.  Max = -21525480603644.984375.
yt : [WARNING  ] 2017-03-30 05:58:40,815 Switching to linear colorbar scaling.
yt : [INFO     ] 2017-03-30 05:58:41,192 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:41,193 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:41,196 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:41,197 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:41,198 Making a fixed resolution buffer of (('gas', 'numdens')) 800 by 800
yt : [INFO     ] 2017-03-30 05:58:41,649 Making a fixed resolution buffer of (velx) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:41,810 Making a fixed resolution buffer of (vely) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:41,949 Making a fixed resolution buffer of (velz) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:42,124 Making a fixed resolution buffer of (dens) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:42,205 Making a fixed resolution buffer of (gpot) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:42,260 Making a fixed resolution buffer of (thermal_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:42,338 Making a fixed resolution buffer of (magnetic_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:42,732 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:42,732 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:42,735 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:42,736 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:42,737 Making a fixed resolution buffer of (('flash', 'velx')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:42,925 Plot image for field ('flash', 'velx') has both positive and negative values. Min = -26107056.000001, Max = 40873084.000001.
yt : [WARNING  ] 2017-03-30 05:58:42,926 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2017-03-30 05:58:43,191 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:43,191 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:43,195 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:43,195 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:43,197 Making a fixed resolution buffer of (('flash', 'vely')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:43,380 Plot image for field ('flash', 'vely') has both positive and negative values. Min = -33082184.000001, Max = 38171231.470971.
yt : [WARNING  ] 2017-03-30 05:58:43,381 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2017-03-30 05:58:43,977 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:43,978 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:43,981 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:43,981 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:43,983 Making a fixed resolution buffer of (('flash', 'gpot')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:44,190 Plot image for field ('flash', 'gpot') has no positive values.  Max = -21513778495484.988281.
yt : [WARNING  ] 2017-03-30 05:58:44,190 Switching to linear colorbar scaling.
yt : [INFO     ] 2017-03-30 05:58:44,430 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:44,431 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:44,434 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:44,434 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:44,436 Making a fixed resolution buffer of (('gas', 'numdens')) 800 by 800
yt : [INFO     ] 2017-03-30 05:58:44,944 Making a fixed resolution buffer of (velx) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,037 Making a fixed resolution buffer of (vely) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,124 Making a fixed resolution buffer of (velz) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,220 Making a fixed resolution buffer of (dens) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,317 Making a fixed resolution buffer of (gpot) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,370 Making a fixed resolution buffer of (thermal_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,456 Making a fixed resolution buffer of (magnetic_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:45,684 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:45,685 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:45,688 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:45,688 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:45,690 Making a fixed resolution buffer of (('flash', 'velx')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:45,901 Plot image for field ('flash', 'velx') has both positive and negative values. Min = -22627464.000001, Max = 49533188.000002.
yt : [WARNING  ] 2017-03-30 05:58:45,901 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2017-03-30 05:58:46,851 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:46,852 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:46,855 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:46,856 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:46,857 Making a fixed resolution buffer of (('flash', 'vely')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:47,065 Plot image for field ('flash', 'vely') has both positive and negative values. Min = -31364576.000001, Max = 25449842.000001.
yt : [WARNING  ] 2017-03-30 05:58:47,065 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2017-03-30 05:58:47,846 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:47,846 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:47,850 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:47,850 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:47,851 Making a fixed resolution buffer of (('flash', 'gpot')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:48,147 Plot image for field ('flash', 'gpot') has no positive values.  Max = -21768437759996.953125.
yt : [WARNING  ] 2017-03-30 05:58:48,148 Switching to linear colorbar scaling.
yt : [INFO     ] 2017-03-30 05:58:48,409 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:48,409 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:48,413 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:48,413 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:48,415 Making a fixed resolution buffer of (('gas', 'numdens')) 800 by 800
yt : [INFO     ] 2017-03-30 05:58:49,007 Making a fixed resolution buffer of (velx) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:49,115 Making a fixed resolution buffer of (vely) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:49,224 Making a fixed resolution buffer of (velz) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:49,341 Making a fixed resolution buffer of (dens) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:49,453 Making a fixed resolution buffer of (gpot) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:49,509 Making a fixed resolution buffer of (thermal_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:50,912 Making a fixed resolution buffer of (magnetic_energy) 1024 by 1024
yt : [INFO     ] 2017-03-30 05:58:51,231 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:51,232 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:51,235 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:51,235 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:51,237 Making a fixed resolution buffer of (('flash', 'velx')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:51,545 Plot image for field ('flash', 'velx') has both positive and negative values. Min = -147272726.715681, Max = 92215173.895858.
yt : [WARNING  ] 2017-03-30 05:58:51,546 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2017-03-30 05:58:51,838 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:51,839 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:51,842 xlim = -88392327679572443136.000000 2911607672320427556864.000000
yt : [INFO     ] 2017-03-30 05:58:51,843 ylim = -2635563720923045429248.000000 364436279076954832896.000000
yt : [INFO     ] 2017-03-30 05:58:51,844 Making a fixed resolution buffer of (('flash', 'vely')) 800 by 800
yt : [WARNING  ] 2017-03-30 05:58:52,149 Plot image for field ('flash', 'vely') has both positive and negative values. Min = -92118324.667611, Max = 106212398.149214.
yt : [WARNING  ] 2017-03-30 05:58:52,149 Switching to symlog colorbar scaling unless linear scaling is specified later
In [45]:
vcut = 100.0 * km/second


##############################################################################################################################
#                             t = 2 Myr
##############################################################################################################################

slc_frb_vx1_t2_fast = {"info":"this dictionary contains the fast moving gas for the quiver plot."}
slc_frb_vx1_t2_slow = {"info":"this dictionary contains the slow moving gas for the quiver plot."}

vx1_fast = np.zeros_like(slc_frb_vx1_t2["velx"].in_units("km/s")) 
vy1_fast = np.zeros_like(slc_frb_vx1_t2["velx"].in_units("km/s")) 

vx1_slow = np.zeros_like(slc_frb_vx1_t2["velx"].in_units("km/s")) 
vy1_slow = np.zeros_like(slc_frb_vx1_t2["velx"].in_units("km/s")) 

for i in range(np.shape(slc_frb_vx1_t2["velx"])[0]):
    for j in range(np.shape(slc_frb_vx1_t2["velx"])[1]):
        vmag_here = np.sqrt(slc_frb_vx1_t2["velx"][i][j].in_units("km/s")**2 + slc_frb_vx1_t2["velx"][i][j].in_units("km/s")**2)
        if vmag_here > vcut:
            vx1_fast[i][j] = slc_frb_vx1_t2["velx"][i][j].in_units("km/s")
            vy1_fast[i][j] = slc_frb_vx1_t2["vely"][i][j].in_units("km/s")
        else:
            vx1_slow[i][j] = slc_frb_vx1_t2["velx"][i][j].in_units("km/s")
            vy1_slow[i][j] = slc_frb_vx1_t2["vely"][i][j].in_units("km/s")


slc_frb_vx1_t2_fast["velx"] = vx1_fast
slc_frb_vx1_t2_fast["vely"] = vy1_fast

slc_frb_vx1_t2_slow["velx"] = vx1_slow
slc_frb_vx1_t2_slow["vely"] = vy1_slow

##############################################################################################################################
#                             t = 4 Myr
##############################################################################################################################


slc_frb_vx1_t4_fast = {"info":"this dictionary contains the fast moving gas for the quiver plot."}
slc_frb_vx1_t4_slow = {"info":"this dictionary contains the slow moving gas for the quiver plot."}

vx1_fast = np.zeros_like(slc_frb_vx1_t4["velx"].in_units("km/s")) 
vy1_fast = np.zeros_like(slc_frb_vx1_t4["velx"].in_units("km/s")) 

vx1_slow = np.zeros_like(slc_frb_vx1_t4["velx"].in_units("km/s")) 
vy1_slow = np.zeros_like(slc_frb_vx1_t4["velx"].in_units("km/s")) 

for i in range(np.shape(slc_frb_vx1_t4["velx"])[0]):
    for j in range(np.shape(slc_frb_vx1_t4["velx"])[1]):
        vmag_here = np.sqrt(slc_frb_vx1_t4["velx"][i][j].in_units("km/s")**2 + slc_frb_vx1_t4["velx"][i][j].in_units("km/s")**2)
        if vmag_here > vcut:
            vx1_fast[i][j] = slc_frb_vx1_t4["velx"][i][j].in_units("km/s")
            vy1_fast[i][j] = slc_frb_vx1_t4["vely"][i][j].in_units("km/s")
        else:
            vx1_slow[i][j] = slc_frb_vx1_t4["velx"][i][j].in_units("km/s")
            vy1_slow[i][j] = slc_frb_vx1_t4["vely"][i][j].in_units("km/s")


slc_frb_vx1_t4_fast["velx"] = vx1_fast
slc_frb_vx1_t4_fast["vely"] = vy1_fast

slc_frb_vx1_t4_slow["velx"] = vx1_slow
slc_frb_vx1_t4_slow["vely"] = vy1_slow

##############################################################################################################################
#                             t = 6 Myr
##############################################################################################################################


slc_frb_vx1_t6_fast = {"info":"this dictionary contains the fast moving gas for the quiver plot."}
slc_frb_vx1_t6_slow = {"info":"this dictionary contains the slow moving gas for the quiver plot."}

vx1_fast = np.zeros_like(slc_frb_vx1_t6["velx"].in_units("km/s")) 
vy1_fast = np.zeros_like(slc_frb_vx1_t6["velx"].in_units("km/s")) 

vx1_slow = np.zeros_like(slc_frb_vx1_t6["velx"].in_units("km/s")) 
vy1_slow = np.zeros_like(slc_frb_vx1_t6["velx"].in_units("km/s")) 

for i in range(np.shape(slc_frb_vx1_t6["velx"])[0]):
    for j in range(np.shape(slc_frb_vx1_t6["velx"])[1]):
        vmag_here = np.sqrt(slc_frb_vx1_t6["velx"][i][j].in_units("km/s")**2 + slc_frb_vx1_t6["velx"][i][j].in_units("km/s")**2)
        if vmag_here > vcut:
            vx1_fast[i][j] = slc_frb_vx1_t6["velx"][i][j].in_units("km/s")
            vy1_fast[i][j] = slc_frb_vx1_t6["vely"][i][j].in_units("km/s")
        else:
            vx1_slow[i][j] = slc_frb_vx1_t6["velx"][i][j].in_units("km/s")
            vy1_slow[i][j] = slc_frb_vx1_t6["vely"][i][j].in_units("km/s")


slc_frb_vx1_t6_fast["velx"] = vx1_fast
slc_frb_vx1_t6_fast["vely"] = vy1_fast

slc_frb_vx1_t6_slow["velx"] = vx1_slow
slc_frb_vx1_t6_slow["vely"] = vy1_slow
yt : [INFO     ] 2017-03-30 05:58:52,525 Making a fixed resolution buffer of (velx) 16 by 16
yt : [INFO     ] 2017-03-30 05:58:52,538 Making a fixed resolution buffer of (vely) 16 by 16
yt : [INFO     ] 2017-03-30 05:58:53,126 Making a fixed resolution buffer of (velx) 16 by 16
yt : [INFO     ] 2017-03-30 05:58:53,139 Making a fixed resolution buffer of (vely) 16 by 16
yt : [INFO     ] 2017-03-30 05:58:53,734 Making a fixed resolution buffer of (velx) 16 by 16
yt : [INFO     ] 2017-03-30 05:58:53,750 Making a fixed resolution buffer of (vely) 16 by 16
In [49]:
fig = plt.figure(figsize=(14,8))

asp_ratio = 14.0 / 8.0

plt_xsize = 0.28
plt_ysize = 0.28 * asp_ratio

y00 = 0.06

scale_slow_vecs = 200
scale_fast_vecs = 1000

##############################################################################################################################
#                             t = 2 Myr
##############################################################################################################################

ax0 = fig.add_axes([0.05, y00, plt_xsize, plt_ysize], frameon=True)

pot_pos = np.zeros_like(total_energy_t2)

for ii in range(np.shape(total_energy_t2)[0]):
    for jj in range(np.shape(total_energy_t2)[1]): 
        if abs(float(total_energy_t2[ii][jj])) > 0:
            pot_pos[ii][jj] = total_energy_t2[ii][jj]
        else:
            pot_pos[ii][jj] = np.nan
            
pot_neg = np.zeros_like(total_energy_t2)
            
for ii in range(np.shape(total_energy_t2)[0]):
    for jj in range(np.shape(total_energy_t2)[1]):
        if float(total_energy_t2[ii][jj]) < 0:
            pot_neg[ii][jj] = abs(float(total_energy_t2[ii][jj]))
        else:
            pot_neg[ii][jj] = np.nan

my_cmap = copy.copy(plt.cm.get_cmap('Reds'))
my_cmap.set_bad(alpha=0) # set how the colormap handles 'bad' values

my_cmap2 = copy.copy(plt.cm.get_cmap('Blues'))
my_cmap2.set_bad(alpha=0) # set how the colormap handles 'bad' values

vmin = 1.0e-15
vmax = 1.0e-7

cax00_pos = ax0.imshow(pot_pos, origin='center', extent = extent_t2, aspect='equal' , cmap=my_cmap,  vmin=vmin, vmax=vmax,  norm=LogNorm())
cax00_neg = ax0.imshow(pot_neg, origin='center', extent = extent_t2, aspect='equal' , cmap=my_cmap2, vmin=vmin, vmax=vmax,  norm=LogNorm())

ax0.set_title("t=2 Myr", fontsize=25)

ax0.set_aspect('auto')

ax0.set_xlabel('x [pc]', {'fontsize': 15})
for tick in ax0.xaxis.get_major_ticks():
    tick.label.set_fontsize(12) 

ax0.set_ylabel('y [pc]', {'fontsize': 15}, labelpad=0)
for tick in ax0.yaxis.get_major_ticks():
    tick.label.set_fontsize(12) 
    
#ax0.contour(slc_frb["numdens"].value, extent = extent, levels=[1000],  alpha=1, linewidths=1,   colors='k', data_source=my_cloud)

ax0.contour(slc_frb_ndens_t2["numdens"].value, extent = extent_t2, levels=[100],  alpha=1, linewidths=2,   colors='k')
ax0.contour(slc_frb_ndens_t2["numdens"].value, extent = extent_t2, levels=[10],   alpha=1, linewidths=2.0, linestyles='dotted', colors='k')


# Plot the velocity vectors
# --------------------------------------------------------------------------------------
#ax0.quiver(slc_frb_vx1_t2["x"].in_units("pc").value, slc_frb_vx1_t2["y"].in_units("pc").value, 
#           slc_frb_vx1_t2["velx"].in_units("km/s").value, slc_frb_vy1_t2["vely"].in_units("km/s").value, units='inches', scale=100)

# Plot the velocity vectors - white if velocity is above 100 km/s
# --------------------------------------------------------------------------------------
ax0.quiver(slc_frb_vx1_t2["x"].in_units("pc").value, slc_frb_vx1_t2["y"].in_units("pc").value, slc_frb_vx1_t2_fast["velx"].value, 
           slc_frb_vx1_t2_fast["vely"].value, units='inches', scale=scale_fast_vecs, color="white", width=0.02)

# Plot the velocity vectors - black if velocity is below 100 km/s
# --------------------------------------------------------------------------------------
ax0.quiver(slc_frb_vx1_t2["x"].in_units("pc").value, slc_frb_vx1_t2["y"].in_units("pc").value, 
           slc_frb_vx1_t2_slow["velx"].value, slc_frb_vx1_t2_slow["vely"].value, units='inches', scale=scale_slow_vecs, color="black")


##############################################################################################################################
#                             t = 4 Myr
##############################################################################################################################

ax1 = fig.add_axes([0.038*2 + 0.005 + plt_xsize , y00, plt_xsize, plt_ysize], frameon=True)


pot_pos = np.zeros_like(total_energy_t4)

for ii in range(np.shape(total_energy_t4)[0]):
    for jj in range(np.shape(total_energy_t4)[1]): 
        if abs(float(total_energy_t4[ii][jj])) > 0:
            pot_pos[ii][jj] = total_energy_t4[ii][jj]
        else:
            pot_pos[ii][jj] = np.nan
            
pot_neg = np.zeros_like(total_energy_t4)
            
for ii in range(np.shape(total_energy_t4)[0]):
    for jj in range(np.shape(total_energy_t4)[1]):
        if float(total_energy_t4[ii][jj]) < 0:
            pot_neg[ii][jj] = abs(float(total_energy_t4[ii][jj]))
        else:
            pot_neg[ii][jj] = np.nan

my_cmap = copy.copy(plt.cm.get_cmap('Reds'))
my_cmap.set_bad(alpha=0) # set how the colormap handles 'bad' values

my_cmap2 = copy.copy(plt.cm.get_cmap('Blues'))
my_cmap2.set_bad(alpha=0) # set how the colormap handles 'bad' values

vmin = 1.0e-15
vmax = 1.0e-7

cax00_pos = ax1.imshow(pot_pos, origin='center', extent = extent_t4, aspect='equal' , cmap=my_cmap,  vmin=vmin, vmax=vmax,  norm=LogNorm())
cax00_neg = ax1.imshow(pot_neg, origin='center', extent = extent_t4, aspect='equal' , cmap=my_cmap2, vmin=vmin, vmax=vmax,  norm=LogNorm())

ax1.set_title("t=4 Myr", fontsize=25)

ax1.set_aspect('auto')

ax1.set_xlabel('x [pc]', {'fontsize': 15})
for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(12) 

ax1.set_ylabel('y [pc]', {'fontsize': 0})
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(12) 

ax1.contour(slc_frb_ndens_t4["numdens"].value, extent = extent_t4, levels=[100],  alpha=1, linewidths=2,   colors='k', data_source=cloud_t4)
ax1.contour(slc_frb_ndens_t4["numdens"].value, extent = extent_t4, levels=[10],   alpha=1, linewidths=2.0, linestyles='dotted', colors='k')

# Plot the velocity vectors
# --------------------------------------------------------------------------------------
#ax1.quiver(slc_frb_vx1_t4["x"].in_units("pc").value, slc_frb_vx1_t4["y"].in_units("pc").value, 
#           slc_frb_vx1_t4["velx"].in_units("km/s").value, slc_frb_vy1_t4["vely"].in_units("km/s").value, units='inches', scale=100)


# Plot the velocity vectors - white if velocity is above 100 km/s
# --------------------------------------------------------------------------------------
ax1.quiver(slc_frb_vx1_t4["x"].in_units("pc").value, slc_frb_vx1_t4["y"].in_units("pc").value, slc_frb_vx1_t4_fast["velx"].value, 
           slc_frb_vx1_t4_fast["vely"].value, units='inches', scale=scale_fast_vecs, color="white", width=0.02)

# Plot the velocity vectors - black if velocity is below 100 km/s
# --------------------------------------------------------------------------------------
ax1.quiver(slc_frb_vx1_t4["x"].in_units("pc").value, slc_frb_vx1_t4["y"].in_units("pc").value, 
           slc_frb_vx1_t4_slow["velx"].value, slc_frb_vx1_t4_slow["vely"].value, units='inches', scale=scale_slow_vecs, color="black")


##############################################################################################################################
#                             t = 6 Myr
##############################################################################################################################

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\#
# Black Box behind the velocity vector                     #
ax5 = fig.add_axes([0.9, 0.785, 0.05, 0.015], frameon=True)  #
ax5.axes.get_xaxis().set_visible(False)                    #
ax5.axes.get_yaxis().set_visible(False)                    #
ax5.patch.set_facecolor("black")                           #
# /////////////////////////////////////////////////////////#

ax2 = fig.add_axes([0.04*2 + 0.03 + plt_xsize*2 , y00, plt_xsize, plt_ysize], frameon=True)


pot_pos = np.zeros_like(total_energy_t6)

for ii in range(np.shape(total_energy_t6)[0]):
    for jj in range(np.shape(total_energy_t6)[1]): 
        if abs(float(total_energy_t6[ii][jj])) > 0:
            pot_pos[ii][jj] = total_energy_t6[ii][jj]
        else:
            pot_pos[ii][jj] = np.nan
            
pot_neg = np.zeros_like(total_energy_t6)
            
for ii in range(np.shape(total_energy_t6)[0]):
    for jj in range(np.shape(total_energy_t6)[1]):
        if float(total_energy_t6[ii][jj]) < 0:
            pot_neg[ii][jj] = abs(float(total_energy_t6[ii][jj]))
        else:
            pot_neg[ii][jj] = np.nan

my_cmap = copy.copy(plt.cm.get_cmap('Reds'))
my_cmap.set_bad(alpha=0) # set how the colormap handles 'bad' values

my_cmap2 = copy.copy(plt.cm.get_cmap('Blues'))
my_cmap2.set_bad(alpha=0) # set how the colormap handles 'bad' values

vmin = 1.0e-15
vmax = 1.0e-7

cax00_pos = ax2.imshow(pot_pos, origin='center', extent = extent_t6, aspect='equal' , cmap=my_cmap,  vmin=vmin, vmax=vmax,  norm=LogNorm())
cax00_neg = ax2.imshow(pot_neg, origin='center', extent = extent_t6, aspect='equal' , cmap=my_cmap2, vmin=vmin, vmax=vmax,  norm=LogNorm())

ax2.set_title("t=6 Myr", fontsize=25)

ax2.set_aspect('auto')

ax2.set_xlabel('x [pc]', {'fontsize': 15})
for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(12) 
    
ax2.tick_params(axis='y', which='both', labelleft='on', labelright='on')
ax2.yaxis.set_label_position("right")

ax2.set_ylabel('y [pc]', {'fontsize': 15})
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(12) 

ax2.contour(slc_frb_ndens_t6["numdens"], extent = extent_t6, levels=[100],  alpha=1, linewidths=2,   colors='k')
ax2.contour(slc_frb_ndens_t6["numdens"], extent = extent_t6, levels=[10],   alpha=1, linewidths=2.0, linestyles='dotted', colors='k')
#ax2.contour(slc_frb_ndens_t6["numdens"], extent = extent_t6, levels=[1],    alpha=1,  linewidths=1.0, linestyles='dotted', colors='k')

# Plot the velocity vectors
# --------------------------------------------------------------------------------------
#vax = ax2.quiver(slc_frb_vx1_t6["x"].in_units("pc").value, slc_frb_vx1_t6["y"].in_units("pc").value, 
#           slc_frb_vx1_t6["velx"].in_units("km/s").value, slc_frb_vy1_t6["vely"].in_units("km/s").value, units='inches', scale=100)


# Plot the velocity vectors - white if velocity is above 100 km/s
# --------------------------------------------------------------------------------------
vax = ax2.quiver(slc_frb_vx1_t6["x"].in_units("pc").value, slc_frb_vx1_t6["y"].in_units("pc").value, slc_frb_vx1_t6_fast["velx"].value, 
           slc_frb_vx1_t6_fast["vely"].value, units='inches', scale=scale_fast_vecs, color="white", width=0.02)

# Plot the velocity vectors - black if velocity is below 100 km/s
# --------------------------------------------------------------------------------------
vax1 = ax2.quiver(slc_frb_vx1_t6["x"].in_units("pc").value, slc_frb_vx1_t6["y"].in_units("pc").value, 
           slc_frb_vx1_t6_slow["velx"].value, slc_frb_vx1_t6_slow["vely"].value, units='inches', scale=scale_slow_vecs, color="black")


##################################################################################################################
#  COLORBAR
##################################################################################################################

cax_red = fig.add_axes([0.3, plt_ysize + 0.33, 0.55, 0.04])

cbar = fig.colorbar(cax00_pos, ticks=[1.0e-14, 1.0e-12, 1.0e-10, 1.0e-8], orientation='horizontal', cax=cax_red)
cbar.ax.set_xticklabels(['10$^{-14}$','10$^{-12}$', '10$^{-10}$', '10$^{-8}$'], fontsize=15)
cbar.ax.invert_xaxis()
cbar.set_label("Kinetic, Thermal, Magnetic dominated", fontsize=15)

cax_blue = fig.add_axes([0.3, plt_ysize + 0.20, 0.55, 0.04])

cbar = fig.colorbar(cax00_neg, ticks=[1.0e-14, 1.0e-12, 1.0e-10, 1.0e-8], orientation='horizontal', cax=cax_blue)
cbar.ax.set_xticklabels(['-10$^{-14}$','-10$^{-12}$', '-10$^{-10}$', '-10$^{-8}$'], fontsize=15) 
cbar.ax.invert_xaxis()
cbar.set_label("Self-gravity dominated", fontsize=15)


########################################################################################################################
# 3D Figure
########################################################################################################################

# my position in z now.
znow = CM_t4[2]/ppc

# Now I have to find out in which indez in the zrange am I.
indexz = int((znow - minzUG) *znump / zrangeUG)

stop = 0
i    = 0
while stop == 0:
    i +=1
    max_here = np.max(UG["numdens"][-1*i,:,indexz])/ppc
    if max_here != 0:
        stop = 1
    if i == np.shape(UG["x"])[1]:
        stop = 1
xhere = UG["x"][-1*i,0,indexz]/ppc

stop = 0
i    = 0
while stop == 0:
    i +=1
    max_here = np.max(UG["numdens"][:,i,indexz])/ppc
    if max_here != 0:
        stop = 1
    if i == np.shape(UG["x"])[0]:
        stop = 1
yhere = UG["y"][0,i,indexz]/ppc


ax3 = fig.add_axes([0.05,0.63,0.20,0.20*asp_ratio], projection='3d')

#Xp, Yp = np.mgrid[extent_t4[0]:extent_t4[1]:2j, extent_t4[2]:extent_t4[3]:2j]
#ax3.plot_surface(Xp, Yp, [znow,znow], alpha=1.0, color="grey", shade=True)


ax3.plot_trisurf(verts[:,0], verts[:,1], faces, verts[:,2],
                color='red', lw=0, antialiased=False, shade=True)

if znow >= minzUG and znow <= maxzUG:
    # line 1
    ax3.plot([minxUG, maxxUG], [minyUG, minyUG], [znow, znow], "-k", linewidth=2)
    # line 2
    ax3.plot([maxxUG, maxxUG], [minyUG, maxyUG], [znow, znow], "-k", linewidth=2)
    # line 3
    ax3.plot([xhere, maxxUG], [maxyUG, maxyUG], [znow, znow], "-k", linewidth=2)
    # line 4
    ax3.plot([minxUG,minxUG], [minyUG, yhere], [znow, znow], "-k", linewidth=2)

ax3.view_init(elev=35, azim=-58)    

#ax3.set_xlabel("X", fontsize = 15)
#ax3.set_ylabel("Y", fontsize = 15)
ax3.set_zlabel("Z", fontsize = 15)



##########################################################
# Add an arrow fr the observed line of sight.
##########################################################
ax4= fig.add_axes([0.1, 0.9, 0.1, 0.1], frameon=False)
ax4.arrow( 0.5, 0.85, 0.0, -0.4, fc="k", ec="k", head_width=0.04, head_length=0.1 )

ax4.axes.get_xaxis().set_visible(False)
ax4.axes.get_yaxis().set_visible(False)

ax4.text(0, 0.5, "Line of\nsight", fontsize=11)


#############################################
# Title and other labels
#############################################

fig.suptitle("Cloud %s  Slice at d$_{CM}$(z)=%2.1f pc \nTotal energy density [erg cm$^{-3}$]. " %(cloud_plt_name, dcm), fontsize=20)
#ax.set_title("Total Energy in and around cloud %s, slice z = %2.1f pc" %(Cloud_name, i*stepsize - rad), fontsize=25)

ax2.quiverkey(vax1, 0.9, 1.59,  50,     '50 km/s',fontproperties={'weight': 'bold'})

ax2.quiverkey(vax,  0.91, 1.49,  500,  '500 km/s', fontproperties={'weight': 'bold'})


#############################################
# Add The lines of the Contours

ax5= fig.add_axes([0.90, 0.6, 0.1, 0.15], frameon=False)

ax5.axes.get_xaxis().set_visible(False)
ax5.axes.get_yaxis().set_visible(False)

ax5.text(0, 0.6, "n=100 cm$^{-3}$ \nContour", fontsize=12)
ax5.plot([0, 0.6], [0.92, 0.92], "-k", linewidth=3)

ax5.text(0, 0.0, "n=10  cm$^{-3}$ \nContour", fontsize=12)
ax5.plot([0, 0.6], [0.38, 0.38], ":k", linewidth=3)

ax5.set_ylim(0, 1)
ax5.set_xlim(0,1)

fig.show()

save_dir = "/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/RedVsBlue/%s" %Cloud_name

#fig.savefig("%s/virial_plot_%s.pdf"%(save_dir, Cloud_name), format="pdf", dpi=100)
if Cloud_name == "M3e3":
    fig.savefig("%s/virial_plot_%s.png"%(save_dir, Cloud_name), format="png", dpi=100)
    print("%s/virial_plot_%s.png"%(save_dir, Cloud_name))
else:
    fig.savefig("%s/virial_plot_%s.pdf"%(save_dir, Cloud_name), format="pdf", dpi=100)
    print("%s/virial_plot_%s.pdf"%(save_dir, Cloud_name))
/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/RedVsBlue/M3e3/virial_plot_M3e3.png
In [47]:
#fig.savefig("%s/virial_plot2_%s.png"%(save_dir, Cloud_name), format="png", dpi=200)
In [48]:
print "Hola"
Hola