In [ ]:
%matplotlib inline
import numpy as np
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 200
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from mpl_toolkits.axes_grid1 import make_axes_locatable
import yt
yt.mylog.setLevel("WARNING")
from scipy import stats
from particles.particle_filters import *
In [ ]:
fname = '/d/d12/ychen/2016_production_runs/1212_L45_M10_b1_h0_10Myr/data/MHD_Jet_10Myr_hdf5_part_1364'
ds = yt.load(fname)
print(ds.current_time.in_units('Myr'))
ds.add_particle_filter("metal")
#ds.add_particle_filter('jetp')
In [ ]:
f0 = '/d/d12/ychen/2016_production_runs/1212_L45_M10_b1_h0_10Myr/data/MHD_Jet_hdf5_part_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 [ ]:
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(pad=0)
plt.xlim(-100,100)
plt.ylim(-100,100)
plt.xlabel('y (kpc)')
plt.ylabel('z (kpc)')
#plt.colormaps('initial radius (kpc)')
In [ ]:
ls = {(0,10): ['solid', 2],
(10,20): ['dotted', 2],
(20,40): ['dashed', 1],
(40,65): ['solid', 1],
(65,100): ['dotted', 1]}
rmin, rmax = -10, 100
ds.add_particle_filter("metal")
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])
rr = ad[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
drr = rr-rr0
fig = plt.figure(figsize=(4.5,3.5))
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
for rlim in ls.keys():
mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
#ax.scatter(drr[mask], rr0[mask], s=1, lw=0)
ax.plot(drr[mask], rr0[mask], '.', markersize=2, markeredgewidth=0, alpha=0.7, rasterized=True)
null = ax2.hist(drr[mask], bins=rmax-rmin, range=(rmin, rmax),
normed=True, histtype='step',
linestyle=ls[rlim][0], linewidth=ls[rlim][1],
cumulative=True,
label=r'%2i < $r_0$ < %2i' % rlim)
ax.text(rmax-20, rmax-15, '%.0f Myr' % ds.current_time.in_units('Myr'))
ax.set_ylabel(r'initial radius $R_0$ [kpc]')
ax.set_ylim(0, rmax)
ax.axvline(x=0, lw=1, ls=':', color='k')
ax.tick_params(labelbottom=False, direction='in')
ax.set_yticks([20, 40, 60, 80, 100])
ax.set_yticklabels([20, 40, 60, 80, 100])
#ax.set_aspect(1)
handles, labels = ax2.get_legend_handles_labels()
new_handles = [Line2D([], [], c=h.get_edgecolor(),
ls=h.get_linestyle(),
lw=h.get_linewidth())
for h in handles]
ax2.legend(new_handles[::-1], labels[::-1],
ncol=2, fontsize=8, columnspacing=0.5,
loc=4)
ax2.axvline(x=0, lw=1, ls=':', color='k')
ax2.grid(ls='--', alpha=0.5)
ax2.set_xlabel(r'radial displacement $\Delta R$ [kpc]')
ax2.set_xlim(rmin, rmax)
ax2.set_ylabel('cumulative fraction')
ax2.set_ylim(0.6, 1.0)
plt.tight_layout()
fig.savefig('particles_dr_%s.pdf' % ds.basename[-4:], bbox_inches='tight')
In [ ]:
ls = {(0,10): ['solid', 1.5],
(10,20): ['dotted', 1.5],
(20,40): ['dashed', 1],
(40,65): ['solid', 1],
(65,100): ['dotted', 1]}
rmin, rmax = -10, 100
ds.add_particle_filter("metal")
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])
rr = ad[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
drr = rr-rr0
fig = plt.figure(figsize=(4.5,5.5))
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
ax3 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
for rlim in ls.keys():
mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
#ax.scatter(drr[mask], rr0[mask], s=1, lw=0)
ax.plot(drr[mask], rr0[mask], '.', markersize=2, markeredgewidth=0, alpha=0.7, rasterized=True)
hist, bin_edges = np.histogram(drr[mask], bins=int((rmax-rmin)*4), range=(rmin, rmax), density=True)
bin_center = (bin_edges[:-1]+bin_edges[1:])/2
ax2.plot(bin_center, np.cumsum(hist)*(bin_edges[1:]-bin_edges[:-1]), linestyle=ls[rlim][0], linewidth=ls[rlim][1],
label=r'%2i < $r_0$ < %2i' % rlim)
null = ax3.hist(drr[mask], bins=int((rmax-rmin)), range=(rmin, rmax),
normed=True, histtype='step', bottom=1E-6,
linestyle=ls[rlim][0], linewidth=ls[rlim][1],
label=r'%2i < $r_0$ < %2i' % rlim)
ax.text(rmax-20, rmax-15, '%.0f Myr' % ds.current_time.in_units('Myr'))
ax.set_ylabel(r'initial radius $R_0$ [kpc]')
ax.set_ylim(0, rmax)
ax.axvline(x=0, lw=0.5, ls='-', color='k')
ax.tick_params(labelbottom=False, direction='in')
ax.set_yticks([20, 40, 60, 80, 100])
ax.set_yticklabels([20, 40, 60, 80, 100])
#ax.set_aspect(1)
ax2.tick_params(labelbottom=False, direction='in')
#handles, labels = ax2.get_legend_handles_labels()
#new_handles = [Line2D([], [], c=h.get_edgecolor(),
# ls=h.get_linestyle(),
# lw=h.get_linewidth())
# for h in handles]
#ax3.legend(new_handles[::-1], labels[::-1],
ax2.legend(ncol=2, fontsize=8, columnspacing=0.5, loc=4)
ax2.axvline(x=0, lw=0.5, ls='-', color='k')
ax2.grid(ls='--', alpha=0.5)
ax2.set_ylabel('cumulative fraction')
ax2.set_ylim(0.0, 1.0)
#ax2.set_yticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax3.axvline(x=0, lw=0.5, ls='-', color='k')
ax3.grid(ls='--', alpha=0.5)
ax3.semilogy()
ax3.set_ylabel('fraction')
ax3.set_ylim(1E-4, 0.8)
ax3.set_xlabel(r'radial displacement $\Delta R$ [kpc]')
ax3.set_xlim(rmin, rmax)
plt.tight_layout()
fig.savefig('particles_dr_histogram_%s.pdf' % ds.basename[-4:], bbox_inches='tight')
In [ ]:
ls = {(0,10): ['solid', 1.5],
(10,20): ['dotted', 1.5],
(20,40): ['dashed', 1],
(40,65): ['solid', 1],
(65,100): ['dotted', 1]}
rmin, rmax = -10, 100
bw = 1
ds.add_particle_filter("metal")
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])
rr = ad[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
drr = rr-rr0
fig = plt.figure(figsize=(4.5,5.5))
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
ax3 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
for rlim in ls.keys():
mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
#ax.scatter(drr[mask], rr0[mask], s=1, lw=0)
ax.plot(drr[mask], rr0[mask], '.', markersize=2, markeredgewidth=0, alpha=0.7, rasterized=True)
std = np.std(drr[mask].v)
#kde = stats.gaussian_kde(drr[mask].v, bw_method='scott')
#print(rlim, std, kde.factor, std*kde.factor)
kde = stats.gaussian_kde(drr[mask].v, bw_method=bw/std)
print(rlim, std, kde.factor, std*kde.factor)
rrlin = np.linspace(rmin,rmax,num=220)
cum_frac = [kde.integrate_box_1d(-20, r) for r in rrlin]
ax2.plot(rrlin, cum_frac, linestyle=ls[rlim][0], linewidth=ls[rlim][1], label=r'%2i < $R_0$ < %2i' % rlim)
prob = kde.pdf(rrlin)
line, = ax3.plot(rrlin, prob, linestyle=ls[rlim][0], linewidth=ls[rlim][1])
ax.text(rmax-20, rmax-10, '%.0f Myr' % ds.current_time.in_units('Myr'))
ax.set_ylabel(r'initial radius $R_0$ [kpc]')
ax.set_ylim(0, rmax)
ax.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax.tick_params(labelbottom=False, direction='in')
ax.set_yticks([20, 40, 60, 80, 100])
ax.set_yticklabels([20, 40, 60, 80, 100])
#ax.set_aspect(1)
handles, labels = ax2.get_legend_handles_labels()
new_handles = [Line2D([], [], c=h.get_color(),
ls=h.get_linestyle(),
lw=h.get_linewidth())
for h in handles]
ax2.legend(new_handles[::-1], labels[::-1],
ncol=2, fontsize=8, columnspacing=0.5,
loc=4)
ax2.tick_params(labelbottom=False, direction='in')
ax2.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax2.grid(ls='--', alpha=0.5)
ax2.set_ylabel('cumulative fraction')
ax2.set_ylim(0.0, 1.0)
ax3.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax3.grid(ls='--', alpha=0.5)
ax3.semilogy()
ax3.set_ylabel('probability density')
ax3.set_ylim(1E-4, 0.7)
ax2.set_yticks([0,0.25,0.5,0.75,1])
ax2.set_yticklabels(['0','1/4','1/2','3/4','1'])
ax3.errorbar(53, 0.2, xerr=2.35482*bw/2, elinewidth=1, capsize=1, color='k')
ax3.annotate('KDE Gaussian FWHM', (53, 0.2), (58, 0.18),
horizontalalignment='left',
verticalalignment='center')
ax3.set_xlabel(r'radial displacement $\Delta R$ [kpc]')
ax3.set_xlim(rmin, rmax)
plt.tight_layout()
fig.savefig('particles_dr_kde_bw%i_%s.pdf' % (bw, ds.basename[-4:]), bbox_inches='tight')
In [ ]:
ls = {(0,10): ['solid', 1.5],
(10,20): ['dotted', 1.5],
(20,40): ['dashed', 1],
(40,65): ['solid', 1],
(65,100): ['dotted', 1]}
rmin, rmax = -10, 100
bw = 4
ds.add_particle_filter("metal")
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])
rr = ad[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
drr = rr-rr0
fig = plt.figure(figsize=(4.5,5.5))
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
ax3 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
drr = np.array([35, 65, 65, 65, 85])
mask = [True, True, True, True, True]
ax.plot(drr[mask], rr0[mask], '.', markersize=2, markeredgewidth=0, alpha=0.7, rasterized=True)
std = np.std(drr[mask])
kde = stats.gaussian_kde(drr[mask], bw_method=bw/std)
print(rlim, std, kde.factor)
rrlin = np.linspace(rmin,rmax,num=110)
cum_frac = [kde.integrate_box_1d(-20, r) for r in rrlin]
ax2.plot(rrlin, cum_frac, linestyle=ls[rlim][0], linewidth=ls[rlim][1], label=r'%2i < $r_0$ < %2i' % rlim)
prob = kde.pdf(rrlin)
ax3.plot(rrlin, prob, linestyle=ls[rlim][0], linewidth=ls[rlim][1])
ax.text(rmax-20, rmax-15, '%.0f Myr' % ds.current_time.in_units('Myr'))
ax.set_ylabel(r'initial radius $R_0$ [kpc]')
ax.set_ylim(0, rmax)
ax.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax.tick_params(labelbottom=False, direction='in')
ax.set_yticks([20, 40, 60, 80, 100])
ax.set_yticklabels([20, 40, 60, 80, 100])
#ax.set_aspect(1)
handles, labels = ax2.get_legend_handles_labels()
new_handles = [Line2D([], [], c=h.get_color(),
ls=h.get_linestyle(),
lw=h.get_linewidth())
for h in handles]
ax2.legend(new_handles[::-1], labels[::-1],
ncol=2, fontsize=8, columnspacing=0.5,
loc=4)
ax2.tick_params(labelbottom=False, direction='in')
ax2.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax2.grid(ls='--', alpha=0.5)
ax2.set_ylabel('estimated CDF')
ax2.set_ylim(0.0, 1.0)
ax3.axvline(x=0, lw=0.5, ls='-', color='k', alpha=0.7)
ax3.grid(ls='--', alpha=0.5)
#ax3.semilogy()
ax3.set_ylabel('estimated PDF')
ax3.set_ylim(1E-4, 0.06)
ax3.errorbar(65, 0.03, xerr=2.35482*bw/2, elinewidth=1, capsize=1, color='k')
ax3.annotate('kde bandwidth', (65, 0.2), (70, 0.2),
horizontalalignment='left',
verticalalignment='center')
ax3.set_xlabel(r'radial displacement $\Delta R$ [kpc]')
ax3.set_xlim(rmin, rmax)
#plt.tight_layout()
#fig.savefig('particles_dr_kde_bw%i_%s.pdf' % (bw, ds.basename[-4:]), bbox_inches='tight')
In [ ]:
ls = {(0,10): ['solid', 2],
(10,20): ['dotted', 2],
(20,40): ['dashed', 1],
(40,65): ['solid', 1],
(65,120): ['dotted', 1]}
rmin, rmax = 0, 120
bw = 1
ds.add_particle_filter("metal")
ad = ds.all_data()
arr = np.argsort(ad[('metal', 'particle_tag')])
rr = ad[('metal','particle_position_spherical_radius')].in_units('kpc')[arr]
fig = plt.figure(figsize=(4.5,5.5))
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
ax2 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
ax3 = divider.append_axes("bottom", 1.0, pad=0.0, sharex=ax)
for rlim in ls.keys():
mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
#ax.scatter(drr[mask], rr0[mask], s=1, lw=0)
ax.plot(rr[mask], rr0[mask], '.', markersize=2, markeredgewidth=0, alpha=0.7, rasterized=True)
std = np.std(drr[mask].v)
kde = stats.gaussian_kde(rr[mask].v, bw_method=bw/std)
print(rlim, np.std(rr[mask]), kde.factor)
rrlin = np.linspace(rmin,rmax,num=110)
cum_frac = [kde.integrate_box_1d(-20, r) for r in rrlin]
ax2.plot(rrlin, cum_frac, linestyle=ls[rlim][0], linewidth=ls[rlim][1], label=r'%2i < $r_0$ < %2i' % rlim)
prob = kde.pdf(rrlin)
ax3.plot(rrlin, prob, linestyle=ls[rlim][0], linewidth=ls[rlim][1])
ax.text(rmax-20, rmax-15, '%.0f Myr' % ds.current_time.in_units('Myr'))
ax.plot(rrlin, rrlin, c='k', alpha=0.7)
ax.set_ylabel(r'initial radius $R_0$ [kpc]')
ax.set_ylim(0, rmax)
ax.axvline(x=0, lw=1, ls=':', color='k')
ax.tick_params(labelbottom=False, direction='in')
ax.set_yticks([20, 40, 60, 80, 100])
ax.set_yticklabels([20, 40, 60, 80, 100])
#ax.set_aspect(1)
handles, labels = ax2.get_legend_handles_labels()
new_handles = [Line2D([], [], c=h.get_color(),
ls=h.get_linestyle(),
lw=h.get_linewidth())
for h in handles]
ax2.legend(new_handles[::-1], labels[::-1],
ncol=2, fontsize=8, columnspacing=0.5,
loc=4)
ax2.tick_params(labelbottom=False, direction='in')
ax2.axvline(x=0, lw=1, ls=':', color='k')
ax2.grid(ls='--', alpha=0.5)
ax2.set_ylabel('cumulative fraction')
ax2.set_ylim(0.0, 1.0)
ax3.axvline(x=0, lw=1, ls=':', color='k')
ax3.grid(ls='--', alpha=0.5)
ax3.semilogy()
ax3.set_ylabel('fraction')
ax3.set_ylim(1E-4, 0.2)
ax3.set_xlabel(r'current radius $R$ [kpc]')
ax3.set_xlim(rmin, rmax)
plt.tight_layout()
fig.savefig('particles_r_kde_%s.pdf' % ds.basename[-4:], bbox_inches='tight')
In [ ]:
from sklearn.neighbors import KernelDensity
for rlim in ls.keys():
mask = np.logical_and(rr0 > rlim[0], rr0 < rlim[1])
# instantiate and fit the KDE model
kde = KernelDensity(bandwidth=2.0, kernel='gaussian')
kde.fit(drr[mask,None])
# score_samples returns the log of the probability density
rrlin = np.linspace(rmin,rmax,num=200)
logprob = kde.score_samples(rrlin[:,None])
plt.plot(rrlin, np.exp(logprob))
plt.semilogy()
plt.ylim(1E-4, 0.5)
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')