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

nus = [(1500, 'MHz'), (150, 'MHz')]
proj_axis = 'x'

prj = {}
frb_I, frb_Q, frb_U = {}, {}, {}
I_bin, Q_bin, U_bin = {}, {}, {}
psi, frac = {}, {}

nu0 = nus[0]
nu1 = nus[1]
nus_str = ['%.1f %s' % nu0, '%.1f %s' % nu1]

for nu in nus:
    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


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

        north_vector = [0,0,1]
    if proj_axis == 'x':
        prj[nu] = yt.ProjectionPlot(ds, proj_axis, fields, width=(80, 'kpc'))
    else:
        prj[nu] = yt.OffAxisProjectionPlot(ds, proj_axis, fields,width=(60, 'kpc'),
                                   north_vector=north_vector)
    frb_I[nu] = prj[nu].frb.data[fields[0]].v
    frb_Q[nu] = prj[nu].frb.data[fields[1]].v
    frb_U[nu] = prj[nu].frb.data[fields[2]].v


    

In [ ]:

    

-


    

In [ ]:

    

# Binning
for nu in nus:
    factor = 1
    nx = 800//factor
    ny = 800//factor

    I_bin[nu] = frb_I[nu].reshape(nx, factor, ny, factor).sum(3).sum(1)
    Q_bin[nu] = frb_Q[nu].reshape(nx, factor, ny, factor).sum(3).sum(1)
    U_bin[nu] = frb_U[nu].reshape(nx, factor, ny, factor).sum(3).sum(1)

    psi[nu] = 0.5*np.arctan2(U_bin[nu], Q_bin[nu])
    frac[nu] = np.sqrt(Q_bin[nu]**2+U_bin[nu]**2)/I_bin[nu]


    

In [ ]:

    

fig = plt.figure(figsize=(16,4))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)

for nu in reversed(nus):
    nu_str = '%.0f%s' % nu
    fig = plot_polarization_histogram(frac[nu], psi[nu], I_bin[nu], fig=fig, label=nu_str)


    

In [ ]:

    

ns, bins = {}, {}
for nu in nus:
    frac_nu = frac[nu]
    ns[nu], bins[nu], patches = plt.hist(frac_nu[frac_nu.nonzero()]*100, range=(0,80), bins=40, \
                              weights=I_bin[nu][frac_nu.nonzero()], normed=True)

plt.cla()
plt.scatter((bins[nu0][1:]+bins[nu0][:-1])/2, (ns[nu0]-ns[nu1])/(ns[nu0]+ns[nu1]))
plt.hlines(0, 0, 80, linestyle=':')
plt.ylim(-0.3,1)


    

In [ ]:

    

fig = plt.figure(figsize=(20,8))
ax1 = fig.add_subplot(121)
ax1.imshow(frac[nu0], vmin=0, vmax=0.7)
plt.title(nus_str[0])
ax2 = fig.add_subplot(122)
img2 = ax2.imshow(frac[nu1], vmin=0, vmax=0.7)
plt.title(nus_str[1])
cb = plt.colorbar(img2, pad=0)
cb.set_label('polarization fraction')


    

In [ ]:

    

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
img = ax.imshow((frac[nu0]-frac[nu1]), cmap='seismic', vmin=-0.5, vmax=0.5)
cb = plt.colorbar(img, pad=0)

cb.set_label('polarization fraction difference (%s - %s)' % tuple(nus_str))


    

In [ ]:

    

fig = plt.figure(figsize=(22,8))
ax1 = fig.add_subplot(121)
ax1.imshow(psi[nu0], vmin=0, vmax=np.pi/2)
plt.title(nus_str[0])
ax2 = fig.add_subplot(122)
img2 = ax2.imshow(psi[nu1], vmin=0, vmax=np.pi/2)
plt.title(nus_str[1])
cb = plt.colorbar(img2)
cb.set_label('polarization angle')


    

In [ ]:

    

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
img = plt.imshow((psi[nu0]-psi[nu1]), cmap='seismic', vmin=-1.5, vmax=1.5)
cb = plt.colorbar(img, pad=0)

cb.set_label('polarization angle difference (%s - %s)' % tuple(nus_str))


    

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