In [ ]:

    

%matplotlib inline
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 200
import matplotlib.pyplot as plt
import os
import yt
yt.mylog.setLevel("WARNING")
import numpy as np
from yt_synchrotron_emissivity import *
from scipy.ndimage import gaussian_filter
import pyfits


    

In [ ]:

    

dir = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/'
proj_axis = 'x'
zoom_fac = 6
nus = [(150, 'MHz'), (1400, 'MHz')]

ptype = 'lobe'
extend_cells = 32
res = (256, 128)
sigma = 1


ds = yt.load(dir + 'data/MHD_Jet_10Myr_hdf5_plt_cnt_0910')
print(ds.basename, ds.current_time.in_units('Myr'))
width = ds.domain_width[[2,1]]/zoom_fac

dirs = ['fits_linear', 'fits_log', 'fits_emis', 'fits_gamc']

fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/%s/synchrotron_lobe_0910.fits'
for dir in dirs:
    hdulist = pyfits.open(fitsname % dir)
    frb_I = {}
    for nu in nus:
        stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
        frb_I[nu] = hdulist[stokes.I[1]].data
    header = hdulist[stokes.I[1]].header
    xr = -header['CRPIX1']*header['CDELT1'] + header['CRVAL1']
    xl = (header['NAXIS1'] - header['CRPIX1'])*header['CDELT1'] + header['CRVAL1']
    yr = -header['CRPIX2']*header['CDELT2'] + header['CRVAL2']
    yl = (header['NAXIS2'] - header['CRPIX2'])*header['CDELT2'] + header['CRVAL2']
    ext = ds.arr([yr, yl, xr, xl], input_units='cm').in_units('kpc')

    nu1, nu2 = nus

    I1 = gaussian_filter(frb_I[nu1], sigma)
    I2 = gaussian_filter(frb_I[nu2], sigma)
    alpha = np.log10(I2/I1)/np.log10(nu2[0]/nu1[0])
    alpha = np.ma.masked_where(I2<2E-3, np.array(alpha))
    #ext = ds.arr([-0.5*width[0], 0.5*width[0], -0.5*width[1], 0.5*width[1]]).in_units('kpc')

    fig = plt.figure(figsize=(8,3))
    ax = fig.add_subplot(111)
    cmap = plt.cm.jet
    im = ax.imshow(alpha.transpose(), cmap=cmap, vmin=-2, vmax=-0.5, extent=ext, origin='lower', aspect='equal')
    ax.set_facecolor('navy')
    cbar = fig.colorbar(im)
    cbar.ax.tick_params(direction='in')
    #cbar.ax.set_yticks([-0.5, -0])


    #ax.annotate(labels[i % len(dir)], (0.65, 0.75) , xycoords='axes fraction', color='white')

    timestamp = '%.1f Myr' % ds.current_time.in_units('Myr')
    ax.annotate(timestamp, (0.04, 0.75) , xycoords='axes fraction', color='white')
    ax.text(40, 25, dir.split('_')[-1], color='white')


    

In [ ]:

    

dir = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/'
ds = yt.load(dir + 'data/MHD_Jet_10Myr_hdf5_plt_cnt_0910')
#dir = '/d/d5/ychen/2015_production_runs/0529_h1_00/'
#ds = yt.load(dir + 'data/MHD_Jet_hdf5_plt_cnt_0600')

proj_axis = 'x'
zoom_fac = 6
nus = [(150, 'MHz'), (1400, 'MHz')]

ptype = 'lobe'
extend_cells = 8
res = (512, 256)
sigma = 1


print(ds.basename, ds.current_time.in_units('Myr'))
width = ds.domain_width[[2,1]]/zoom_fac

fitsname = synchrotron_fits_filename(ds, dir, ptype, proj_axis)
hdulist = pyfits.open(fitsname)
frb_I = {}
for nu in nus:
    stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
    frb_I[nu] = hdulist[stokes.I[1]].data
header = hdulist[stokes.I[1]].header
xr = -header['CRPIX1']*header['CDELT1'] + header['CRVAL1']
xl = (header['NAXIS1'] - header['CRPIX1'])*header['CDELT1'] + header['CRVAL1']
yr = -header['CRPIX2']*header['CDELT2'] + header['CRVAL2']
yl = (header['NAXIS2'] - header['CRPIX2'])*header['CDELT2'] + header['CRVAL2']
ext = ds.arr([yr, yl, xr, xl], input_units='cm').in_units('kpc')

#         ds_sync = yt.load(synchrotron_filename(ds, extend_cells=extend_cells))

#         # Setting up units and coordinates (we want z-y figures)
#         ds_sync.field_list
#         ds_sync.coordinates.x_axis['x'] = 2
#         ds_sync.coordinates.x_axis[0] = 2
#         ds_sync.coordinates.y_axis['x'] = 1
#         ds_sync.coordinates.y_axis[0] = 1
#         frb_I = {}
#         for nu in nus:
#             stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
#             if stokes.I not in ds_sync.field_list: continue
#             if proj_axis in ['x','y','z']:
#                 p = yt.ProjectionPlot(ds_sync, proj_axis, stokes.I, center=[0,0,0], width=width, max_level=6)
#                 frb_I[nu] = p.frb.data[stokes.I].v
#             else:
#                 p = yt.OffAxisProjectionPlot(ds_sync, proj_axis, stokes.I, center=[0,0,0], width=width, north_vector=[0,1,0])
#                 frb_I[nu] = p.frb.data[stokes.I].v

nu1, nu2 = nus

I1 = gaussian_filter(frb_I[nu1], sigma)
I2 = gaussian_filter(frb_I[nu2], sigma)
alpha = np.log10(I2/I1)/np.log10(nu2[0]/nu1[0])
alpha = np.ma.masked_where(I2<1E-7, np.array(alpha))
#ext = ds.arr([-0.5*width[0], 0.5*width[0], -0.5*width[1], 0.5*width[1]]).in_units('kpc')

fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(111)
cmap = plt.cm.jet
im = ax.imshow(alpha.transpose(), cmap=cmap, vmin=-2, vmax=-0.5, extent=ext, origin='lower', aspect='equal')
ax.set_facecolor('navy')
cbar = fig.colorbar(im)
cbar.ax.tick_params(direction='in')
#cbar.ax.set_yticks([-0.5, -0])


#ax.annotate(labels[i % len(dir)], (0.65, 0.75) , xycoords='axes fraction', color='white')

timestamp = '%.1f Myr' % ds.current_time.in_units('Myr')
ax.annotate(timestamp, (0.04, 0.75) , xycoords='axes fraction', color='white')


    

In [ ]:

    

ptype = 'lobe'
proj_axis = 'x'
nu = (150, 'MHz')
zoom_fac = 6
stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
field = stokes.I[1]

dirs = ['fits_linear', 'fits_log', 'fits_emis', 'fits_gamc']

fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/%s/synchrotron_lobe_0910.fits'
for dir in dirs:
    ds_fits = yt.load(fitsname % dir)

    
    ds_fits.coordinates.x_axis['z'] = 1
    ds_fits.coordinates.x_axis[2] = 1
    ds_fits.coordinates.y_axis['z'] = 0
    ds_fits.coordinates.y_axis[2] = 0
    p = yt.SlicePlot(ds_fits, 'z', field)
    p.set_log(field, True)
    p.set_zlim(field, 1E-2, 1E2)
    p.annotate_text((40, 25), dir.split('_')[-1], coord_system='plot', text_args={'color':'black'})
    p.show()


    

In [ ]:

    

ptype = 'lobe'
proj_axis = 'x'
nu = (150, 'MHz')
zoom_fac = 6
stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
field = stokes.I[1]
fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/fits_linear/synchrotron_lobe_0910.fits'

ds_fits = yt.load(fitsname)
ds_fits.coordinates.x_axis['z'] = 1
ds_fits.coordinates.x_axis[2] = 1
ds_fits.coordinates.y_axis['z'] = 0
ds_fits.coordinates.y_axis[2] = 0
p = yt.SlicePlot(ds_fits, 'z', field)
p.set_log(field, True)
p.set_zlim(field, 1E-2, 1E2)
p.show()


    

In [ ]:

    

ptype = 'lobe'
proj_axis = 'x'
nu = (150, 'MHz')
zoom_fac = 6
stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
field = stokes.I[1]
fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/fits/synchrotron_lobe_0910.fits'

ds_fits = yt.load(fitsname)
ds_fits.coordinates.x_axis['z'] = 1
ds_fits.coordinates.x_axis[2] = 1
ds_fits.coordinates.y_axis['z'] = 0
ds_fits.coordinates.y_axis[2] = 0
p = yt.SlicePlot(ds_fits, 'z', field)
p.set_log(field, True)
p.set_zlim(field, 1E-3, 1E-1)
p.show()


    

In [ ]: