In [ ]:

    

import yt
from synchrotron.yt_synchrotron_emissivity import StokesFieldName, setup_part_file
from particle_filters import *
ptype = 'lobe'
nu = (100, 'MHz')
ds = yt.load('/d/d5/ychen/2015_production_runs/0529_h1_00/data/MHD_Jet_hdf5_plt_cnt_0600_synchrotron_peak_gamc_gc32')
print(setup_part_file(ds))
ds.add_particle_filter(ptype)
stokes = StokesFieldName(ptype, nu, 'x', field_type='flash')
#res = ds.domain_dimensions[1:]*ds.refine_by**ds.index.max_level//zoom_fac//2
kpc = yt.units.kpc


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', stokes.I, center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])
plot.set_log(stokes.I, True)
plot.set_zlim(stokes.I, 1E-28, 1E-20)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
plot.annotate_grids(linewidth=0.5, alpha=0.5)#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

plot.annotate_clear()
plot.show()


    

In [ ]:

    

import yt
from synchrotron.yt_synchrotron_emissivity import StokesFieldName, setup_part_file
from particle_filters import *
ptype = 'lobe'
nu = (100, 'MHz')
#ds = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0910_synchrotron_peak_gamc_gc32')
ds = yt.load('/d/d5/ychen/2015_production_runs/0529_h1_00/data/MHD_Jet_hdf5_plt_cnt_0600_synchrotron_peak_gamc_gc32')
print(setup_part_file(ds))
ds.add_particle_filter(ptype)
stokes = StokesFieldName(ptype, nu, 'x', field_type='flash')
#res = ds.domain_dimensions[1:]*ds.refine_by**ds.index.max_level//zoom_fac//2
kpc = yt.units.kpc
plot = yt.SlicePlot(ds, 'x', 'lobe_nnw_gamc_dtau', center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])


    

In [ ]:

    

plot.set_zlim('lobe_nnw_gamc_dtau', 1E3, 1E6)
#plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
plot.show()


    

In [ ]:

    

plot.set_zlim(stokes.I, 1E-28, 1E-22)


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', stokes.I, center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (60, 'kpc')])
plot.set_log(stokes.I, True)
plot.set_zlim(stokes.I, 1E-28, 1E-6)
#plot.set_buff_size((200,200))1
plot.set_axes_unit('kpc')
plot.annotate_grids(linewidth=0.5, alpha=0.5)#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', stokes.I, center=ds.arr([0.01,0.0,0.0], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])
plot.set_log(stokes.I, True)
plot.set_zlim(stokes.I, 1E-28, 1E-20)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
plot.annotate_grids(linewidth=0.5, alpha=0.5)#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', stokes.I, center=ds.arr([0.1,0.1,40], input_units='kpc'), width=[(10, 'kpc'), (10, 'kpc')],)
plot.set_log(stokes.I, True)
plot.set_zlim(stokes.I, 1E-30, 1E-20)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
plot.annotate_grids()#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', stokes.I, center=ds.arr([0.1,0.1,40], input_units='kpc'), width=[(10, 'kpc'), (10, 'kpc')],)
plot.set_log(stokes.I, True)
plot.set_zlim(stokes.I, 1E-30, 1E-20)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
plot.annotate_grids()#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

import yt
from synchrotron.yt_synchrotron_emissivity import setup_part_file
ds = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0910')
setup_part_file(ds)
plot = yt.SlicePlot(ds, 'x', 'magnetic_field_strength', center=ds.arr([0.01,0.01,0.01], input_units='kpc'), width=[(60, 'kpc'), (120, 'kpc')])
plot.show()


    

In [ ]:

    

plot.annotate_clear()
plot.set_zlim('magnetic_field_strength', 1E-7, 2E-5)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
#plot.annotate_grids(linewidth=0.5, alpha=0.2)#, draw_ids=True)
#plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')


    

In [ ]:

    

from synchrotron.yt_synchrotron_emissivity import _jet_volume_fraction
ds = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0910')
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-5, 1E0)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
#plot.annotate_grids(linewidth=0.5, alpha=0.2)#, draw_ids=True)
#plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')


    

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-5, 1E0)
plot.show()


    

In [ ]:

    

import matplotlib.pyplot as plt
ad = ds.all_data()
plt.hist(np.log10(ad['jet_volume_fraction']), bins='auto', range=(-5, 0))


    

In [ ]:

    

import matplotlib.pyplot as plt
ad = ds.all_data()
ptype = 'lobe'
print(ds.add_particle_filter(ptype))
tag = ad[ptype, 'particle_tag']
filter = np.abs(ad[ptype, 'particle_position_x']) < 0.25*yt.units.kpc.in_units('cm').v
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=(10,10))
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=5,c=c,linewidth=0,cmap='arbre', vmin=2,vmax=3.5,alpha=0.8)
cb=plt.colorbar(sc)
cb.set_label('gamc')
plt.show()


    

In [ ]:

    

print(ad['all', 'particle_tag'].shape)
print(ad['jet', 'particle_tag'].shape)
print(ad['jetp', 'particle_tag'].shape)
print(ad['lobe', 'particle_tag'].shape)


    

In [ ]:

    

ds_flash = yt.load('/home/ychen/data/0only_0529_h1/MHD_Jet_hdf5_plt_cnt_0500', 
             particle_filename='/home/ychen/data/0only_0529_h1/MHD_Jet_hdf5_part_0500_updated' )

kpc = yt.units.kpc
plot = yt.SlicePlot(ds_flash, 'x', 'velocity_magnitude', center=ds.arr([0,0,0], input_units='kpc'),
                    width=[(30, 'kpc'), (60, 'kpc')])
plot.set_axes_unit('kpc')
plot.annotate_grids(linewidth=0.5, alpha=0.5)#, draw_ids=True)
#plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

ad = ds_flash.all_data()
lobe_ad = ad.cut_region(["obj['velocity_magnitude'] < 3.0e8"])
plot = yt.SlicePlot(ds_flash, 'x', 'velocity_magnitude', center=ds.arr([0,0,0], input_units='kpc'),
                    width=[(30, 'kpc'), (60, 'kpc')], data_source=lobe_ad)

plot.show()


    

In [ ]:

    

ad = ds_flash.all_data()
lobe1_ad = ad.cut_region(["obj['velocity_magnitude'] < 1.5e9"])
plot = yt.SlicePlot(ds_flash, 'x', 'velocity_magnitude', center=ds.arr([0,0,0], input_units='kpc'),
                    width=[(30, 'kpc'), (60, 'kpc')], data_source=lobe1_ad)
plot.show()


    

In [ ]:

    

from synchrotron.yt_synchrotron_emissivity import _jet_volume_fraction
ds_flash.add_field(('gas', 'jet_volume_fraction'), function=_jet_volume_fraction,
              display_name="Jet Volume Fraction", sampling_type='cell')
plot = yt.SlicePlot(ds_flash, 'x', 'jet ', center=ds.arr([0,0,0], input_units='kpc'),
                    width=[(30, 'kpc'), (60, 'kpc')])
plot.set_zlim('jet ', 1E-2, 1)
plot.show()


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', vstokes.I, center=ds.arr([0,0,0], input_units='kpc'), width=(2, 'kpc'))
#plot.set_log(stokes.I, True)
#plot.set_zlim(stokes.I, 1E-28, 1E-20)
plot.set_axes_unit('kpc')
#plot.set_buff_size((200,200))
plot.annotate_grids(linewidth=0.5, alpha=0.5, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

ds.field_list


    

In [ ]:

    

g = ds.index.grids[119376]
g['lobe', 'particle_gamc']


    

In [ ]:

    

plot_zoom = yt.SlicePlot(ds, 'y', stokes.I, center=ds.arr([0,1-0.15,2.5-0.15], input_units='kpc'), width=(2, 'kpc'))
plot_zoom.set_log(stokes.I, True)
plot_zoom.set_zlim(stokes.I, 1E-28, 1E-20)
plot.set_buff_size((200,200))
plot_zoom.set_axes_unit('kpc')
plot_zoom.annotate_clear()
plot_zoom.annotate_grids(linewidth=0.5, alpha=0.5, draw_ids=True)
plot_zoom.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype=ptype)
#plot.annotate_cell_edges(alpha=0.1)
plot_zoom.show()


    

In [ ]:

    

extend_cells=32
g = ds.index.grids[5852]
print(g)
print(g.LeftEdge.in_units('kpc'))
print(g.RightEdge.in_units('kpc'))
print(g.dds.in_units('kpc'))
stokes = StokesFieldName(ptype, (150, 'MHz'), 'x', field_type='flash')
I = g[stokes.I]
ad = ds.all_data()
for j, pos in enumerate(ad[ptype, 'particle_position']):
    continue_loop = False
    for i in [0,1,2]:
        if pos[i] < g.LeftEdge[i] - extend_cells*g.dds[i] or pos[i] > g.RightEdge[i] + extend_cells*g.dds[i]:
            continue_loop = True
    if continue_loop:
        continue
    else:
        print(j, pos.in_units('kpc'))
print(g[ptype, 'particle_den0'])


    

In [ ]:

    

print(g[stokes.I])


    

In [ ]:

    

import yt
from synchrotron.yt_synchrotron_emissivity import add_synchrotron_dtau_emissivity, StokesFieldName
ptype = 'lobe'
nu = (100, 'MHz')
ds_flash = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0910',
             particle_filename='/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_part_0910_updated') 
stokesd = StokesFieldName(ptype, nu, 'x')
add_synchrotron_dtau_emissivity(ds_flash, ptype=ptype, nu=nu, proj_axis='x', extend_cells=8)


    

In [ ]:

    

plot = yt.SlicePlot(ds_flash, 'x', 'lobe_nnw_gamc-').zoom(8)
plot.show()


    

In [ ]:

    

g = ds_flash.index.grids[21004]
print(g.id)
left_edge = g.LeftEdge - 64*g.dds
right_edge = g.RightEdge + 64*g.dds
box = ds_flash.box(left_edge, right_edge)

box['lobe', 'particle_sync_spec_100.0MHz']
#g['deposit', 'lobe_nn_sync_spec_150.0MHz']


    

In [ ]:

    

g['deposit', 'lobe_nnw_sync_spec_100.0MHz']
#g['deposit', 'nn_emissivity_i_lobe_100.0MHz_x']


    

In [ ]:

    

g = ds.index.grids[21004]
print(g.id)
g['flash', 'nn_emissivity_i_lobe_100.0MHz_x']


    

In [ ]:

    

np.exp(np.log(np.float32(np.finfo(np.float64).tiny)))


    

In [ ]:

    

plot = yt.SlicePlot(ds_flash, 'z', stokesd.I, center=ds_flash.arr([0,1-0.15,2.5-0.15], input_units='kpc'), width=(2, 'kpc'))
plot.set_log(stokesd.I, True)
plot.set_zlim(stokesd.I, 1E-28, 1E-20)
#plot.set_buff_size((200,200))
plot.set_axes_unit('kpc')
plot.annotate_grids(linewidth=0.5, alpha=0.5)#, draw_ids=True)
plot.annotate_particles(p_size=4, width=(0.25, 'kpc'), ptype='lobe')
#plot.annotate_cell_edges(alpha=0.1)
plot.show()


    

In [ ]:

    

fields = ['pressure', 'magnetic_pressure']
p = yt.SlicePlot(ds_flash, 'z', fields, center=ds.arr([0,1-0.15,2.5-0.15], input_units='kpc'), width=(2, 'kpc'))
p.show()


    

In [ ]:

    

ad = ds_flash.all_data()
#print(ad[ptype, 'particle_gamc'])
deposit_field = 'particle_sync_spec_%s' % stokesd.nu_str
#print(ad[ptype, deposit_field])
field_name = "%s_nn_%s" % (ptype, deposit_field.replace('particle_', ''))
#print(ad['deposit', field_name])


    

In [ ]:

    

p = yt.SlicePlot(ds_flash, 'z', field_name, center=ds_flash.arr([0,1-0.15,2.5-0.15], input_units='kpc'), width=(2, 'kpc'))
p.annotate_particles(p_size=10, width=(2, 'kpc'), ptype=ptype)
p.annotate_grids(alpha=0.5, draw_ids=True)
p.show()


    

In [ ]:

    

import yt.geometry.particle_deposit as particle_deposit
cls = particle_deposit.deposit_nearest
g_flash = ds_flash.index.grids[5852]
op = cls(tuple(g_flash.ActiveDimensions[::-1]), 'cubic')


    

In [ ]:

    

import numpy as np
ds_flash.add_particle_filter(ptype)
ad = ds_flash.all_data()
fields = [ad[ptype, deposit_field]]
fields = [np.ascontiguousarray(f) for f in fields]


    

In [ ]:

    

g_flash = ds_flash.index.grids[5852]
print(g_flash)
op.initialize()
op.process_grid(g_flash, ad[ptype, 'particle_position'], fields, extend_cells=32)


    

In [ ]:

    

print(g_flash.LeftEdge)
center = (g_flash.LeftEdge + g_flash.dds*7.5)
print(g_flash.RightEdge)
print(g_flash.dds, np.sum(8*8*g_flash.dds**2))
dist = center.v - np.array([-1.5e+21, 9.5e+20, 4.3e+21])
print(np.sum(dist*dist))


    

In [ ]:

    

import matplotlib.pyplot as plt
vals = op.finalize()
null = plt.hist(np.log10(vals.flatten()), bins=50)


    

In [ ]:

    

g_flash = ds_flash.index.grids[5851]
print(g_flash)
tuple(g.ActiveDimensions[::-1])


    

In [ ]:

    

print(g_flash[ptype, 'particle_gamc'])
deposit_field = 'particle_sync_spec_%s' % stokesd.nu_str
print(g_flash[ptype, deposit_field])
field_name = "%s_nn_%s" % (ptype, deposit_field.replace('particle_', ''))
print(g_flash[stokesd.I].v)


    

In [ ]:

    

data = g
los = np.array([1,0,0])
los = los/np.sqrt(np.sum(los*los))
Bvec = np.array([data[(ptype, 'particle_magx')],\
                 data[(ptype, 'particle_magy')],\
                 data[(ptype, 'particle_magz')]])/np.sqrt(4.0*np.pi)
cross = np.cross(los, Bvec, axisb=0)
Bsina = np.sqrt(np.sum(cross*cross, axis=-1))
Bsina = data.apply_units(Bsina, 'gauss')
e = yt.physical_constants.elementary_charge
me = yt.physical_constants.mass_electron
c = yt.physical_constants.c
den1 = data[(ptype, 'particle_den1')]
dtau = data[(ptype, 'particle_dtau')]
gamc = (np.abs(data[(ptype, 'particle_dens')] / den1))**(1./3.)/dtau
print(gamc)
nuc = 3.0*gamc**2*e*Bsina/(4.0*np.pi*me*c) 
print(nuc.in_units('MHz'))


    

In [ ]:

    

import matplotlib.pyplot as plt
import numpy as np
ad = ds.all_data()
null =plt.hist(np.log10(ad['particle_den0']), bins=40, range=(-25.85, -25.7))
from tools import calcDen0
print(calcDen0(g))


    

In [ ]:

    

g = ds.index.grids[-1]


    

In [ ]:

    

ds.index.num_grids


    

In [ ]: