In [ ]:
%matplotlib inline
import os
import yt
yt.mylog.setLevel("WARNING")
import numpy as np
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
from matplotlib import colors
from mpl_toolkits.axes_grid1 import make_axes_locatable
from conduction.yt_conduction_fields import *
from yt_cluster_ratio_fields import *
In [ ]:
ds = yt.load('/d/d12/ychen/2016_production_runs/1212_L45_M10_b1_h0_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_1047')
sp1 = ds.sphere([0,0,0], (180, 'kpc'))
sp2 = ds.sphere([0,0,0], (0.5, 'kpc'))
sp = sp1 - sp2
ds.current_time.in_units('Myr')
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'jet ')
plot.zoom(4)
plot.set_zlim('jet ', 1E-7, 1)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'xray_emissivity_0.1_100_keV')
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.zoom(4)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'xray_cooling_time', )
plot.set_unit('xray_cooling_time', 'Gyr')
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.zoom(4)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'temperature')
plot.set_unit('temperature', 'keV', equivalency='thermal')
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.zoom(4)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'entropy', data_source=sp)
plot.set_zlim('entropy', 30, 250)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.zoom(4)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'H_nuclei_density')
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.zoom(4)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'temperature_gradient_magnitude', data_source=sp, width=(240, 'kpc'))
plot.set_unit('temperature_gradient_magnitude', 'K/pc')
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'temperature_gradient_x', data_source=sp, width=(240, 'kpc'))
plot.set_unit('temperature_gradient_x', 'K/pc')
plot.set_cmap('temperature_gradient_x', 'seismic')
plot.set_zlim('temperature_gradient_x', -1E3, 1E3)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True, plot_args={'color': 'white'})
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_conduction_coefficient', data_source=sp, width=(240, 'kpc'))
plot.set_zlim('spitzer_conduction_coefficient', 1E29, 1E32)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_heat_flux_x', data_source=sp, width=(240, 'kpc'))
plot.set_log('spitzer_heat_flux_x' , True, linthresh=1E-3)
plot.set_cmap('spitzer_heat_flux_x', 'seismic')
plot.set_zlim('spitzer_heat_flux_x', -1E-1, 1E-1)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_heat_flux_divergence', data_source=sp, width=(240, 'kpc'))
plot.set_log('spitzer_heat_flux_divergence' , True, linthresh=1E-25)
plot.set_cmap('spitzer_heat_flux_divergence', 'seismic_r')
plot.set_zlim('spitzer_heat_flux_divergence', -1E-23, 1E-23)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_heating_rate', data_source=sp, width=(240, 'kpc'))
plot.set_log('spitzer_heating_rate' , True, linthresh=1E40)
plot.set_cmap('spitzer_heating_rate', 'seismic')
plot.set_zlim('spitzer_heating_rate', -1E42, 1E42)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_heating_rate', data_source=sp, width=(240, 'kpc'))
plot.set_log('spitzer_heating_rate' , True, linthresh=1E40)
plot.set_cmap('spitzer_heating_rate', 'seismic')
plot.set_zlim('spitzer_heating_rate', -1E42, 1E42)
plot.annotate_contour('entropy_ratio', ncont=2, clim=(0.9, 1.0), take_log=False)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'total_heating_rate', data_source=sp, width=(240, 'kpc'))
plot.set_log('total_heating_rate' , True, linthresh=1E40)
plot.set_cmap('total_heating_rate', 'seismic')
plot.set_zlim('total_heating_rate', -1E42, 1E42)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'y', 'total_cooling_time', data_source=sp, width=(240, 'kpc'))
plot.set_log('total_cooling_time' , True, linthresh=1E2)
plot.set_cmap('total_cooling_time', 'seismic')
plot.set_zlim('total_cooling_time', -1E4, 1E4)
plot.annotate_contour('jet ', ncont=2, clim=(1E-3, 1E-2), take_log=True)
plot.show()
In [ ]:
extrema = {'entropy': (20, 250), 'spherical_radius': ((10, 'kpc'), (150, 'kpc'))}
fields = ['spitzer_heating_rate', 'xray_luminosity_0.1_100_keV', 'cell_mass']
prof = yt.create_profile(sp, ['entropy', 'spherical_radius'], fields, weight_field=None,
extrema=extrema)
In [ ]:
pp = yt.PhasePlot.from_profile(prof, figure_size=6)
pp.set_unit('spherical_radius', 'kpc')
pp.set_cmap('spitzer_heating_rate', 'RdBu_r')
pp.set_log('spitzer_heating_rate', False)
pp.set_zlim('spitzer_heating_rate', -1E44, 1E44)
pp.set_ylim(20, 150)
pp.set_xlim(20, 250)
plot2 = pp['spitzer_heating_rate']
r_ticks = [30, 40, 50, 60, 80, 100, 150]
plot2.axes.set_yticks(r_ticks)
#plot2.axes.yaxis.set_ticklabels(r_ticks)
plot2.axes.set_yticklabels(r_ticks)
plot2.axes.set_yticklabels([], minor=True)
#print(plot2.axes.yaxis.get_ticklabels())
entropy_ticks = [20, 30, 40, 60, 80, 100, 150, 200]
plot2.axes.set_xticks(entropy_ticks)
plot2.axes.set_xticklabels(entropy_ticks)
plot2.save('heating_rate_entropy_r_100Myr_poloidal.pdf')
In [ ]:
prof.y
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.step(prof.y.in_units('kpc'), prof['spitzer_heating_rate'].T.sum(axis=1),
label='Spitzer Cond. (heating)', linewidth=1, color='C1')
ax.step(prof.y.in_units('kpc'), -prof['spitzer_heating_rate'].T.sum(axis=1),
label='Spitzer Cond. (cooling)', linestyle='dotted', linewidth=1, color='C1')
ax.step(prof.y.in_units('kpc'), prof['xray_luminosity_0.1_100_keV'].T.sum(axis=1),
label='X-ray Cooling', linestyle='dashed', color='C0')
ax.step(prof.y.in_units('kpc'), 0.1*prof['spitzer_heating_rate'].T.sum(axis=1),
label='10% Spitzer Cond. (heating)', linewidth=2, color='C3')
ax.step(prof.y.in_units('kpc'), -0.1*prof['spitzer_heating_rate'].T.sum(axis=1),
label='10% Spitzer Cond. (cooling)', linestyle='dotted', linewidth=2, color='C3')
ax.semilogx()
ax.semilogy()
ax.grid()
In [ ]:
fig = plt.figure(figsize=(4.5,5.5))
axMesh = fig.add_subplot(111)
divider = make_axes_locatable(axMesh)
axHist = divider.append_axes("bottom", size=1.9, pad=0, sharex=axMesh)
cax = divider.append_axes("right", size="4%", pad=0)
im = axMesh.pcolormesh(prof.x, prof.y.in_units('kpc'), prof['spitzer_heating_rate'].T.v*0.1,
cmap='RdBu_r', norm=colors.SymLogNorm(vmin=-1E43, vmax=1E43, linthresh=1E41))
axMesh.set_xscale('log')
axMesh.set_yscale('log')
r_ticks = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150]
r_ticklabels = [10, 20, 30, 40, 50, 60, '', 80, '', 100, '', '', '', '', 150]
axMesh.set_ylim(10,150)
axMesh.set_yticks([100])
axMesh.set_yticklabels([])
axMesh.set_yticks(r_ticks, minor=True)
axMesh.set_yticklabels(r_ticklabels, minor=True)
axMesh.set_ylabel('radius [kpc]')
axMesh.tick_params(axis='x', which='both', direction='in', labelbottom=False, top='on')
axMesh.tick_params(axis='y', which='both', direction='in')
axMesh.text(23, 115, '%.0f Myr' % ds.current_time.in_units('Myr'))
cb = plt.colorbar(im, cax=cax)
cb.set_ticks([-1E43, -1E42, -1E41, 0, 1E41, 1E42, 1E43])
#cax.set_yticks([np.linspace(-9E43, -2E42, 8)], minor=True)
cax.tick_params(axis='y', which='both', direction='in')
cax.set_ylabel('10% Spitzer Heating Rate [erg/s]')
axHist.step(prof.x, 0.1*prof['spitzer_heating_rate'].T.sum(axis=0),
label='10% Spitzer Cond. (heating)', linewidth=2, color='C3')
axHist.step(prof.x, -0.1*prof['spitzer_heating_rate'].T.sum(axis=0),
label='10% Spitzer Cond. (cooling)', linestyle='dotted', linewidth=2, color='C3')
axHist.step(prof.x, prof['xray_luminosity_0.1_100_keV'].T.sum(axis=0),
label='X-ray Cooling', linestyle='dashed', color='C0')
axHist.step(prof.x, 0.01*prof['spitzer_heating_rate'].T.sum(axis=0),
label='1% Spitzer Cond. (heating)', linewidth=1, color='C1')
axHist.step(prof.x, -0.01*prof['spitzer_heating_rate'].T.sum(axis=0),
label='1% Spitzer Cond. (cooling)', linestyle='dotted', linewidth=1, color='C1')
axHist.set_ylim(3E40, 5E44)
axHist.set_yscale('log')
yticks = [1E41, 1E42, 1E43, 1E44]
axHist.set_yticks(yticks)
axHist.set_yticklabels(yticks)
axHist.set_ylabel('heating/cooling rate [erg/s]')
entropy_ticks = [20, 30,40, 50, 60, 70, 80, 90, 100,
110, 120, 130, 140, 150, 160, 170, 180, 190, 200,
210, 220, 230, 240]
entropy_ticklabels = [20, '', 40, '', 60, '', 80, '', 100,
'', '', '', '', 150, '', '', '', '', 200,
'', '', '', '']
axHist.semilogx()
axHist.set_xlim(20, 250)
axHist.tick_params(axis='x', which='both', direction='in')
axHist.tick_params(axis='y', which='both', direction='in')
axHist.set_xticks([100, 200])
axHist.set_xticklabels([])
axHist.set_xticks(entropy_ticks, minor=True)
axHist.set_xticklabels(entropy_ticklabels, minor=True)
axHist.set_xlabel(r'entropy [keV cm$^2$]')
axHist.axhline(0, ls='-', lw=1, color='k')
axHist.legend(frameon=False, loc=2, fontsize=8, ncol=2, columnspacing=0.7)
fig.savefig('heating_rate_100Myr_poloidal.pdf', bbox_inches='tight')
In [ ]:
0.1*prof['spitzer_heating_rate'].T.sum(axis=0)[:41]
In [ ]:
prof.x[39]
In [ ]:
sum(0.1*prof['spitzer_heating_rate'].T.sum(axis=0)[:40])
In [ ]:
sum(prof['xray_luminosity_0.1_100_keV'].T.sum(axis=0)[:41])
In [ ]:
sum(0.1*prof['spitzer_heating_rate'].T.sum(axis=0)[:55])-sum(prof['xray_luminosity_0.1_100_keV'].T.sum(axis=0)[:55])
In [ ]:
extrema = {'entropy': (20, 250), 'spherical_radius': ((10, 'kpc'), (150, 'kpc'))}
logs = {'entropy': False, 'spherical_radius': False}
fields = ['spitzer_heating_rate', 'xray_luminosity_0.1_100_keV', 'cell_mass']
prof_linear = yt.create_profile(sp, ['entropy', 'spherical_radius'], fields, weight_field=None,
logs=logs, extrema=extrema)
In [ ]:
fig = plt.figure(figsize=(4.5,5.5))
axMesh = fig.add_subplot(111)
divider = make_axes_locatable(axMesh)
axHist = divider.append_axes("bottom", size=1.9, pad=0, sharex=axMesh)
cax = divider.append_axes("right", size="4%", pad=0)
im = axMesh.pcolormesh(prof_linear.x, prof_linear.y/3.086E+21, prof_linear['spitzer_heating_rate'].T.v*0.1,
cmap='RdBu_r', norm=colors.SymLogNorm(vmin=-1E43, vmax=1E43, linthresh=1E41))
#axMesh.set_xscale('log')
#axMesh.set_yscale('log')
r_ticks = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150]
r_ticklabels = [10, 20, 30, 40, 50, 60, '', 80, '', 100, '', '', '', '', 150]
axMesh.set_ylim(10,150)
axMesh.set_yticks([100])
axMesh.set_yticklabels([])
axMesh.set_yticks(r_ticks, minor=True)
axMesh.set_yticklabels(r_ticklabels, minor=True)
axMesh.set_ylabel('Radius (kpc)')
axMesh.tick_params(axis='x', which='both', direction='in', labelbottom=False, top='on')
axMesh.tick_params(axis='y', which='both', direction='in')
axMesh.text(23, 115, '%.0f Myr' % ds.current_time.in_units('Myr'))
cb = plt.colorbar(im, cax=cax)
cb.set_ticks([-1E43, -1E42, -1E41, 0, 1E41, 1E42, 1E43])
#cax.set_yticks([np.linspace(-9E43, -2E42, 8)], minor=True)
cax.tick_params(axis='y', which='both', direction='in')
cax.set_ylabel('10% Spitzer Heating Rate (erg/s)')
axHist.step(prof_linear.x, prof_linear['spitzer_heating_rate'].T.sum(axis=0),
label='Spitzer Cond. (heating)', linewidth=1, color='C1')
axHist.step(prof_linear.x, -prof_linear['spitzer_heating_rate'].T.sum(axis=0),
label='Spitzer Cond. (cooling)', linestyle='dotted', linewidth=1, color='C1')
axHist.step(prof_linear.x, prof_linear['xray_luminosity_0.1_100_keV'].T.sum(axis=0),
label='X-ray Cooling', linestyle='dashed', color='C0')
axHist.step(prof_linear.x, 0.1*prof_linear['spitzer_heating_rate'].T.sum(axis=0),
label='10% Spitzer Cond. (heating)', linewidth=2, color='C3')
axHist.step(prof_linear.x, -0.1*prof_linear['spitzer_heating_rate'].T.sum(axis=0),
label='10% Spitzer Cond. (cooling)', linestyle='dotted', linewidth=2, color='C3')
axHist.set_ylim(2E41, 4E45)
axHist.set_yscale('log')
yticks = [1E42, 1E43, 1E44, 1E45]
axHist.set_yticks(yticks)
axHist.set_yticklabels(yticks)
axHist.set_ylabel('Heating/Cooling rate (erg/s)')
entropy_ticks = [20, 30, 40, 50, 60, 70, 80, 90, 100,
110, 120, 130, 140, 150, 160, 170, 180, 190, 200,
210, 220, 230, 240]
entropy_ticklabels = [20, 30, 40, 50, 60, '', 80, '', 100,
'', '', '', '', 150, '', '', '', '', 200,
'', '', '', '']
#axHist.semilogx()
axHist.set_xlim(20, 250)
axHist.tick_params(axis='x', which='both', direction='in')
axHist.tick_params(axis='y', which='both', direction='in')
axHist.set_xticks([100, 200])
axHist.set_xticklabels([])
axHist.set_xticks(entropy_ticks, minor=True)
axHist.set_xticklabels(entropy_ticklabels, minor=True)
axHist.set_xlabel(r'Entropy (cm$^2\ $keV)')
axHist.axhline(0, ls='-', lw=1, color='k')
axHist.legend(frameon=False, loc=2, fontsize=8.5, ncol=2, columnspacing=0.5)
fig.savefig('heating_rate_300Myr_poloidal.pdf', bbox_inches='tight')
In [ ]:
prof.save_as_dataset('100Myr_radius_entropy_profile')
In [ ]:
prof_linear.save_as_dataset('300Myr_radius_entropy_profile_linear')
In [ ]:
plt.step(prof.x, prof['cell_mass'].T.sum(axis=0),
label='cell mass', linewidth=1, color='C1')
plt.semilogx()
In [ ]:
plt.step(prof_linear.x, prof_linear['cell_mass'].T.sum(axis=0),
label='cell mass', linewidth=1, color='C1')
In [ ]:
extrema = {'spherical_radius': ((10, 'kpc'), (150, 'kpc'))}
logs = {'spherical_radius': True}
fields = ['spitzer_heating_rate', 'xray_luminosity_0.1_100_keV', 'cell_mass']
prof_radius = yt.create_profile(sp, 'spherical_radius', fields, weight_field=None,
extrema=extrema, logs=logs)
In [ ]: