In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import yt
import logging
logging.getLogger('yt').setLevel(logging.ERROR)
import numpy as np
from yt import derived_field
nu = (1500, 'MHz')
qfname = 'nn_emissivity_q_lobe_%.1f%s' % nu
qnfname = ('flash', 'nn_emissivity_qn_lobe_%.1f%s' % nu)
@derived_field(name=qnfname, sampling_type='cell')
def _qn_emissivity(field, data):
#re = np.nan_to_num(data[qfname])
re = data[qfname]
re[np.isinf(re)] = 0.0
return re
ufname = 'nn_emissivity_u_lobe_%.1f%s' % nu
unfname = ('flash', 'nn_emissivity_un_lobe_%.1f%s' % nu)
@derived_field(name=unfname, sampling_type='cell')
def _un_emissivity(field, data):
#re = np.nan_to_num(data[ufname])
re = data[ufname]
re[np.isinf(re)] = 0.0
return re
In [ ]:
proj_axis = 'x'
if type(proj_axis) is str:
postfix = ('_synchrotron_%%.1f%%s_%s' % proj_axis) % nu
elif type(proj_axis) is list:
postfix = ('_synchrotron_%%.1f%%s_%.1f_%.1f_%.1f' % tuple(proj_axis)) % nu
dir = '/home/ychen/d9/FLASH4/2015_production_runs/0529_L45_M10_b1_h1'
fname = dir + '/MHD_Jet_hdf5_plt_cnt_0630' + postfix
ds = yt.load(fname)
In [ ]:
ptype = 'lobe'
norm = yt.YTQuantity(*nu).in_units('GHz').value**0.5
fields = []
for pol in ['i', 'qn', 'un']:
field = ('flash', ('nn_emissivity_%s_%s_%%.1f%%s' % (pol, ptype)) % nu)
fields.append(field)
In [ ]:
slc = yt.SlicePlot(ds, 'x', fields[1][1],center=[0,0,0], width=((50,'kpc'),(80,'kpc')))
slc.set_cmap(fields[1][1], cmap='seismic')
slc.show()
In [ ]:
ad = ds.all_data()
ad[fields[1]][np.isinf(ad[fields[1]])]
In [ ]:
ad[fields[1]][np.isinf(ad[fields[1]])] = 0.0
null = plt.hist(ad[fields[1]], bins=100)
plt.ylim(0,100)
In [ ]:
#%pdb
proj = yt.ProjectionPlot(ds, 'x', fields, center=[0,0,0], width=((50,'kpc'),(80,'kpc')))
for field in fields:
if 'nn_emissivity_i' in field[1]:
proj.set_zlim(field, 1E-3/norm, 1E1/norm)
cmap = plt.cm.hot
cmap.set_bad('k')
proj.set_cmap(field, cmap)
else:
cmap = plt.cm.seismic
proj.set_cmap(field, cmap)
proj.set_zlim(field, -1E0/norm, 1E0/norm)
proj.annotate_timestamp(corner='upper_left', time_format="{time:6.3f} {units}",
time_unit='Myr', draw_inset_box=True)
proj.show()
In [ ]:
frb_I = proj.frb.data[fields[0]].v
frb_Q = proj.frb.data[fields[1]].v
frb_U = proj.frb.data[fields[2]].v
proj.annotate_polline(frb_I, frb_Q, frb_U, factor=16)
proj.show()
In [ ]:
#print(np.nonzero(np.nan_to_num(ds.all_data()[fields[1]])))
adq = ds.all_data()[fields[2]]
print(adq.min(), adq.max())
null = plt.hist(adq, bins=100, range=(2E-20, 1E-18))
In [ ]:
print(proj_axis)
north_vector = [0,0,1]
prj = yt.OffAxisProjectionPlot(ds, proj_axis,fields,width=(60, 'kpc'),
north_vector=north_vector)
In [ ]:
prj.set_zlim(fields[0], 1E-3/norm, 1E1/norm)
cmap = plt.cm.hot
cmap.set_bad('k')
prj.set_cmap(fields[0], cmap)
prj.annotate_clear()
prj.set_cmap(fields[1], 'seismic')
prj.set_zlim(fields[1], -10, 10)
prj.set_cmap(fields[2], 'seismic')
prj.set_zlim(fields[2], -10, 10)
prj.show()
In [ ]:
frb_I = prj.frb.data[fields[0]].v
frb_Q = prj.frb.data[fields[1]].v
frb_U = prj.frb.data[fields[2]].v
prj.annotate_clear()
prj.annotate_polline(frb_I, frb_Q, frb_U, factor=16)
prj.show()
In [ ]:
#null = plt.hist(np.log10(arri.flatten()), range=(-15,3), bins=100)
frb_I = prj.frb.data[fields[0]].v
frb_Q = prj.frb.data[fields[1]].v
frb_U = prj.frb.data[fields[2]].v
factor = 1
nx = 800//factor
ny = 800//factor
I_bin = frb_I.reshape(nx, factor, ny, factor).sum(3).sum(1)
Q_bin = frb_Q.reshape(nx, factor, ny, factor).sum(3).sum(1)
U_bin = frb_U.reshape(nx, factor, ny, factor).sum(3).sum(1)
psi = 0.5*np.arctan2(U_bin, Q_bin)
frac = np.sqrt(Q_bin**2+U_bin**2)/I_bin
#print len(psi.flatten()[psi.flatten().nonzero()])
#null = plt.hist(np.abs(psi.flatten())[psi.flatten().nonzero()], bins=50)
#plt.xlim(0, 0.5*np.pi)
In [ ]:
fig = plt.figure(figsize=(16,4))
ax1 = fig.add_subplot(131)
null = ax1.hist(frac[I_bin.nonzero()].flatten()*100, range=(0,80), bins=40)
ax1.set_xlabel('Polarization fraction (%)')
#ax1.set_xlim(0, 80)
ax2 = fig.add_subplot(132)
null = ax2.hist(psi[I_bin.nonzero()].flatten(), bins=50, range=(-0.5*np.pi, 0.5*np.pi))
x_tick = np.linspace(-0.5, 0.5, 5, endpoint=True)
x_label = [r"$-\pi/2$", r"$-\pi/4$", r"$0$", r"$+\pi/4$", r"$+\pi/2$"]
ax2.set_xlim(-0.5*np.pi, 0.5*np.pi)
ax2.set_xticks(x_tick*np.pi)
ax2.set_xticklabels(x_label)
ax2.set_title(ds.basename + ' %.1f %s' % nu)
ax3 = fig.add_subplot(133)
null = ax3.hist(np.abs(psi[I_bin.nonzero()].flatten()), bins=25, range=(0.0, 0.5*np.pi))
ax3.set_xlim(0.0, 0.5*np.pi)
ax3.set_xticks([x_tick[2:]*np.pi])
ax3.set_xticks(x_tick[2:]*np.pi)
ax3.set_xticklabels(x_label[2:])
In [ ]:
frb_I = prj.frb.data[fields[0]].v
frb_Q = prj.frb.data[fields[1]].v
frb_U = prj.frb.data[fields[2]].v
factor = 1
nx = 800//factor
ny = 800//factor
I_bin_150 = frb_I.reshape(nx, factor, ny, factor).sum(3).sum(1)
Q_bin_150 = frb_Q.reshape(nx, factor, ny, factor).sum(3).sum(1)
U_bin_150 = frb_U.reshape(nx, factor, ny, factor).sum(3).sum(1)
In [ ]:
def plot_polarization_histogram(I_bin, Q_bin, U_bin, fig=None, label=None):
psi = 0.5*np.arctan2(U_bin, Q_bin)
frac = np.sqrt(Q_bin**2+U_bin**2)/I_bin
if not fig:
fig = plt.figure(figsize=(16,4))
ax1 = fig.axes[0]
null = ax1.hist(frac[I_bin.nonzero()].flatten()*100, range=(0,80), bins=40, alpha=0.5)
ax1.set_xlabel('Polarization fraction (%)')
#ax1.set_xlim(0, 80)
ax2 = fig.axes[1]
null = ax2.hist(psi[I_bin.nonzero()].flatten(), bins=50, range=(-0.5*np.pi, 0.5*np.pi), alpha=0.5)
x_tick = np.linspace(-0.5, 0.5, 5, endpoint=True)
x_label = [r"$-\pi/2$", r"$-\pi/4$", r"$0$", r"$+\pi/4$", r"$+\pi/2$"]
ax2.set_xlim(-0.5*np.pi, 0.5*np.pi)
ax2.set_xticks(x_tick*np.pi)
ax2.set_xticklabels(x_label)
#ax2.set_title(ds.basename + ' %.1f %s' % nu)
ax3 = fig.axes[2]
null = ax3.hist(np.abs(psi[I_bin.nonzero()].flatten()), bins=25, range=(0.0, 0.5*np.pi), alpha=0.5,\
label=label)
ax3.legend()
ax3.set_xlim(0.0, 0.5*np.pi)
ax3.set_xticks([x_tick[2:]*np.pi])
ax3.set_xticks(x_tick[2:]*np.pi)
ax3.set_xticklabels(x_label[2:])
return fig
In [ ]:
fig = plt.figure(figsize=(16,4))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)
fig = plot_polarization_histogram(I_bin_150, Q_bin_150, U_bin_150, fig=fig, label='150MHz')
fig = plot_polarization_histogram(I_bin, Q_bin, U_bin, fig=fig, label='1500MHz')
In [ ]: