In [ ]:

    

%matplotlib inline
import numpy as np
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
import yt
yt.mylog.setLevel("WARNING")

from particles.particle_filters import *


    

In [ ]:

    

fname = '/home/ychen/data/0only_1212_h0_10Myr_rerun/data/MHD_Jet_10Myr_hdf5_plt_cnt_1050'
ds = yt.load(fname)
ds.add_particle_filter("metal")
#ds.add_particle_filter('jetp')


    

In [ ]:

    

f0 =  '/home/ychen/data/0only_1212_h0_10Myr_rerun/data/MHD_Jet_10Myr_hdf5_plt_cnt_0000'
ds0 = yt.load(f0)
ds0.add_particle_filter("metal")
ad0 = ds0.all_data()
arr0 = np.argsort(ad0[('metal', 'particle_tag')])
rr0 = ad0[('metal','particle_position_spherical_radius')].in_units('kpc')[arr0]


    

In [ ]:

    

ds0.parameter_filename
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])


    

In [ ]:

    

xx=ad[('metal', 'particle_posx')].in_units('kpc')[arr]
yy=ad[('metal', 'particle_posy')].in_units('kpc')[arr]
zz=ad[('metal', 'particle_posz')].in_units('kpc')[arr]

filtr = np.abs(xx)<30

plt.figure(figsize=(12,10))
plt.scatter(yy[filtr][::3], zz[filtr][::3], c=rr0[filtr][::3], s=2, lw=0, vmax=80, vmin=0, cmap='gnuplot2_r')
plt.colorbar() 
plt.xlim(-100,100)
plt.ylim(-100,100)
plt.xlabel('y (kpc)')
plt.ylabel('z (kpc)')
#plt.colormaps('initial radius (kpc)')


    

In [ ]:

    

rr = ds.all_data()[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
drr = rr-rr0


    

In [ ]:

    

from mpl_toolkits.axes_grid1 import make_axes_locatable

ls = {(0,10): ['solid', 2],
      (10,20): ['dotted', 2],
      (20,30): ['dashed', 1],
      (30,60): ['solid', 1],
      (60,90): ['dotted', 1]}

fig = plt.figure(figsize=(5,4))

ax = fig.add_subplot(111)

divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.1, sharex=ax)

for rlim in [(0,10), (10,20), (20,30), (30,60), (60,90)]:
    mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
    ax.scatter(drr[mask], rr0[mask], s=1, lw=0)
    
    null = ax2.hist(drr[mask], bins=65, range=(-5, 60),
                    normed=True, histtype='step',
                    linestyle=ls[rlim][0], linewidth=ls[rlim][1],
                    cumulative=True,
                    label=r'%2i < $r_0$ < %2i' % rlim, alpha=0.9)
ax.text(45, 80, '%.1f Myr' % ds.current_time.in_units('Myr'))
ax.set_ylabel(r'initial radius $r_0$ (kpc)')
ax.set_ylim(0, 90)
ax.axvline(x=0, lw=1, ls=':', color='k')
ax.tick_params(labelbottom=False, direction='in')

handles, labels = ax2.get_legend_handles_labels()
ax2.legend(handles[::-1], labels[::-1], fontsize=8, loc=4)

ax2.set_xlabel(r'$\Delta r$ (kpc)')
ax2.set_xlim(-5, 60)
ax2.set_ylabel('cumulative fraction')
ax2.set_ylim(0.6, 1.0)
plt.tight_layout()
fig.savefig('particle_dr_100Myr.png', dpi=300)


    

In [ ]:

    

arr = np.array(ad[('metal', 'particle_tag')], dtype=int)-1


    

In [ ]:

    

xx[arr]


    

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')