In [ ]:

    

import yt
import numpy as np
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt

#ds = yt.load('/d/d5/ychen/2015_production_runs/0529_h1_00/data/MHD_Jet_hdf5_plt_cnt_0600')
ds = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0910')


    

In [ ]:

    

from synchrotron.yt_synchrotron_emissivity import _jet_volume_fraction
ds.add_field(('gas', 'jet_volume_fraction'), function=_jet_volume_fraction,
             display_name="Jet Volume Fraction", sampling_type='cell')
plot = yt.SlicePlot(ds, 'x', 'jet_volume_fraction', center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])

plot.set_zlim('jet_volume_fraction', 1E-2, 1E0)
plot.show()


    

In [ ]:

    

plot.set_zlim('jet_volume_fraction', 5E-2, 1E0)
plot.show()


    

In [ ]:

    

ad = ds.all_data()
null = plt.hist(np.log10(ad['jet_volume_fraction']), range=(-2,0.00), bins=100, log=True)


    

In [ ]:

    

plt.scatter(np.log10(ad['flash', 'jet ']), np.log10(ad['jet_volume_fraction']), s=1, lw=0)


    

In [ ]:

    

plot.set_zlim('jet_volume_fraction', 1E-1, 1E0)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', 'jet ', center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])
plot.set_zlim('jet ', 1E-10, 1E0)
plot.set_log('jet ', True)
#plot.set_cmap('jet ', 'RdBu_r')
#plot.set_buff_size((200,200))


    

In [ ]:

    

data = ds.all_data()

from yt.units import g, cm, Kelvin
mH = yt.utilities.physical_constants.mass_hydrogen
k  = yt.utilities.physical_constants.boltzmann_constant

rhoCore = ds.parameters['sim_rhocore']*g/cm**3
rCore   = ds.parameters['sim_rcore']*cm
densitybeta = ds.parameters['sim_densitybeta']
Tout    = ds.parameters['sim_tout']*Kelvin
Tcore   = ds.parameters['sim_tcore']*Kelvin
rCoreT  = ds.parameters['sim_rcoret']*cm
gammaICM= ds.parameters['sim_gammaicm']
mu      = ds.parameters['sim_mu']

r = data['index', 'spherical_radius']

density0 = rhoCore*(1.0 + (r/rCore)**2)**(-1.5*densitybeta)
T0 = Tout*(1.0+(r/rCoreT)**3)/(Tout/Tcore+(r/rCoreT)**3)
P0 = density0/mu/mH*k*T0
icm_mass_fraction = 1.0 - data['flash', 'jet ']
P = data['gas', 'pressure']
density = data['gas', 'density']

icm_volume_fraction = (P0/P)**(1/gammaICM)*icm_mass_fraction*density/density0

#icm_volume_fraction = np.where(icm_volume_fraction < 1.0, icm_volume_fraction, 1.0)


    

In [ ]:

    

null = plt.hist(np.log10(1-icm_volume_fraction), range=(-3,0.05), bins=100, log=True)


    

In [ ]:

    

data = ds.all_data()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)
null = ax1.hist(np.log10(data['flash', 'jet ']), range=(-10.02,0), bins=100, log=True)
null = ax2.hist(np.log10(data['flash', 'jet ']), range=(-9.9,0), bins=100, log=False)


    

In [ ]:

    

#mask = data['flash', 'jet '] < 0.0
max(data['flash', 'jet '])


    

In [ ]:

    

plt.scatter(np.log10(data['flash', 'jet ']), np.log10(1-icm_volume_fraction), s=1, lw=0)


    

In [ ]:

    

mask = icm_mass_fraction == 1.00
icm_mass_fraction[mask]


    

In [ ]:

    

np.where(icm_volume_fraction > 1.0)


    

In [ ]: