In [ ]:

    

%matplotlib inline
#%pdb
import numpy as np
import matplotlib.pyplot as plt
import yt
from yt.analysis_modules.spectral_integrator.api import \
     add_xray_emissivity_field

dir ='/home/ychen/d9/FLASH4/stampede/0529_L45_M10_b1_h1/'
fname = dir+'MHD_Jet_hdf5_plt_cnt_0620'

energy_ranges = [(0.3, 1.5), (1.5, 3.5), (3.5, 7.0)]

ds = yt.load(fname)

for energy_range in energy_ranges:
    print add_xray_emissivity_field(ds, *energy_range, with_metals=False, constant_metallicity=0.02)


    

In [ ]:

    

fnames = ['xray_emissivity_%.1f_%.1f_keV' % er for er in energy_ranges]
proj = ds.proj(['xray_emissivity_%.1f_%.1f_keV' % er for er in energy_ranges], 'y')


    

In [ ]:

    

ext = ds.arr([-1.544E23, 1.544E23, -7.72E22, 7.72E22], input_units='code_length')
frb = proj.to_frb(ext[1]-ext[0], (1024,512), height=(ext[3]-ext[2]))
#image = np.array([frb[fname]/frb[fname].max() for fname in fnames]).transpose([1,2,0])
image = np.array([frb[fname] for fname in fnames]).transpose([1,2,0])

print image.shape
print image[:,:,0].max()
print image[:,:,1].max()
print image[:,:,2].max()
image = image/image.max()


    

In [ ]:

    

import matplotlib.pyplot as plt
plt.figure(figsize = (20,10))
plt.imshow(image, origin='lower', extent=ext.in_units('kpc'))
plt.xlabel('z (kpc)')
plt.ylabel('x (kpc)')

plt.show()


    

In [ ]:

    

fnames = ['xray_emissivity_%.1f_%.1f_keV' % er for er in energy_ranges]
proj = ds.proj(['xray_emissivity_%.1f_%.1f_keV' % er for er in energy_ranges], 'z')


    

In [ ]:

    

ext = ds.arr([-7.72E22, 7.72E22, -7.72E22, 7.72E22], input_units='code_length')
frb = proj.to_frb(ext[1]-ext[0], (512,512), height=(ext[3]-ext[2]))
#image = np.array([frb[fname]/frb[fname].max() for fname in fnames]).transpose([1,2,0])
image = np.array([frb[fname] for fname in fnames]).transpose([1,2,0])
print image.shape
image = image/image.max()

plt.figure(figsize = (10,10))
plt.imshow(image, origin='lower', extent=ext.in_units('kpc'))
plt.xlabel('z (kpc)')
plt.ylabel('x (kpc)')
plt.show()


    

In [ ]:

    

oaprojplot = yt.OffAxisProjectionPlot(ds, (0.0, 0.5, 1.0), 'xray_emissivity_0.3_1.5_keV', max_level=4)


    

In [ ]:

    

oaprojplot.zoom(10)
oaprojplot.show()


    

In [ ]:

    

from yt.visualization.plot_window import get_oblique_window_parameters, OffAxisProjectionDummyDataSource, PWViewerMPL
from yt.visualization.fixed_resolution import OffAxisProjectionFixedResolutionBuffer
normal = (0.0, 0.5, 1.0)
center = (0.0, 0.0, 0.0)
width = (40, 'kpc')
depth = (40, 'kpc')
max_level=4
fields = ['xray_emissivity_%.1f_%.1f_keV' % er for er in energy_ranges]

(bounds, center_rot) = get_oblique_window_parameters(normal,center,width,ds,depth=depth)


    

In [ ]:

    

axes_unit=None; weight_field=None
north_vector=None; volume=None; no_ghost=False,
le=None; re=None; interpolated=False; method="integrate"
oap_width = ds.arr((bounds[1] - bounds[0], bounds[3] - bounds[2], bounds[5] - bounds[4]))
OffAxisProj = OffAxisProjectionDummyDataSource(center_rot, ds, normal, oap_width, fields, interpolated, \
                                               weight=weight_field,  volume=volume, no_ghost=no_ghost, \
                                               le=le, re=re, north_vector=north_vector, method=method)
OffAxisProj.dd.max_level = max_level


    

In [ ]:

    

OffAxisProj


    

In [ ]:

    

#%pdb
class OffAxisProjection(PWViewerMPL):
    _plot_type = 'OffAxisSlice'
    _frb_generator = OffAxisProjectionFixedResolutionBuffer

plot = OffAxisProjection(OffAxisProj, bounds, fields=fields, origin='center-window',\
                    periodic=False, oblique=True, fontsize=18)


    

In [ ]:

    

p = plot['xray_emissivity_0.3_1.5_keV']


    

In [ ]:

    

r = plot['xray_emissivity_0.3_1.5_keV'].image._A
g = plot['xray_emissivity_1.5_3.5_keV'].image._A
b = plot['xray_emissivity_3.5_7.0_keV'].image._A
image = np.array([r,g,b]).transpose([1,2,0])
image = image/image.max()


    

In [ ]:

    

print image[:,:,0].max()
print image[:,:,1].max()
print image[:,:,2].max()


    

In [ ]:

    

plt.figure(figsize = (10,10))
plt.imshow(image, origin='lower')
#plt.xlabel('z (kpc)')
#plt.ylabel('x (kpc)')
plt.show()


    

In [ ]: