In [ ]:

    

import yt
import logging
import numpy as np
import matplotlib
matplotlib.rcParams['savefig.dpi'] = 150
import matplotlib.pyplot as plt
import multiprocessing
logging.getLogger('yt').setLevel(logging.ERROR)
dss = yt.load('/home//ychen/d9/FLASH4/2016_test/1028_metal_test/MHD_Jet_hdf5_plt_cnt*')


    

In [ ]:

    

end=100
ds = dss[end]
ad = ds.all_data()
print ad.center
print ad.right_edge
print ad.left_edge
print ds.field_list
#for field in ds.derived_field_list:
#    if True:
#        print field


    

In [ ]:

    

center = [0.0,0.0,0.0]
radius = (1.5E23, 'code_length')
sp = ds.sphere(center, radius)
#print sqrt(sp['velocity_x']**2+sp['velocity_y']**2+sp['velocity_z']**2)
#print sp['velocity_magnitude']
plt.hist(np.log2(sp['cell_volume']/sp['cell_volume'].min()), bins=8)


    

In [ ]:

    

#print sp['thermal_energy']
print sp['kinetic_energy']/sp['density']
print
#print sp['thermal_energy'] + sp['kinetic_energy']/sp['density']
#print sp['total_energy']
print sp['pressure']
print sp['pressure']/sp['density']/(1.6666-1)


    

In [ ]:

    

print sum(sp['thermal_energy']*sp['density']*sp['cell_volume']) + sum(sp['kinetic_energy']*sp['cell_volume'])
print sum(sp['thermal_energy']*sp['cell_mass']) + sum(sp['kinetic_energy']*sp['cell_volume'])
print sum(sp['thermal_energy']*sp['cell_mass']) + sum(sp['kinetic_energy']/sp['density']*sp['cell_mass'])


    

In [ ]:

    

Eth = sum(sp['pressure']/sp['density']/(1.6666-1)*sp['cell_mass'])
Ekin = sum(sp['kinetic_energy']/sp['density']*sp['cell_mass'])

Etot = Eth + Ekin
print Etot


    

In [ ]:

    

plot = yt.SlicePlot(dss[120], 'x', 'jet ', center=(0,0,0))
plot.annotate_sphere(center, radius)
plot.zoom(8)
plot.show()


    

In [ ]:

    

ts = []
Ekins, Eths, Etots, Emags = [], [], [], []


    

In [ ]:

    

for ds in dss[0:13]:
    print ds
    sp = ds.sphere(center, radius)
    
    
    Eth = sum(sp['pressure']/(1.6666-1)*sp['cell_volume']*sp['jet '])
    Ekin = sum(sp['kinetic_energy']*sp['cell_volume']*sp['jet '])
    Emag = sum(sp['magnetic_pressure']*sp['cell_volume']*sp['jet '])
    Etot = Eth + Ekin
    t = ds.current_time
    ts.append(t)
    Eths.append(Eth)
    Ekins.append(Ekin)
    Etots.append(Etot)
    Emags.append(Emag)


    

In [ ]:

    

print ts


    

In [ ]:

    

plt.plot(np.array(ts)/3.154E13, (np.array(Etots)-np.array(Etots[0]))/np.array(Emags), label='total E')
plt.plot(np.array(ts)/3.154E13, (np.array(Ekins)-np.array(Ekins[0]))/np.array(Emags), label='kinetic E')
plt.plot(np.array(ts)/3.154E13, (np.array(Eths)-np.array(Eths[0]))/np.array(Emags), label='thermal E')
#plt.plot(ts/3.154E13, (np.array(Emags)-np.array(Emags[0]))*30, label='magnetic E')
plt.legend(loc=0)
plt.xlim(0,25)
plt.xlabel('t (Myr)')
plt.ylabel('E/Emag')
#plt.title('h1')


    

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 [ ]: