In [ ]:
#%pdb
%matplotlib inline
import numpy as np
import matplotlib
matplotlib.rcParams['savefig.dpi'] = 100
import matplotlib.pyplot as plt
from yt_emissivity import *
In [ ]:
fname = '/d/d4/ychen/2016_production_runs/1111_L45_M10_b01_h1/MHD_Jet_hdf5_plt_cnt_0200'
ds = yt.load(fname)
#print ds.derived_field_list
In [ ]:
nu = (150, 'MHz')
ptype = 'lobe'
pars = add_emissivity(ds, ptype=ptype, nu=nu)
print pars
norm = yt.YTQuantity(*nu).in_units('GHz').value**0.5
field = ('deposit', ('nn_emissivity_%s_%%.1f%%s' % ptype) % nu)
proj_axis = 'x'
plot = yt.SlicePlot(ds, proj_axis, field, center=(0,0,0))
#plot.set_cmap(field, 'Blue-Red')
plot.set_zlim(field, 1E-25/norm, 1E-21/norm)
plot.zoom(20)
plot.show()
In [ ]:
plot = yt.ProjectionPlot(ds, proj_axis, field, center=(0,0,0), max_level=2)
#plot.set_zlim(field, 1E-3/norm, 1E1/norm)
#plot.zoom(20)
plot.show()
In [ ]:
projs = {}
proj_axis = 'x'
ptype = 'jnsp'
for nu in [(150, 'MHz'), (1.4, 'GHz')]:
pars = add_emissivity(ds, ptype=ptype, nu=nu)
field = ('deposit', ('nn_emissivity_%s_%%.1f%%s' % ptype) % nu)
projs[nu] = ds.proj(field, proj_axis, center=[0,0,0])
In [ ]:
ext = ds.arr([-7.72E22, 7.72E22, -1.544E23, 1.544E23], input_units='code_length')
frb1 = projs[(1.4, 'GHz')].to_frb(ext[1]-ext[0], (1600,800), height=(ext[3]-ext[2]))
frb2 = projs[(150, 'MHz')].to_frb(ext[1]-ext[0], (1600,800), height=(ext[3]-ext[2]))
S1 = frb1[('deposit', 'nn_emissivity_%s_1.4GHz' % ptype)]
S2 = frb2[('deposit', 'nn_emissivity_%s_150.0MHz' % ptype)]
alpha = np.log(S1/S2)/np.log(1400/150)
alpha = np.ma.masked_where(S1 < 1E-3, np.array(alpha))
In [ ]:
x=450
y=800
print S1[y,x]
print alpha[y,x]
print alpha.mask[y,x]
print np.array(alpha)[y,x]
In [ ]:
fig = plt.figure(figsize=(8,12), dpi=150)
plt.imshow(alpha, vmin=-2, vmax=-0.5, extent=ext.in_units('kpc'), origin='lower', interpolation='nearest')
cb = plt.colorbar()
In [ ]:
fig = plt.figure(figsize=(8,12), dpi=150)
ims = plt.imshow(np.array(alpha), vmin=-2, vmax=0, extent=ext.in_units('kpc'), origin='lower', interpolation='nearest', aspect='auto')
plt.xlabel('y (kpc)')
plt.ylabel('z (kpc)')
plt.clim(-2.0,-0.5)
cb = plt.colorbar()
cb.set_label('Spectral Index (1.4GHz/150MHz)')
plt.savefig('alpha.png')
In [ ]:
from yt.visualization.plot_window import PlotWindow
plot = PlotWindow(ds.all_data(), ext)
In [ ]:
dir = '/home/ychen/data/0only_0529_h1/'
dfn = 'MHD_Jet_hdf5_plt_cnt_0040'
import shelve
ds = yt.load(dir+dfn)
saved_fn = shelve.open(dir+'projs/'+dfn+'_projs.cpkl')
saved_fn.items
In [ ]:
field = ('gas', 'density')
L = [1,0,1]
north_vector = [0,0,1]
proj = yt.OffAxisProjectionPlot(ds, L, field, center=(0,0,0), north_vector=north_vector)
In [ ]:
proj.set_zlim(field, 2.2E-2, 3.1E-2)
proj.set_log(field, False)
proj.show()
In [ ]:
def _deposit_average_filling(field, grid):
#data._debug()
ptype = 'jetp'
deposit_field = 'particle_emissivity'
uq = grid['gas', 'pressure'].uq
if isinstance(grid, FieldDetector): return grid[ptype, deposit_field]*uq
if len(grid[ptype, deposit_field]) > 0: return grid['gas', 'pressure']*grid[ptype, deposit_field].mean()
#else: return data['zeros']*uq
elif grid.Level == 0: return grid['zeros']*uq
else:
pfield = np.zeros(0)
grids_to_go = grid.Parent.Children
#for i, gr in enumerate(grids_to_go):
# if len(gr.Children) > 0:
# del grids_to_go[i]
# grids_to_go.extend(gr.Children)
for gr in grids_to_go:
pfield = np.concatenate([pfield, gr[ptype, deposit_field]])
if len(pfield) == 0:
return grid['zeros']*uq
else:
return grid['gas', 'pressure']*pfield.mean()
ds.add_field(('deposit', 'avgfill_emissivity'), function=_deposit_average_filling, validators=[ValidateGridType()],
display_name="%s Emissivity" % nu, units='erg/s/cm**3/Hz', take_log=True, force_override=True)
In [ ]:
logging.getLogger('yt').setLevel(logging.INFO)
slice = yt.SlicePlot(ds, 'x', ('deposit', 'avgfill_emissivity'), center=(0,0,0))
slice.set_zlim(('deposit', 'avgfill_emissivity'), 1E-36, 1E-32)
slice.annotate_grids()
slice.zoom(10)
In [ ]:
logging.getLogger('yt').setLevel(logging.INFO)
proj = yt.ProjectionPlot(ds, 'x', ('deposit', 'avgfill_emissivity'), center=(0,0,0))
proj.set_zlim(('deposit', 'avgfill_emissivity'), 1E-36, 1E-32)
#proj.annotate_grids()
proj.zoom(10)
proj.show()
In [ ]:
proj.set_zlim(('deposit', 'avgfill_emissivity'), 1E-14, 1E-10)
proj.save('/home/ychen/d9/FLASH4/stampede/0529_L45_M10_b1_h1/0620_avgfill_deposited_emissivity2.png')
In [ ]:
#logging.getLogger('yt').setLevel(logging.INFO)
#proj = yt.ProjectionPlot(ds, 'x', ('deposit', 'emissivity_nn'), center=(0,0,0))
#proj.set_zlim(('deposit', 'avgfill_emissivity'), 1E-36, 1E-32)
proj.set_zlim(('deposit', 'emissivity_nn'), 1E-14, 1E-10)
#proj.annotate_grids()
#proj.zoom(8)
proj.show()
In [ ]:
ds.add_deposited_particle_field(('io', 'particle_emissivity'), 'nearest')
def _synchrotron_emissivity_pres(field, data):
return data['deposit', 'io_nn_emissivity']*data['gas', 'pressure']
fn = ds.add_field(('deposit', 'emissivity_nn'), function=_synchrotron_emissivity_pres,
display_name='Emissivity', units='erg/s/cm**3/Hz', take_log=True, force_override=True)
In [ ]:
fn = 'emissivity_nn'
slice = yt.SlicePlot(ds, 'x', fn, center=(0,0,0))
slice.set_zlim(fn, 1E-36, 1E-32)
slice.annotate_grids()
slice.zoom(10)