In [1]:
%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

In [2]:
from yt_conduction_fields import *

In [3]:
ds = yt.load('/home/ychen/data/0only_1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_1380')
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')


Out[3]:
299.95294319868765 Myr

In [4]:
plot = yt.SlicePlot(ds, 'y', 'jet ')
plot.zoom(4)
plot.set_zlim('jet ', 1E-7, 1)
plot.annotate_contour('jet ', ncont=2, clim=(1E-4, 1E-3), take_log=True)
plot.show()




In [5]:
plot = yt.SlicePlot(ds, 'y', 'xray_emissivity_0.1_100_keV')
plot.zoom(4)
plot.show()




In [6]:
plot = yt.SlicePlot(ds, 'y', 'xray_cooling_time', )
plot.set_unit('xray_cooling_time', 'Gyr')
plot.zoom(4)
plot.show()




In [16]:
plot = yt.SlicePlot(ds, 'y', 'temperature')
plot.set_unit('temperature', 'keV', equivalency='thermal')
plot.zoom(4)
plot.show()




In [8]:
plot = yt.SlicePlot(ds, 'y', 'entropy', data_source=sp)
plot.set_zlim('entropy', 30, 250)
plot.zoom(4)
plot.show()




In [13]:
plot = yt.SlicePlot(ds, 'y', 'H_nuclei_density')
plot.zoom(4)
plot.show()




In [10]:
plot = yt.SlicePlot(ds, 'y', 'temperature_gradient_magnitude', data_source=sp, width=(240, 'kpc'))
plot.set_unit('temperature_gradient_magnitude', 'K/pc')
plot.show()




In [11]:
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.show()


yt : [WARNING  ] 2018-05-16 19:40:57,740 Plot image for field ('gas', 'temperature_gradient_x') has both positive and negative values. Min = -0.000000, Max = 0.000000.
yt : [WARNING  ] 2018-05-16 19:40:57,745 Switching to symlog colorbar scaling unless linear scaling is specified later


In [17]:
plot = yt.SlicePlot(ds, 'y', 'spitzer_conduction_coefficient', data_source=sp, width=(240, 'kpc'))
plot.set_zlim('spitzer_conduction_coefficient', 1E29, 1E32)
plot.show()




In [6]:
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.show()


yt : [WARNING  ] 2018-06-20 12:36:02,876 Plot image for field ('gas', 'spitzer_heat_flux_x') has both positive and negative values. Min = -0.059746, Max = 0.048133.
yt : [WARNING  ] 2018-06-20 12:36:02,880 Switching to symlog colorbar scaling unless linear scaling is specified later


In [7]:
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.show()


yt : [WARNING  ] 2018-06-20 12:38:01,457 Plot image for field ('gas', 'spitzer_heat_flux_divergence') has both positive and negative values. Min = -0.000000, Max = 0.000000.
yt : [WARNING  ] 2018-06-20 12:38:01,462 Switching to symlog colorbar scaling unless linear scaling is specified later


In [4]:
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.show()


yt : [WARNING  ] 2018-08-05 14:27:38,291 Plot image for field ('gas', 'spitzer_heating_rate') has both positive and negative values. Min = -125325531416275362071510391439947378720768.000000, Max = 190128496943913322050048100259300504502272.000000.
yt : [WARNING  ] 2018-08-05 14:27:38,296 Switching to symlog colorbar scaling unless linear scaling is specified later


In [5]:
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.show()


yt : [WARNING  ] 2018-08-05 14:29:36,532 Plot image for field ('gas', 'total_heating_rate') has both positive and negative values. Min = -125336774535115103741047678339331320184832.000000, Max = 190110057323424332180498954269088278380544.000000.
yt : [WARNING  ] 2018-08-05 14:29:36,536 Switching to symlog colorbar scaling unless linear scaling is specified later


In [28]:
plot = yt.SlicePlot(ds, 'y', 'total_cooling_time', data_source=sp, width=(240, 'kpc'))

plot.set_log('total_cooling_time' , False)
plot.set_cmap('total_cooling_time', 'seismic')
plot.set_zlim('total_cooling_time', -1E2, 1E2)
plot.show()




In [6]:
extrema = {'entropy': (20, 250)}
logs = {'entropy': True}
fields = ['spitzer_heating_rate', 'xray_luminosity_0.1_100_keV', 'cell_mass']
prof_entropy = yt.create_profile(sp, 'entropy', fields, weight_field=None, 
                         extrema=extrema, logs=logs)

In [89]:
fig = plt.figure(figsize=(6,4))
plt.step(prof_entropy.x, -prof_entropy['xray_luminosity_0.1_100_keV'],
         label='X-ray Cooling', linestyle='dotted')
plt.step(prof_entropy.x, prof_entropy['spitzer_heating_rate'],
         label='Spitzer Thermal Conduction', linestyle='dashed')
plt.step(prof_entropy.x, prof_entropy['spitzer_heating_rate']-prof_entropy['xray_luminosity_0.1_100_keV'], 
         label='Total $f_{\mathrm{Sp}}$ = 1', linewidth=1)
plt.step(prof_entropy.x, 0.1*prof_entropy['spitzer_heating_rate']-prof_entropy['xray_luminosity_0.1_100_keV'], 
         label=r'Total $f_{\mathrm{Sp}}$ = 0.1', linewidth=2)
plt.ylim(-1E45, 1E45)
plt.yscale('symlog', linthreshy=1E42)
yticks = [-1E45, -1E44, -1E43, -1E42, 0, 1E42, 1E43, 1E44, 1E45]

plt.yticks(yticks, yticks)
plt.ylabel('Heating(+)/Cooling(-) rate (erg/s)')
plt.semilogx()
plt.xlim(20, 250)
entropy_ticks = [20, 30, 40, 60, 80, 100, 150, 200]
plt.xticks(entropy_ticks, entropy_ticks)
plt.xlabel(r'Entropy (cm$^2\ $ keV)')
plt.axhline(0, ls='-', lw=1, color='k')
plt.text(170, 1E44, '%.0f Myr' % ds.current_time.in_units('Myr'))
plt.legend(loc=3, fontsize=9, ncol=2, columnspacing=0.5)


Out[89]:
<matplotlib.legend.Legend at 0x7fef01622358>

In [8]:
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 [81]:
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)
plot2 = pp['spitzer_heating_rate']

r_ticks = [30, 40, 50, 60, 80, 100, 150]
pp.set_ylim(20, 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


Out[81]:

In [52]:
prof.field_data


Out[52]:
{('gas',
  'cell_mass'): YTArray([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        ..., 
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           1.19024069e+43,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           1.92838119e+44,   1.13162302e+43,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   2.00533186e+44,   1.89781368e+43]]) g,
 ('gas',
  'spitzer_heating_rate'): YTArray([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        ..., 
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
          -1.04430214e+42,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
          -1.77563490e+43,  -1.78034445e+42,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,  -2.27244710e+43,  -1.81766184e+42]]) erg/s,
 ('gas',
  'xray_luminosity_0.1_100_keV'): YTArray([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        ..., 
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           4.07876239e+41,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           6.46341409e+42,   3.71263215e+41,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   6.42142093e+42,   5.93276637e+41]]) erg/s}

In [45]:
plot2


Out[45]:

In [16]:
from yt.visualization.plot_container import FieldTransform
from yt.visualization.tick_locators import LogLocator

linthresh = 1E43

symlog_transform = FieldTransform('symlog', None, LogLocator()) 
pp._field_transform['spitzer_heating_rate'] = symlog_transform
pp._field_transform['spitzer_heating_rate'].func = linthresh
pp._setup_plots()
pp.annotate_text(12, 200, '%.0f Myr' % ds.current_time.in_units('Myr'))
print(pp._field_transform['spitzer_heating_rate'])
pp.set_zlim('spitzer_heating_rate', -5E44, 5E44)
plot2 = pp['spitzer_heating_rate']
entropy_ticks = [20, 30, 40, 50, 60, 70, 80, 90, 100, 200]
plot2.axes.set_yticks(entropy_ticks)
plot2.axes.set_yticklabels(entropy_ticks)
r_ticks = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
plot2.axes.set_xticks(r_ticks)
plot2.axes.set_xticklabels(r_ticks)
plot2


<yt.visualization.plot_container.FieldTransform object at 0x7f869b026080>
/home/ychen/src/yt-git/yt/visualization/plot_container.py:118: RuntimeWarning: invalid value encountered in log10
  expA = np.floor(np.log10(vmin))
Out[16]:

In [70]:
keV = yt.units.keV
cm = yt.units.cm
kappa = 1*(keV**(5/2)*cm**3)
kappa.convert_to_base()


Out[70]:
1.0274860530720377e-22 cm**8*g**(5/2)/s**5

In [35]:
heating_rate = prof['spitzer_heating_rate'].flatten()
null = plt.hist(np.log10(heating_rate[heating_rate>0]), bins=20, histtype='step', label='heating')
null = plt.hist(np.log10(-heating_rate[heating_rate<0]), bins=20, histtype='step', label='cooling')
plt.legend()


Out[35]:
<matplotlib.legend.Legend at 0x7fbead792c88>

In [ ]: