In [ ]:

    

%matplotlib inline
import numpy as np
import os
import util
import multiprocessing
import matplotlib
matplotlib.rcParams['savefig.dpi'] = 150
import matplotlib.pyplot as plt
import yt
import logging
from yt.utilities.file_handler import HDF5FileHandler
from yt.funcs import mylog
#logging.getLogger('yt').setLevel(logging.WARNING)
print(yt.__version__)

from particle_filters import *

e  = yt.utilities.physical_constants.elementary_charge #4.803E-10 esu
me = yt.utilities.physical_constants.mass_electron #9.109E-28
c  = yt.utilities.physical_constants.speed_of_light #2.998E10

def setup_part_file(ds):
    filename = os.path.join(ds.directory,ds.basename)
    ds._particle_handle = HDF5FileHandler(filename.replace('plt_cnt', 'part')+'_updated')
    ds.particle_filename = filename.replace('plt_cnt', 'part')+'_updated'
    mylog.info('Changed particle files to:' + ds.particle_filename)


    

In [ ]:

    

#fname = '/home/ychen/d9/FLASH4/stampede/1022_L45_M10_b1_h1_10Myr/MHD_Jet_10Myr_hdf5_plt_cnt_0630'
#fname = '/home/ychen/d9/FLASH4/2015_production_runs/0529_L45_M10_b1_h1/MHD_Jet_hdf5_plt_cnt_0800'
#fname = '/home/ychen/d9/FLASH4/2016_test/1106_L45_M3_b01_h1/MHD_Jet_hdf5_plt_cnt_0002'
#fname = '/home/ychen/Mount/stampede/1106_L45_M3_b1_h1/MHD_Jet_hdf5_plt_cnt_0059'
#fname = '/home/ychen/d9/FLASH4/2016_test/1028_metal_test/MHD_Jet_hdf5_plt_cnt_0002'
fname = '/d/d5/ychen/2015_production_runs/1110_h0_00/data/MHD_Jet_hdf5_plt_cnt_0100'
ds = yt.load(fname)
setup_part_file(ds)

ad = ds.all_data()


    

In [ ]:

    

#plt.hist(ad['particle_dtau'], bins=100, range=(-1E-5,8E-4))
mask = ad['particle_den1'] < 0
print(ad['particle_dtau'][mask])
#gamc = (ad['particle_dens']/ad['particle_den1'])**(1./3.)/ ad['particle_dtau']
#dens = ad['particle_dens']
#den1 = ad['particle_den1']
#print(gamc.max())
#print(gamc.min())
#print((dens)[np.isnan(gamc)])


    

In [ ]:

    

ptype = 'jet'
ds.add_particle_filter(ptype)
tag = ad[ptype, 'particle_tag']
#filter = np.logical_and((tag % 1 == 0), (ad['all', 'particle_shok']==0))
filter = tag % 5 == 0
y = ad[ptype, 'particle_position_y'][filter]/3.08567758E21
z = ad[ptype, 'particle_position_z'][filter]/3.08567758E21
age = (ds.current_time.v - np.abs(ad[ptype, 'particle_tadd'][filter]))/ds.current_time.v
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-20,20)
ax.set_ylim(-40,40)

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=age,linewidth=0,cmap='arbre_r', vmin=0,vmax=1,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('normalized age')
plt.show()


    

In [ ]:

    

field = 'particle_nuc_dtau'

ptype = 'lobe'
ds.add_particle_filter(ptype)
tag = ad[ptype, 'particle_tag']
#filter = np.logical_and((tag % 1 == 0), (ad['all', 'particle_shok']==0))
filter = tag % 5 == 0
x = ad[ptype, 'particle_position_x'][filter]/3.08567758E21
z = ad[ptype, 'particle_position_z'][filter]/3.08567758E21

if 'dtau' in field:
    gamc = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1'])**(1./3.) \
           / ad[ptype, 'particle_dtau']


B = np.sqrt(ad[(ptype, 'particle_magx')][filter]**2
           +ad[(ptype, 'particle_magy')][filter]**2
           +ad[(ptype, 'particle_magz')][filter]**2)*np.sqrt(4.0*np.pi)
B = ad.apply_units(B, 'gauss')
# Cutoff frequency
fdata = np.log10(3.0*gamc[filter]**2*e*B/(4.0*np.pi*me*c))
vmin=7.5; vmax=12; cmap='algae'
cblabel=u'log $\\nu_c$'

if 'dtau' in field:
    cblabel=cblabel+' (dtau)'


age = (ds.current_time.v - np.abs(ad[ptype, 'particle_tadd'][filter]))/ds.current_time.v
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-20,20)
ax.set_xlabel('x')
ax.set_ylim(-40,40)
ax.set_ylabel('z')

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(x,z,s=1,c=age,linewidth=0,cmap='arbre_r', vmin=0,vmax=1,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label(cblabel)
plt.show()


    

In [ ]:

    

ptype = 'lobe'
ds.add_particle_filter(ptype)
tag = ad[ptype, 'particle_tag']
#filter = np.logical_and((tag % 1 == 0), (ad['all', 'particle_shok']==0))
filter = tag % 5 == 0
y = ad[ptype, 'particle_position_y'][filter]/3.08567758E21
z = ad[ptype, 'particle_position_z'][filter]/3.08567758E21
age = (ds.current_time.v - np.abs(ad[ptype, 'particle_tadd'][filter]))/ds.current_time.v
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-20,20)
ax.set_ylim(-40,40)

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=age,linewidth=0,cmap='arbre_r', vmin=0,vmax=1,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('normalized age')
plt.show()


    

In [ ]:

    

ptype = 'lobe'
print ds.add_particle_filter(ptype)
tag = ad[ptype, 'particle_tag']
filter = (tag % 5 == 0)
y = ad[ptype, 'particle_position_y'][filter]/3.08567758E21
z = ad[ptype, 'particle_position_z'][filter]/3.08567758E21
c = np.log10(ad[ptype, 'particle_gamc'][filter])
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-20,20)
ax.set_ylim(-40,40)

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=c,linewidth=0,cmap='arbre', vmin=2,vmax=5,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('gamc')
plt.show()


    

In [ ]:

    

ds.add_particle_filter('test')


    

In [ ]:

    

('particle_shok') in ds.derived_field_list


    

In [ ]:

    

p = yt.SlicePlot(ds,'x', 'density', width=((20,'kpc'),(20,'kpc')), )
p.annotate_particles((2, 'kpc'), ptype="metal")
p.annotate_grids(linewidth=0.1)
#p = yt.ParticleProjectionPlot(ds, 'x', 'particle_tag', width=((500,'kpc'),(1000,'kpc')), depth=(500,'kpc'))
#p.set_log('particle_posz', False)
p.show()


    

In [ ]:

    

ptPerGrid = np.zeros(ds.index.num_grids, dtype=int)
ptPerMass = np.zeros(ds.index.num_grids)
ptPerVolume = np.zeros(ds.index.num_grids)
level = np.zeros(ds.index.num_grids, dtype=int)
rr = np.zeros(ds.index.num_grids)
for i in range(ds.index.num_grids):
    grid = ds.index.grids[i]
    grid.set_field_parameter('center', ds.arr([0,0,0]))
    ptPerGrid[i] = grid.NumberOfParticles
    ptPerMass[i] = grid.NumberOfParticles/np.sum(grid['cell_mass'])
    ptPerVolume[i] = grid.NumberOfParticles/np.sum(grid['cell_volume'].in_units('kpc**3'))
    level[i] = grid.Level
    rr[i] = grid['radius'][4,4,4].in_units('kpc')


    

In [ ]:

    

#print ptPerMass
#plt.scatter(level, ptPerGrid, s=1, lw=0)
#plt.scatter(level, ptPerMass*1E42, s=1, lw=0, c='r')
plt.scatter(rr, ptPerVolume, s=1, lw=0, c='r')
plt.xlim(0,500)
plt.xlabel('r (kpc)')
plt.ylabel('# particles per kpc^3')


    

In [ ]:

    

kpc=yt.units.kpc.in_units('cm')
ad = ds.all_data()
sphere=ds.sphere(np.array([0,0,0]), 150*kpc)


    

In [ ]:

    

print len(ad['particle_posz'])
numParticle = len(sphere['particle_posz'])
print numParticle
volume = 3./4.*np.pi*100**3
print volume
print numParticle/volume
print 'mass', sum(sphere['cell_mass'])
print 'volume', sum(sphere['cell_volume'])


    

In [ ]:

    

ad = ds.all_data()
ad.set_field_parameter('center', ds.arr([0,0,0]))
kpc = yt.units.kpc
plt.hist(ad[('io', 'particle_radius')]/kpc.in_units('cm').v, bins=100
)#, cumulative=True)
#plt.scatter(ad['particle_tag'], ad['particle_blk'], c=ad['particle_proc'], s=1, lw=0)
plt.xlabel('kpc')


    

In [ ]:

    

ma = int(max(box['particle_tag']))
alltag = np.array(ad['particle_tag'], dtype=int)


    

In [ ]:

    

for i in range(1,ma+1):
    if i != alltag[i-1]:
        print i, alltag[i-1], ad['particle_posx'][i]


    

In [ ]:

    

print alltag


    

In [ ]:

    

ad = ds.all_data()
tag = ad['all', 'particle_tag']
filter = np.logical_and((tag % 2 == 0), (ad['all', 'particle_shok']==0))
x = ad['all', 'particle_tadd'][filter]/3.15569E13
y = np.log10(ad['all', 'particle_dens'][filter]/ad['all', 'particle_den0'][filter])
#y = np.sqrt(ad['all', 'particle_posx'][filter]**2+ad['all', 'particle_posy'][filter]**2)
#pos = ad['all', 'particle_position'][filter]
#y = ds.find_field_values_at_points('jet ', pos)
z = np.log10(ad['all', 'particle_gamc'][filter])

#B = np.log10(ad['all', 'particle_magx'][filter]**2\
#             +ad['all', 'particle_magy'][filter]**2\
#             +ad['all', 'particle_magz'][filter]**2)

fig=plt.figure(figsize=(12,8))
sc=plt.scatter(x,y,c=z,s=1,linewidth=0,vmin=2,vmax=5)
plt.xlim(-0.1,10)
plt.xlabel('t (Myr)')
cb=plt.colorbar(sc)


    

In [ ]:

    

print z.min()
print z.max()


    

In [ ]:

    

ad = ds.all_data()
tag = ad['all', 'particle_tag']
filter = np.logical_and((tag % 5 == 0), (ad['all', 'particle_shok']==0))
y = ad['all', 'particle_position_y'][filter]/3.08567758E21
z = ad['all', 'particle_position_z'][filter]/3.08567758E21
age = (ds.current_time.v - np.abs(ad['all', 'particle_tadd'][filter]))/ds.current_time.v
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-30,30)
ax.set_ylim(-60,60)

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=age,linewidth=0,cmap='arbre_r', vmin=0,vmax=1,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('normalized age')
plt.show()


    

In [ ]:

    

from tools import calcDen0

ad = ds.all_data()
tag = ad['all', 'particle_tag']
filter = np.logical_and((ad['all', 'particle_shok']==0), (ad['all', 'particle_type']==1.0))
den0core = calcDen0(ad, ptype='io')
den0 = ad['all', 'particle_den0'][filter]
n,bins,patches = plt.hist(den0, bins=100, range=(1.7315E-26,1.733E-26))
for delta in [5E-31, 1E-30, 3E-30]:
    den0min = den0core - delta #1.732E-26
    den0max = den0core + delta #1.7325E-26

    plt.vlines(den0min, 0, max(n))
    plt.vlines(den0max, 0, max(n))
    
    corefil = np.logical_and((ad['all', 'particle_den0']>den0min), 
                     (ad['all', 'particle_den0']<den0max))
    jetfil = np.logical_and((ad['all', 'particle_shok']==0), (ad['all', 'particle_type']==1.0))
    fil = np.logical_and(corefil, jetfil)

    njetp = len(ad['all', 'particle_dens'][jetfil]) 
    njetpc = len(ad['all', 'particle_dens'][fil])

    print('Number of Particles in the jet core:', njetpc)
    print('Number of Particles:', njetp)
    print('Fraction:', float(njetpc)/float(njetp))
#plt.vlines(den0core, 0, max(n))
#plt.semilogy()


    

In [ ]:

    

filter = np.logical_and((ad['all', 'particle_shok']==0), (ad['all', 'particle_type']==1.0))
ds.add_particle_filter('jet')
den0core = calcDen0(ad, ptype='jet')
den0 = ad['all', 'particle_den0'][filter]
print(np.mean(den0))
print(np.mean(den0core))
n,bins,patches = plt.hist(den0, bins=100, range=(1.69E-26, 1.6905E-26), alpha=0.3
n,bins,patches = plt.hist(den0core, bins=100, range=(1.69E-26, 1.6905E-26), alpha=0.3)
ad.keys()


    

In [ ]:

    

corefil = np.logical_and((ad['all', 'particle_den0']>den0min), 
                     (ad['all', 'particle_den0']<den0max))
jetfil = (ad['all', 'particle_shok']==0)
fil = np.logical_and(corefil, jetfil)

njetp = len(ad['all', 'particle_dens'][jetfil]) 
njetpc = len(ad['all', 'particle_dens'][fil])

print njetpc
print njetp
print float(njetpc)/float(njetp)


    

In [ ]:

    

ad = ds.all_data()
tag = ad['all', 'particle_tag']
filter = np.logical_and((tag % 10 == 0), fil)
y = ad['all', 'particle_position_y'][filter]/3.08567758E21
z = ad['all', 'particle_position_z'][filter]/3.08567758E21
c = np.log10(ad['all', 'particle_gamc'][filter])
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(5,8))
ax=fig.add_subplot(111)
ax.set_xlim(-30,30)
ax.set_ylim(-60,60)

ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=c,linewidth=0,cmap='arbre', vmin=2,vmax=5,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('gamc')
plt.show()


    

In [ ]:

    

#ds.add_particle_filter('xcenter')
#ds.add_deposited_particle_field(('jetp', 'particle_gamc'), 'cic')
ds.add_deposited_particle_field(('jetp', 'particle_gamc'), 'cic')
#ds.add_deposited_particle_field(('jetp', 'particle_ones'), 'cic')
def _weighted_cic(field, data):
    return data['deposit', 'jetp_cic_gamc']/data['deposit', 'jetp_count']
ds.add_field(('deposit', 'jetp_weighted_cic_gamc'), function=_weighted_cic)


    

In [ ]:

    

#reg = ds.region(center=[ 0.,  0.,  0.], left_edge=[ -5.0e21,  -1.2e+22,  -2.4e+22], right_edge=[  5.0e21,   1.2e+22,   2.4e+22])
#print reg['all', 'particle_position_x']
#print ad.left_edge
#print ad.right_edge

fields = 'jetp_weighted_cic_gamc'
slice = yt.SlicePlot(ds, 'x', fields=fields)#, center=(0.0,0.0,6.17E22), width=(40,'kpc'))
slice.zoom(12)
slice.set_log(fields, True)
slice.set_zlim(fields, 1E2, 1E5)
slice.annotate_grids()
slice.set_origin('native')
slice.show()


    

In [ ]:

    

#finfo = ds.field_info[('io', 'particle_ones')]
#finfo.take_log
fn = ds.add_deposited_particle_field(('jetp', 'particle_gamc'), 'nearest') 
print fn


    

In [ ]:

    

#reg = ds.region(center=[ 0.,  0.,  0.], left_edge=[ -5.0e21,  -1.2e+22,  -2.4e+22], right_edge=[  5.0e21,   1.2e+22,   2.4e+22])
#print reg['all', 'particle_position_x']
#print ad.left_edge
#print ad.right_edge


slice = yt.SlicePlot(ds, 'x', fields=fn)#, center=(0.0,0.0,6.17E22), width=(40,'kpc'))
slice.zoom(10)
slice.set_log(fn, True)
slice.set_zlim(fn, 1E2, 1E5)
slice.annotate_grids()
slice.set_origin('native')
slice.show()


    

In [ ]:

    

ad = ds.all_data()
grids = ds.index.grids
print len(grids)
null = plt.hist(ds.index.grid_particle_count, bins=100, range=(0,500), log=True)
#print grids


    

In [ ]:

    

ad = ds.all_data()

shokfil = np.logical_and((ad['io', 'particle_shok']==1), (ad['io', 'particle_gamc']>0))
jetfil = np.logical_and((ad['io', 'particle_shok']==0), (ad['io', 'particle_gamc']>0))

shokpart = (ad['io', 'particle_gamc'][shokfil])
jetpart = (ad['io', 'particle_gamc'][jetfil])
n,bins,patches = plt.hist(np.log10(shokpart), bins=500, alpha=0.5, range=(0,50))
n,bins,patches = plt.hist(np.log10(jetpart), bins=500, alpha=0.5, range=(0,50))
#plt.semilogy()


    

In [ ]:

    

p = yt.ParticlePlot(ds, ('xcenter', 'particle_posy'), ('xcenter', 'particle_posz'), ('xcenter', 'particle_gamc'))
p.set_zlim(('xcenter', 'particle_gamc'), 1E2, 1E6)
p.show()


    

In [ ]:

    

#obj = ds.covering_grid(level=0, left_edge=ds.domain_left_edge,\
#                                      dims=ds.domain_dimensions)
obj = ds.region(
len(obj['io','particle_dens'])
#slice = yt.SlicePlot(obj, 'x', fields='pressure')
#print obj.fcoords.shape


    

In [ ]:

    




    

In [ ]:

    

ds.add_deposited_particle_field(('jetp', 'particle_gamc'), 'cic')
dim=128
zoom=1
obj = ds.arbitrary_grid(ds.domain_left_edge/zoom, ds.domain_right_edge/zoom, dims=[dim,dim,dim*2])
#print obj['jetp_cic_gamc'][dim/2,:,:].shape
#print obj.fcoords.shape
fig = plt.figure(figsize=(10,15),dpi=200)
ext = np.array([-9.65E22, 9.65E22, -1.93E23, 1.93E23])/3.08567758E21
ims = plt.imshow(np.log10(obj['deposit', 'jetp_cic_gamc'][dim/2,:,:]).transpose(), extent=ext, \
                 origin='lower', cmap='algae', vmin=2, vmax=5, interpolation='nearest')

cb=plt.colorbar(ims)
cb.set_label(u'log $\gamma_c$')
plt.xlabel('y (kpc)')
plt.ylabel('z (kpc)')


    

In [ ]:

    

def _gamc_interp(field, data):
    dim=512
    obj = data.ds.arbitrary_grid([-9.65E22, -9.65E22, -1.93E23], [9.65E22, 9.65E22, 1.93E23], dims=[dim, dim, dim*2])
    x0 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,0].min()
    x1 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,0].max()
    y0 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,1].min()
    y1 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,1].max()
    z0 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,2].min()
    z1 = obj.fcoords.reshape([dim,dim,dim*2,3])[:,:,:,2].max()
    boundaries = [x0,x1,y0,y1,z0,z1]
    interp = yt.utilities.linear_interpolators.TrilinearFieldInterpolator(obj['deposit', 'jetp_cic_gamc'], \
                                                                 boundaries, ['x','y','z'], True)
    field_data = interp(data)

    return field_data

ds.add_field(('deposit', 'gamc_interp'), function=_gamc_interp)


    

In [ ]:

    

#ad=ds.all_data()
#gamc_interp = ad['gamc_interp']
#print gamc_interp.shape
print gamc_interp.min()


    

In [ ]:

    

slice = yt.SlicePlot(ds, 'x', fields='gamc_interp')
slice.set_zlim('gamc_interp', 1E2, 1E5)
slice.zoom(32)
slice.annotate_grids()
slice.set_origin('native')
slice.show()


    

In [ ]:

    

ptype = 'xcenter'
ad = ds.all_data()
tag = ad[ptype, 'particle_tag']
#filter = np.logical_and((tag % 20 == 0), (ad['all', 'particle_shok']==0))
filter = np.logical_and.reduce(((tag % 1 == 0), (abs(ad[ptype, 'particle_posx']) < 3.08567758E21)))
y = ad[ptype, 'particle_position_y'][filter]/3.08567758E21
z = ad[ptype, 'particle_position_z'][filter]/3.08567758E21
gamc = np.log10(np.abs(ad[ptype, 'particle_gamc'][filter]))
#plt.hist(tag, bins=140)
#plt.scatter(ad['all', 'particle_tadd'], ad['all', 'particle_tag'], s=1,linewidth=0)
fig=plt.figure(figsize=(6.2,10),dpi=200)
ax=fig.add_subplot(111)
ext = np.array([-30, 30, -60, 60])
ax.set_xlim(ext[0],ext[1])
ax.set_ylim(ext[2],ext[3])
ax.annotate('%6.3f Myr' % (float(ds.current_time)/3.15569E13),\
            (1,0), xytext=(0.05, 0.96),  textcoords='axes fraction',\
            horizontalalignment='left', verticalalignment='center')
sc=ax.scatter(y,z,s=1,c=gamc,linewidth=0,cmap='algae',vmin=2,vmax=5,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label(u'log $\gamma_c$')
plt.xlabel('y (kpc)')
plt.ylabel('z (kpc)')
plt.show()


    

In [ ]:

    

ad = ds.all_data()
pcount = ad['deposit', 'jetp_count']


    

In [ ]:

    

n,bins,patches = plt.hist(pcount, bins=13, range=(-0.5,12.5), log=True)
#print n, bins
bin_centers = (bins[1:] + bins[:-1])/2.0
print sum(n[:]*bin_centers)
print sum(n[1:])
print sum(n[0])/sum(n[:]), sum(n[1:])/sum(n[:])
plt.xlim(-0.5,12.5)
plt.xlabel('# of particles per cell')
plt.ylabel('count of cells')
#plt.text(7, 2E6, 'jet component >= 0.01')


    

In [ ]:

    

from matplotlib.colors import LogNorm
counts, xedges, yedges, Image = \
    plt.hist2d(np.log10(ad['jet ']), allcount, bins=(200, 13), range=[[-11, 0.0], [-0.5,12.5]], norm=LogNorm())
plt.colorbar()
plt.xlabel('log jet')
plt.ylabel('# particles per cell')


    

In [ ]:

    

ds.particle_unions['all'].name
ad = ds.all_data()
ad['particle_position_z']
slab_width = ds.domain_width.value[0]
tag = ad['particle_tag']
tag.argmin()


    

In [ ]:

    

from yt.fields.particle_fields import add_nearest_neighbor_field

fn, = add_nearest_neighbor_field("all", "particle_position", ds, 4)
print fn
dd = ds.all_data()
print dd[fn]


    

In [ ]: