In [ ]:
%matplotlib inline
import util
import numpy as np
from matplotlib import pyplot as plt
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import os
import yt
import logging
logging.getLogger('yt').setLevel(logging.ERROR)
from scipy.signal import find_peaks_cwt
from tools import setup_cl
In [ ]:
# Scan for files
#
dirs = ['/home/ychen/data/0only_0525_hinf_nojiggle/',\
'/home/ychen/data/0only_0314_h1_nojiggle/',\
'/home/ychen/data/0only_0330_h0_nojiggle/',\
'/home/ychen/data/0only_0518_hydro_nojiggle/'
]
colors, labels = setup_cl(dirs)
def rescan(dir, printlist=False):
files = util.scan_files(dir, '*hdf5_plt_cnt_[0-9][0-9][0-9][0-9]', walk=True, printlist=printlist)
return files
for dir in dirs:
files = rescan(dir, False)
print(files[-1])
#ds = yt.load(files[21].fullpath)
#print(ds.basename)
In [ ]:
fno = '0260'
profs = {}
#zlim = ((-40, 'kpc'), (40, 'kpc'))
fields = ['velocity_z', 'density', 'pressure', 'cylindrical_radial_velocity', 'magnetic_pressure']
#figdir = os.join(dir, 'profiles/')
#for f in files[50:51]:
for dir in dirs:
ds = yt.load(os.path.join(dir, 'MHD_Jet_nojiggle_hdf5_plt_cnt_%s' % fno))# load data
print(dir, ds, ds.current_time.in_units('Myr'))
disk = ds.disk([0,0,0], [0,0,1], (0.25, 'kpc'), (80, 'kpc'))
profs[dir] = yt.create_profile(disk, 'z', fields, n_bins=400, logs={'z': False}, weight_field='cell_mass')
In [ ]:
for dir, prof in profs.items():
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(prof.x.in_units('kpc'), np.abs(prof['velocity_z'])/3E9, '-', ms=3, label='|vz|/0.1c')
ax.plot(prof.x.in_units('kpc'), prof['cylindrical_radial_velocity']/3E8, '--', ms=3, label='|vr|/0.01c')
p0 = 1.17E-9
ax.plot(prof.x.in_units('kpc'), prof['pressure']/p0, '-', ms=3, label='pressure/(%.2e dyne/cm**2)' % p0)
rho0 = 1.69E-26
ax.plot(prof.x.in_units('kpc'), prof['density']/rho0, '-x', ms=3, label='density/(%.2e g/cm**3)' % rho0)
#ax.plot(prof.x.in_units('kpc')[peaks], prof['density'][peaks]/1E-26, 'x', ms=15, label='peaks')
#plt.plot(zs, vzs, '-', label=label, \
# alpha=float(ds.basename[-4:])/1400)#, c='b')
plt.ylim(-0.5,2.0)
plt.hlines(0, -40, 40, linestyle=':')
plt.xlim(-40, 40)
plt.legend(loc=2, fontsize='small')
peaks = find_peaks_cwt(prof['density'], np.arange(1,10))
print(prof.x.in_units('kpc')[peaks])
plt.vlines(prof.x.in_units('kpc')[peaks], -0.5, 2.0, color='grey', linestyle=':')
plt.title('%s %.2f Myr' % (labels[dir], ds.current_time.in_units('Myr')))
plt.show()
#plt.annotate('%.2f Myr' % ds.current_time.in_units('Myr'), (1,1), xytext=(0.80, 0.9), textcoords='axes fraction')
#plt.savefig(figdir+'vz_p_jet_%iMyr_cell_mass.png' % ds.current_time.in_units('Myr'))
In [ ]:
from scipy.signal import find_peaks_cwt
peaks = find_peaks_cwt(prof['density'], np.arange(1,15))
print(peaks)
print(prof.x.in_units('kpc')[peaks])
In [ ]:
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(prof.x.in_units('kpc'), np.abs(prof['velocity_z'])/3E9, '-x', ms=3, label='|vz|/0.1c')
ax.plot(prof.x.in_units('kpc'), prof['cylindrical_radial_velocity']/3E8, '--', ms=3, label='|vr|/0.01c')
ax.plot(prof.x.in_units('kpc'), prof['pressure']/1E-9, '-^', ms=3, label='pressure/(1E-9 dyne/cm**2)')
ax.plot(prof.x.in_units('kpc'), prof['density']/1E-26, 'o-', ms=3, label='density/(1E-26 g/cm**3)')
#ax.plot(prof.x.in_units('kpc')[peaks], prof['density'][peaks]/1E-26, 'x', ms=15, label='peaks')
#plt.plot(zs, vzs, '-', label=label, \
# alpha=float(ds.basename[-4:])/1400)#, c='b')
plt.ylim(-0.5,1.8)
plt.hlines(0, -40, 40, linestyle=':')
plt.xlim(-40, 40)
plt.legend(loc=2, fontsize='small')
plt.vlines(prof.x.in_units('kpc')[peaks], -0.5, 1.8, color='grey', linestyle=':')
#plt.text(40,1.1, '%.2f Myr\ncell mass weighted' % ds.current_time.in_units('Myr'))
#plt.annotate('%.2f Myr' % ds.current_time.in_units('Myr'), (1,1), xytext=(0.80, 0.9), textcoords='axes fraction')
#plt.savefig(figdir+'vz_p_jet_%iMyr_cell_mass.png' % ds.current_time.in_units('Myr'))
In [ ]:
fno = '0260'
profs = {}
#zlim = ((-40, 'kpc'), (40, 'kpc'))
fields = ['density', 'magnetic_pressure', 'magnetic_field_poloidal', 'magnetic_field_toroidal', 'magnetic_field_z']
#figdir = os.join(dir, 'profiles/')
#for f in files[50:51]:
for dir in dirs:
ds = yt.load(os.path.join(dir, 'data/MHD_Jet_nojiggle_hdf5_plt_cnt_%s' % fno))# load data
print(dir, ds, ds.current_time.in_units('Myr'))
disk = ds.disk([0,0,0], [0,0,1], (0.25, 'kpc'), (100, 'kpc'))
if 'hydro' in dir:
profs[dir] = yt.create_profile(disk, 'z', ['density'], n_bins=400, logs={'z': False}, weight_field='cell_volume')
else:
profs[dir] = yt.create_profile(disk, 'z', fields, n_bins=400, logs={'z': False}, weight_field='cell_volume')
In [ ]:
sim_rhoCore = 9.6E-26 # g/cm3 = N particles/cm3 * 1.67E-24 = 1.91 * N electrons/cm3 * 1.67E-24
sim_densityBeta = 0.53 # for the beta-model density profile
sim_rCore = 8.1E22 # 1.26 arcmin *20.81 kpc/arcmin
kpc = 3.086e+21
distance = np.linspace(0.1*kpc, 150*kpc, 150)
densityBG = sim_rhoCore*(1.0 + (distance/sim_rCore)**2)**(-1.5*sim_densityBeta)
fig, ax = plt.subplots(figsize=(6,4))
for dir, prof in profs.items():
ax.plot(prof.x.in_units('kpc'), np.abs(prof['density']), '-', label=labels[dir], c=colors[dir])
print(dir, np.average(prof['density'][200:320]))
ax.plot(distance/kpc, densityBG, ls=':', label='ICM background')
plt.semilogy()
plt.ylim(1E-28, 1E-25)
plt.ylabel(r'density (g/cm$^3$)')
plt.xlim(0, 45)
plt.xlabel('z (kpc)')
plt.legend(loc=4)
plt.figtext(0.15, 0.80, '%.2f Myr' % (ds.current_time.in_units('Myr')))
plt.savefig('line_profile_nojiggle_260.pdf')
plt.show()
In [ ]:
sim_rhoCore = 9.6E-26 # g/cm3 = N particles/cm3 * 1.67E-24 = 1.91 * N electrons/cm3 * 1.67E-24
sim_densityBeta = 0.53 # for the beta-model density profile
sim_rCore = 8.1E22 # 1.26 arcmin *20.81 kpc/arcmin
kpc = 3.086e+21
fig, ax = plt.subplots(figsize=(6,4))
for dir, prof in profs.items():
densityBG = sim_rhoCore*(1.0 + (prof.x.v/sim_rCore)**2)**(-1.5*sim_densityBeta)
ax.plot(prof.x.in_units('kpc'), np.abs(prof['density']/densityBG), '-', label=labels[dir], c=colors[dir])
print(dir, np.average((prof['density']/densityBG)[200:320]))
plt.semilogy()
#plt.ylim(1E-28, 1E-25)
plt.ylabel(r'density ratio')
plt.xlim(0, 45)
plt.xlabel('z (kpc)')
plt.legend(loc=4)
plt.figtext(0.15, 0.80, '%.2f Myr' % (ds.current_time.in_units('Myr')))
#plt.savefig('line_profile_nojiggle_260.pdf')
plt.show()
In [ ]:
fig, ax = plt.subplots(figsize=(10,6))
for dir, prof in profs.items():
if 'hydro' in dir: continue
ax.plot(prof.x.in_units('kpc'), np.abs(prof['magnetic_pressure']), '-', ms=3, label=labels[dir], c=colors[dir])
plt.hlines(0, -40, 40, linestyle=':')
plt.ylim(1E-13, 1.5E-9)
plt.semilogy()
plt.xlim(0, 45)
plt.legend(loc=2)
plt.title('Magnetic Pressure at %.2f Myr' % (ds.current_time.in_units('Myr')))
plt.show()
In [ ]:
fig, ax = plt.subplots(figsize=(10,6))
for dir, prof in profs.items():
ax.plot(prof.x.in_units('kpc'), np.abs(prof['magnetic_field_toroidal']), '-', ms=3, label=labels[dir], c=colors[dir])
plt.hlines(0, -40, 40, linestyle=':')
plt.ylim(1E-11, 1E-3)
#plt.semilogy()
plt.xlim(-50, 50)
plt.legend(loc=2)
plt.title('Toroidal Magnetic Fields at %.2f Myr' % (ds.current_time.in_units('Myr')))
plt.show()
In [ ]:
fig, ax = plt.subplots(figsize=(10,6))
for dir, prof in profs.items():
ax.plot(prof.x.in_units('kpc'), np.abs(prof['magnetic_field_poloidal']), '-', ms=3,
label=labels[dir], c=colors[dir], alpha=0.7)
plt.hlines(0, -40, 40, linestyle=':')
plt.ylim(1E-11, 1E-3)
plt.semilogy()
plt.xlim(-50, 50)
plt.legend(loc=2)
plt.title('Poloidal Magnetic Fields at %.2f Myr' % (ds.current_time.in_units('Myr')))
plt.show()
In [ ]:
fig, ax = plt.subplots(figsize=(10,6))
for dir, prof in profs.items():
ax.plot(prof.x.in_units('kpc'), np.abs(prof['magnetic_field_z']), '-', ms=3, label=labels[dir], c=colors[dir])
plt.hlines(0, -40, 40, linestyle=':')
plt.ylim(1E-10, 1E-3)
plt.semilogy()
plt.xlim(-50, 50)
plt.legend(loc=2)
plt.title('Bz at %.2f Myr' % (ds.current_time.in_units('Myr')))
plt.show()
In [ ]: