In [ ]:
import yt
import logging
import numpy as np
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
from yt_cluster_ratio_fields import *
logging.getLogger('yt').setLevel(logging.ERROR)
from mpl_toolkits.axes_grid1 import AxesGrid
In [ ]:
ds = yt.load('/home/ychen/Mount/gdrive/2015_production_runs/1022_L45_M10_b1_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0630')
print(ds)
sl = yt.SlicePlot(ds, 'y', 'entropy_ratio').zoom(4)
sl.set_zlim('entropy_ratio', 0.5, 2)
sl.set_cmap('entropy_ratio', 'seismic')
sl.annotate_timestamp(corner='upper_left', time_format='{time:6.0f} {units}',
time_unit='Myr', text_args={'color': 'k'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.8, 0.8),
plot_args={'colors': 'k', 'lw': 1, 'linestyles': ':'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.6, 0.6),
plot_args={'colors': 'grey', 'lw': 1, 'linestyles': '--'})
sl._setup_plots()
pl = sl.plots['entropy_ratio']
pl.cb.set_ticks([0.5,1,2])
pl.cb.set_ticklabels([0.5,1,2])
sl.show()
In [ ]:
ds = yt.load('/home/ychen/Mount/gdrive/2015_production_runs/1022_L45_M10_b1_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_1062')
print(ds)
sl = yt.SlicePlot(ds, 'y', 'entropy_ratio').zoom(4)
sl.set_zlim('entropy_ratio', 0.5, 2)
sl.set_cmap('entropy_ratio', 'seismic')
sl.annotate_timestamp(corner='upper_left', time_format='{time:6.0f} {units}',
time_unit='Myr', text_args={'color': 'k'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.8, 0.8),
plot_args={'colors': 'k', 'lw': 1, 'linestyles': ':'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.6, 0.6),
plot_args={'colors': 'grey', 'lw': 1, 'linestyles': '--'})
sl._setup_plots()
pl = sl.plots['entropy_ratio']
pl.cb.set_ticks([0.5,1,2])
pl.cb.set_ticklabels([0.5,1,2])
sl.show()
In [ ]:
ds = yt.load('/home/ychen/Mount/gdrive/2015_production_runs/1022_L45_M10_b1_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_1380')
print(ds)
sl = yt.SlicePlot(ds, 'y', 'entropy_ratio', width=[(300, 'kpc'), (150, 'kpc')])
sl.set_zlim('entropy_ratio', 0.5, 2)
sl.set_cmap('entropy_ratio', 'seismic')
sl.annotate_timestamp(corner='upper_left', time_format='{time:6.0f} {units}',
time_unit='Myr', text_args={'color': 'k'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.8, 0.8),
plot_args={'colors': 'k', 'lw': 1, 'linestyles': ':'})
sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.6, 0.6),
plot_args={'colors': 'grey', 'lw': 1, 'linestyles': '--'})
sl._setup_plots()
pl = sl.plots['entropy_ratio']
pl.cb.set_ticks([0.5,1,2])
pl.cb.set_ticklabels([0.5,1,2])
sl.show()
In [ ]:
files = ['/home/ychen/temp/1022_L45_M10_b1_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_%04i' % ind
for ind in [630, 1062, 1380]]
files
In [ ]:
center=(0.0,0.0,0.0)
fig = plt.figure()
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
nrows_ncols = (3, 1),
axes_pad = 0.05,
label_mode = "L",
share_all = True,
cbar_location="right",
cbar_mode="single",
cbar_size="1.5%",
cbar_pad="0%")
cbs = []
for i, f in enumerate(files):
# Load the data and create a single plot
ds = yt.load(f)# load data
print(ds)
sl = yt.SlicePlot(ds, 'y', 'entropy_ratio', width=[(260, 'kpc'), (120, 'kpc')])
sl.set_zlim('entropy_ratio', 0.5, 2)
sl.set_cmap('entropy_ratio', 'seismic')
sl.set_buff_size((800,400))
sl.annotate_timestamp(corner='upper_right', time_format='{time:6.0f} {units}',
time_unit='Myr', text_args={'color': 'k'})
sl.annotate_contour('jet ', ncont=1, clim=(1E-3, 1E-3),
plot_args={'colors': 'yellow', 'lw': 0.1, 'alpha': 0.8, 'linestyles': '-'})
#sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.8, 0.8),
# plot_args={'colors': 'k', 'lw': 0.5, 'linestyles': ':'})
#sl.annotate_contour(field='entropy_ratio', ncont=1, clim=(0.6, 0.6),
# plot_args={'colors': 'grey', 'lw': 0.5, 'linestyles': '--'})
sl.set_font_size(9)
pl = sl.plots['entropy_ratio']
#sl.annotate_text((0.05, 0.95), labels[i], coord_system='axis', text_args={'color':'k'})
# This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
pl.figure = fig
pl.axes = grid[i].axes
pl.cax = grid.cbar_axes[i]
cbs.append(pl.cb)
# Finally, this actually redraws the plot.
sl._setup_plots()
#pl.cax.yaxis.set_ticks([0.5,1,2])
#pl.cax.yaxis.set_ticklabels([0.5,1,2])
#sl.plots[field].axes.set_xticks([-10,0,10,20])
#sl.plots[field].axes.set_xticklabels([-10,0,10,20])
In [ ]:
sl._colorbar_label[('gas', 'entropy_ratio')] = [0.5, 1, 2]
sl._colorbar_label[('gas', 'entropy_ratio')]
In [ ]:
for cax in grid.cbar_axes:
print(cax.)
print(cax.get_yticks())
In [ ]:
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(-55, 55)
ax.set_yticks([-50,0,50])
ax.set_yticklabels([-50,0,50])
ax.set_xlim(-130, 130)
ax.set_xticks([-100, -50, 0, 50, 100])
ax.set_xticklabels([-100, -50, 0, 50, 100])
#ax.minorticks_off()
ax.tick_params(axis='x', color='grey')
ax.tick_params(axis='y', color='grey')
ax.grid(ls='--', alpha=0.5)
if i == 0:
ax.set_xlabel('z (kpc)')
ax.set_ylabel('x (kpc)')
else:
ax.set_xlabel('')
ax.set_ylabel('x (kpc)')
fig.set_figwidth(5)
fig.set_figheight(8)
fig.savefig('entropy_ratio_slices.pdf', bbox_inches='tight')
fig
In [ ]:
center = [0.0,0.0,0.0]
radius = (200, 'kpc')
sp = ds.sphere(center, radius)
low_entropy = sp.cut_region(["obj['entropy_ratio'] < 0.8"])
In [ ]:
sl = yt.SlicePlot(ds, 'y', 'entropy_ratio', data_source=low_entropy).zoom(4)
sl.set_zlim('entropy_ratio', 0.5, 2)
sl.set_cmap('entropy_ratio', 'seismic')
sl.annotate_sphere(center, radius)
sl.show()
In [ ]:
sl = yt.SlicePlot(ds, 'y', 'temperature_ratio', data_source=low_entropy).zoom(4)
sl.set_zlim('temperature_ratio', 0.5, 2)
sl.set_cmap('temperature_ratio', 'seismic')
sl.annotate_sphere(center, radius)
sl.show()
In [ ]:
plt.scatter(low_entropy['temperature_ratio'], low_entropy['entropy_ratio'], s=1, lw=0)
plt.xlabel('temperature_ratio')
plt.ylabel('entropy_ratio')
In [ ]:
print(sum(low_entropy['cell_mass'].in_units('Msun')))
In [ ]:
print(sum(sp['cell_mass'].in_units('Msun')))
In [ ]:
'%e' % 3590553354757.157
In [ ]:
mH = yt.physical_constants.mass_hydrogen
kB = yt.physical_constants.kb
In [ ]:
5/2*sum(low_entropy['cell_mass']*low_entropy['temperature']*(1.0/low_entropy['temperature_ratio']-1))/0.61/mH*kB
In [ ]:
table = np.genfromtxt('/home/ychen/data/0only_1022_h1_10Myr/gridanalysis_entropy_ratio.txt')
table.shape
tt = table[:,1]
for i in np.arange(table.shape[1]-1, 1, -1):
if max(table[:,i]) > 0:
plt.fill_between(tt, table[:,i], 1E7, label='x=%.1f' % (i/10))
plt.xlabel('time (Myr)')
plt.ylabel(r'mass with entropy ratio < x (M$_{\odot}$)')
#plt.semilogy()
#plt.ylim(1E7, 1E11)
plt.ylim(0, 3.5E10)
plt.xlim(0,500)
plt.legend(loc=2, fontsize=7)
In [ ]:
ls = {9: 'solid', 8: 'dashed', 7: 'dotted', 6: 'solid', 5: 'dashed', 4: 'dotted', 3: 'solid'}
plt.figure(figsize=(6,3))
E_jet = 1E45*10*yt.units.Myr.in_units('s').v
print(E_jet)
table = np.genfromtxt('/home/ychen/data/0only_1022_h1_10Myr/gridanalysis_entropy_ratio_energy.txt')
tt = table[:,1]
print(table.shape)
for i in np.arange(table.shape[1]-1, 1, -1):
if max(table[:,i]) > 0:
plt.fill_between(tt, table[:,i], 1E-3*E_jet, label='x=%.1f' % (i/10), linestyle=ls[i])
plt.xlabel('time (Myr)')
plt.ylabel('energy budget of low-entropy gas (erg)')
#plt.semilogy()
plt.ylim(1E-3*E_jet, 0.7E0*E_jet)
plt.legend(loc=2)
plt.twinx()
#plt.semilogy()
plt.ylim(1E-3, 0.7E0)
plt.ylabel(r'(E$_{\mathrm{jet}}$)')
plt.xlim(0,300)
In [ ]:
Etots_arr = np.array(Etots)
Emags_arr = np.array(Emags)
Eths_arr = np.array(Eths)
ts_arr = np.array(ts)
Ps = (Etots_arr[1:-1]-Etots_arr[0:-2])/(ts_arr[1:-1]-ts_arr[0:-2])
plt.plot(ts[:-2], Ps)
#plt.plot(ts[:-2], (Eths_arr/Emags_arr/100.0)[:-2])
plt.ylabel('erg')
plt.xlabel('s')
print Ps
In [ ]:
plt.plot(ts, Emags)
In [ ]: