In [ ]:
#%pdb
%matplotlib inline
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
import os
import yt
yt.mylog.setLevel("ERROR")
import numpy as np
from synchrotron.yt_synchrotron_emissivity import *
from mpl_toolkits.axes_grid1 import AxesGrid, ImageGrid
import pyfits
from tools import setup_cl
In [ ]:
dir = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/'
ds = yt.load(dir+'data/MHD_Jet_10Myr_hdf5_plt_cnt_0100')
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 = synchrotron_fits_filename(ds, dir, ptype, proj_axis)
ds = yt.load(os.path.join(dir, 'data/MHD_Jet_10Myr_hdf5_plt_cnt_0100'))
width = ds.domain_width[[2,1]].in_units('kpc')/zoom_fac
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, width=width)
p.set_log(field, True)
p.set_zlim(field, 1E-2, 1E2)
p.show()
In [ ]:
#dirs = ['/home/ychen/data/00only_0605_hinf/',\
# '/home/ychen/data/00only_0529_h1/',\
# '/home/ychen/data/00only_0605_h0/',]
def plot_synchrotron_imagegrid(proj_axis, nu):
dirs = ['/d/d5/ychen/2015_production_runs/0204_hinf_10Myr/',\
'/d/d5/ychen/2015_production_runs/1022_h1_10Myr/',\
'/d/d5/ychen/2016_production_runs/1212_h0_10Myr/',]
filenumbers = [100, 200, 600, 910, 1050]
iterator = []
for filenumber in filenumbers:
for dir in dirs:
iterator.append((filenumber, dir))
labels = ['toroidal', 'helical', 'poloidal']
ptype = 'lobe'
zoom_fac = 6
extend_cells = 32
res = (256, 128)
stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
fig = plt.figure()
grid = ImageGrid(fig, (0.075,0.05,0.85,0.90),
nrows_ncols = (len(filenumbers), len(dirs)),
axes_pad = 0.05,
label_mode = "L",
share_all = True,
cbar_location="right",
cbar_mode="single",
cbar_size="2%",
cbar_pad="0%")
field = stokes.I[1]
for i, (filenumber, dir) in enumerate(iterator):
norm = yt.YTQuantity(*nu).in_units('GHz').value**0.5
# Load the data and create a single plot
ds = yt.load(os.path.join(dir, 'data/MHD_Jet_10Myr_hdf5_plt_cnt_%04d' % filenumber))
print(dir, filenumber, ds.current_time.in_units('Myr'))
width = ds.domain_width[[2,1]].in_units('kpc')/zoom_fac
fitsname = synchrotron_fits_filename(ds, dir, ptype, proj_axis)
if not os.path.isfile(fitsname): continue
ds_fits = yt.load(fitsname)
ds_fits.current_time = ds.current_time
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, width=width)
# ds_sync = yt.load(synchrotron_file_name(ds, extend_cells=extend_cells))
# # Setting up units and coordinates (we want z-y figures)
# ds_sync.field_list
# ds_sync.field_info[field].units = 'Jy/cm/arcsec**2'
# ds_sync.field_info[field].output_units = 'Jy/cm/arcsec**2'
# 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
#
# if field not in ds_sync.field_list: continue
# if proj_axis in ['x','y','z']:
# p = yt.ProjectionPlot(ds_sync, proj_axis, field, center=[0,0,0], width=width, max_level=6)
# else:
# p = yt.OffAxisProjectionPlot(ds_sync, proj_axis, field, center=[0,0,0], width=width, north_vector=[0,1,0])
#p.set_buff_size(res)
p.set_figure_size(14)
p.set_axes_unit('kpc')
p.set_log(field, True)
# Colormap and colorbar limits
# Ensure the colorbar limits match for all plots
p.set_zlim(field, 1E-3/norm, 1E1/norm)
cmap = plt.cm.hot
cmap.set_bad('k')
p.set_cmap(field, cmap)
if i // len(dirs) == 0:
p.annotate_text((0.61, 0.80), labels[i], coord_system='axis', text_args={'color':'w'})
#p.set_font_size(9)
if i % len(dirs) == 0:
p.annotate_timestamp(0.02, 0.80, time_format="{time:6.1f} {units}", time_unit='Myr', text_args={'color':'w'})
#p.annotate_timestamp(corner='upper_left', time_format="{time:6.3f} {units}",
# time_unit='Myr', draw_inset_box=True)
# Annotate polline (needs projection of Stokes IQU)
#p._recreate_frb()
#frb_I = p.frb.data[stokes.I].v
#frb_Q = p.frb.data[stokes.Q].v
#frb_U = p.frb.data[stokes.U].v
#p.annotate_clear()
#p.annotate_polline(frb_I[nu], frb_Q, frb_U, factor=25, scale=15)
# This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
plot = p.plots[field]
plot.figure = fig
plot.axes = grid[i].axes
plot.cax = grid.cbar_axes[i]
# Finally, this actually redraws the plot.
p._setup_plots()
return fig
In [ ]:
fig = plot_synchrotron_imagegrid('x', (150, 'MHz'))
In [ ]:
clabel = r'$\rm{Projected }$ $\rm{Emissivity\ I\ Lobe\ 150Mhz\ x$$\ \ \left(\frac{\rm{Jy}}{\rm{arcsec}^{2}}\right)$'
cax = grid.cbar_axes[0]
cax.set_ylabel(clabel)
for i, ax in enumerate(grid.axes_all):
ax.tick_params(axis='x', color='grey')
ax.tick_params(axis='y', color='grey')
ax.set_ylim(-40, 40)
ax.set_yticks([-25,0,25])
ax.set_yticklabels([-25,0,25])
ax.set_xticks([-75,-50,-25,0,25,50,75])
ax.set_xticklabels(['',-50,'',0,'',50,''])
ax.tick_params(axis='x', color='grey')
ax.tick_params(axis='y', color='grey')
ax.grid(ls='--', alpha=0.5)
if i != 6:
ax.set_ylabel('')
else:
ax.set_ylabel('y (kpc)')
if i != 13:
ax.set_xlabel('')
else:
ax.set_xlabel('z (kpc)')
#print(ax.get_xlabel())
fig.subplots_adjust(left=0.2, bottom=0.2, right=0.5)
fig.set_figwidth(13)
fig.set_figheight(11)
fig.savefig('synchrotron_150MHz_x.pdf', bbox_inches='tight')
fig
In [ ]:
fig = plot_synchrotron_imagegrid([1,0,2], (150, 'MHz'))
In [ ]:
clabel = r'$\rm{Projected }$ $\rm{Emissivity\ I\ Lobe\ 150Mhz\ [1,0,2]}$$\ \ \left(\frac{\rm{Jy}}{\rm{arcsec}^{2}}\right)$'
cax = grid.cbar_axes[0]
cax.set_ylabel(clabel)
for i, ax in enumerate(grid.axes_all):
ax.tick_params(axis='x', color='grey')
ax.tick_params(axis='y', color='grey')
ax.set_xlim(-60, 60)
ax.set_ylim(-30, 30)
ax.set_yticks([-25,0,25])
ax.set_yticklabels([-25,0,25])
ax.set_xticks([-50,-25,0,25,50])
ax.set_xticklabels([-50,'',0,'',50,])
ax.grid(ls='--', alpha=0.5)
if i != 6:
ax.set_ylabel('')
else:
ax.set_ylabel('y (kpc)')
if i != 13:
ax.set_xlabel('')
else:
ax.set_xlabel('z (kpc)')
#print(ax.get_xlabel())
fig.subplots_adjust(left=0.2, bottom=0.2, right=0.5)
fig.set_figwidth(13)
#fig.set_figheight(10)
#fig.savefig('synchrotron_150MHz_1_0_2.pdf', bbox_inches='tight')
fig
In [ ]:
from scipy.ndimage import gaussian_filter
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
sigma = 1
nu1, nu2 = nus
I1 = gaussian_filter(frb_I[nu1], sigma)
I2 = gaussian_filter(frb_I[nu2], sigma)
alpha = np.log10(I2/I1)/np.log10(1400/150)
alpha = np.ma.masked_where(I2<1E-3, np.array(alpha))
ext = ds.arr([-0.5*width[0], 0.5*width[0], -0.5*width[1], 0.5*width[1]])
plt.figure(figsize=(10,4), dpi=150)
cmap = plt.cm.jet
cmap.set_bad('navy')
plt.imshow(alpha, cmap=cmap, vmin=-1.5, vmax=-0.5, extent=ext.in_units('kpc'), origin='lower', aspect='equal')
plt.xlabel('z (kpc)')
plt.ylabel('y (kpc)')
plt.axes().tick_params(direction='in')
cb = plt.colorbar(fraction=0.20, pad=0, aspect=20)
cb.set_label('Spectral Index (%s) (1.4GHz/150MHz)' % ptype)
cb.ax.tick_params(direction='in')
In [ ]:
#dirs = ['/home/ychen/data/00only_0605_hinf/',\
# '/home/ychen/data/00only_0529_h1/',\
# '/home/ychen/data/00only_0605_h0/',]
dirs = ['/d/d5/ychen/2015_production_runs/0204_hinf_10Myr/',\
'/d/d5/ychen/2015_production_runs/1022_h1_10Myr/',\
'/d/d5/ychen/2015_production_runs/0204_h0_10Myr/',]
filenumbers = [100, 200, 600, 910, 1050]
labels = ['toroidal', 'helical', 'poloidal']
ptype = 'lobe'
proj_axis = [1,0,2]
nu = (1400, 'MHz')
extend_cells = 32
res = (256, 128)
stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
field = stokes.I[1]
data = {}
for j, dir in enumerate(dirs):
times = np.zeros(len(filenumbers))
I_total = np.zeros(len(filenumbers))
for i, filenumber in enumerate(filenumbers):
# Load the data and create a single plot
ds = yt.load(os.path.join(dir, 'data/MHD_Jet_10Myr_hdf5_plt_cnt_%04d' % filenumber))
print(dir, filenumber, ds.current_time.in_units('Myr'))
fitsname = synchrotron_fits_filename(ds, dir, ptype, proj_axis)
if not os.path.isfile(fitsname): continue
ds_fits = yt.load(fitsname)
times[i] = ds.current_time.in_units('Myr')
ad = ds_fits.all_data()
I_total[i] = np.sum(ad[field])
data[dir] = (times, I_total)
In [ ]:
data_30 = data.copy()
In [ ]:
colors, labels = setup_cl(dirs)
for dir, (time, I_total) in data_x.items():
plt.plot(time, I_total, 'o-', label=labels[dir]+' x', c=colors[dir])
for dir, (time, I_total) in data_30.items():
plt.plot(time, I_total, '^:', label=labels[dir]+' 30 deg', c=colors[dir])
plt.legend()
plt.semilogy()
plt.xlabel('time (Myr)')
plt.ylabel('Total Synchrotron Intensity')
plt.title('150MHz')
In [ ]:
colors, labels = setup_cl(dirs)
for dir, (time, I_total) in data_x.items():
plt.plot(time, I_total, 'o-', label=labels[dir]+' x', c=colors[dir])
for dir, (time, I_total) in data_30.items():
plt.plot(time, I_total, '^:', label=labels[dir]+' 30 deg', c=colors[dir])
plt.legend()
plt.semilogy()
plt.xlabel('time (Myr)')
plt.ylabel('Total Synchrotron Intensity')
plt.title('1400MHz')
In [ ]: