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

#@yt.derived_field('xray_')
#def _xray_exclude_jet(field, data):


    

In [ ]:

    

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']
zoom_fac = 6
metallicity = 0.5
data = {}
for j, dir in enumerate(dirs[1:2]):
    times = np.zeros(len(filenumbers))
    xray_lum = 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))
        times[i] = ds.current_time.in_units('Myr')
        xray_fields = yt.add_xray_emissivity_field(ds, 0.1, 100, table_type='apec', metallicity=metallicity)
        emis, l, pemis = xray_fields
        box = ds.box(ds.domain_left_edge/zoom_fac, ds.domain_right_edge/zoom_fac)
        print(dir, filenumber, ds.current_time.in_units('Myr'))
        xray_lum[i] = box.quantities.total_quantity(l)


    

In [ ]:

    

times = np.array([1.584457825611463, 3.1689045195099665, 9.506526296590856, 51.82957048885901, 99.36638020798718])
plt.plot(times, xray_lum, 'o-')
plt.xlabel('time (Myr)')
plt.ylabel('X-ray luminosity (erg/s)')


    

In [ ]:

    

metallicity = 0.5
ds = yt.load('/home/ychen/data/0only_0529_h1/data/MHD_Jet_hdf5_plt_cnt_0000')
xray_fields = yt.add_xray_emissivity_field(ds, 0.1, 100, table_type='apec', metallicity=metallicity)
sp = ds.sphere([0,0,0], (100, 'kpc'))
emis, l, pemis = xray_fields


    

In [ ]:

    

mask = sp['jet '] < 0.1


    

In [ ]:

    

kB = yt.physical_constants.kboltz
t_cool = sp['pressure']*sp['cell_volume']/sp[l]


    

In [ ]:

    

t_cool.in_units('Myr')


    

In [ ]:

    

plt.scatter(t_cool.in_units('Myr'), sp['spherical_radius'].in_units('kpc'), c=sp['jet '], s=1, lw=0)


    

In [ ]:

    

plot = yt.SlicePlot(ds, 'x', emis)
plot.zoom(6)
plot.show()


    

In [ ]:

    

ds = yt.load(os.path.join(dirs[1], 'data/MHD_Jet_10Myr_hdf5_plt_cnt_0600'))
xray_fields = yt.add_xray_emissivity_field(ds, 0.1, 100.0, table_type='apec', metallicity=metallicity)
plot = yt.ProjectionPlot(ds, 'x', xray_fields[0])
plot.zoom(6)
plot.show()


    

In [ ]:

    

plot = yt.ProjectionPlot(ds, 'z', emis)
plot.zoom(6)
plot.show()


    

In [ ]:

    

ds = yt.load(os.path.join(dirs[1], 'data/MHD_Jet_10Myr_hdf5_plt_cnt_0600'))
xray_fields = yt.add_xray_emissivity_field(ds, 0.1, 100, table_type='apec', metallicity=metallicity)
emis, l, pemis = xray_fields


    

In [ ]:

    

yt.ProfilePlot(ds.all_data(), 'temperature', l, weight_field=None)


    

In [ ]:

    

yt.ProfilePlot(ds.all_data(), 'jet ', l, weight_field=None, accumulation=True, x_log=False, y_log={l:False})


    

In [ ]:

    

yt.PhasePlot(ds.all_data(), 'spherical_r', l, 'temperature')


    

In [ ]:

    

import h5py
h5f = h5py.File('apec_emissivity_v2.h5', 'r')


    

In [ ]:

    

for k, v in h5f.items():
    print(k, v)


    

In [ ]:

    

E = h5f['E'].value
dE = np.diff(E)
emissivity_primordial = h5f['emissivity_primordial'].value
emissivity_metals =h5f['emissivity_metals'].value
log_T = h5f['log_T'].value
h5f.close()


    

In [ ]:

    

E


    

In [ ]:

    

mask = E[1:] > 0.1
plt.plot(log_T, np.sum(emissivity_primordial[:,mask]*dE[mask], axis=-1), label='primordial')
plt.plot(log_T, np.sum(emissivity_metals[:,mask]*dE[mask], axis=-1), label='metals')
plt.plot(log_T, np.sum((emissivity_primordial[:,mask]+emissivity_metals[:,mask])*dE[mask], axis=-1), label='total (Z=1)')
#plt.plot(log_T, np.sum((emissivity_primordial.value+0.5*emissivity_metals.value)*ebin, axis=-1), label='total (Z=0.5)')
plt.legend()
plt.xlabel('log T')
plt.ylabel(r'x-ray (0.1-100 keV) cooling rate (erg cm$^3$/s)')
plt.semilogy()
plt.ylim(1E-25, 1E-22)


    

In [ ]:

    

import h5py
h5fc = h5py.File('cloudy_emissivity_v2.h5', 'r')


    

In [ ]:

    

for k, v in h5fc.items():
    print(k, v)


    

In [ ]:

    

E = h5fc['E']
ebin = np.diff(E)
emissivity_primordial = h5fc['emissivity_primordial'][-1,:,:]
emissivity_metals =h5fc['emissivity_metals'][-1,:,:]
log_T = h5fc['log_T']
log_nH = h5fc['log_nH']


    

In [ ]:

    

plt.plot(log_T, np.sum(emissivity_primordial*ebin, axis=-1), label='primordial')
plt.plot(log_T, np.sum(emissivity_metals*ebin, axis=-1), label='metals')
plt.plot(log_T, np.sum((emissivity_primordial+emissivity_metals)*ebin, axis=-1), label='total (Z=1)')
#plt.plot(log_T, np.sum((emissivity_primordial.value+0.5*emissivity_metals.value)*ebin, axis=-1), label='total (Z=0.5)')
plt.legend()
plt.xlabel('log T')
plt.semilogy()
plt.ylim(1E-25, 1E-22)


    

In [ ]:

    

log_nH.value


    

In [ ]: