In [ ]:
import yt
import matplotlib.pyplot as plt
import numpy as np
from synchrotron.yt_synchrotron_emissivity import setup_part_file, add_synchrotron_dtau_emissivity
In [ ]:
ds = yt.load('/d/d5/ychen/2015_production_runs/1022_h1_10Myr/data/MHD_Jet_10Myr_hdf5_plt_cnt_0600')
setup_part_file(ds)
In [ ]:
ptype = 'lobe'
ds.add_particle_filter(ptype)
ad = ds.all_data()
In [ ]:
pdata = ad[ptype, 'particle_tau']
null = plt.hist(np.log10(pdata), range=(-6, -2), bins='auto')
plt.title('log particle_tau')
In [ ]:
pdata = ad[ptype, 'particle_gamc']
null = plt.hist(np.log10(pdata), range=(1, 6), bins='auto')
plt.title('log particle_gamc')
In [ ]:
pdata = ad[ptype, 'particle_dtau']
null = plt.hist(np.log10(pdata), range=(-6, -2), bins='auto')
median = np.median(np.log10(pdata))
plt.vlines(median, 0, 3000)
plt.title('log particle_dtau')
In [ ]:
pdata = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1'])**(1./3.)/ ad[ptype, 'particle_dtau']
null = plt.hist(np.log10(pdata), range=(1, 6), bins='auto', label='gamc_dtau (den1)', alpha=0.8)
pdata = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den0'])**(1./3.)/ ad[ptype, 'particle_dtau']
null = plt.hist(np.log10(pdata), range=(1, 6), bins='auto', label='gamc_dtau', alpha=0.8)
pdata = ad[ptype, 'particle_gamc']
null = plt.hist(np.log10(pdata), range=(1, 6), bins='auto', label='gamc', alpha=0.8)
plt.legend()
plt.title('log particle_gamc_dtau')
In [ ]:
pdata = (ad[ptype, 'particle_dens'])
null = plt.hist(np.log10(pdata), bins='auto', label=r'$\rho$', alpha=0.8)
pdata = (ad[ptype, 'particle_den1'])
null = plt.hist(np.log10(pdata), bins='auto', label=r'$\rho_1$', alpha=0.8)
plt.xlim(-28, -24)
#plt.vlines(0, 0, 4000)
plt.legend()
plt.title(r'log $\rho$')
In [ ]:
pdata = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1'])
null = plt.hist(np.log10(pdata), bins='auto', label=r'log $\rho / \rho_1$', alpha=0.8)
pdata = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den0'])
null = plt.hist(np.log10(pdata), bins='auto', label=r'log $\rho / \rho_0$', alpha=0.8)
plt.legend()
#plt.xlim(-3, 2)
plt.axvline(0)
plt.title(r'log $\rho / \rho_1$')
In [ ]:
add_synchrotron_dtau_emissivity(ds, nu=(150, 'MHz'))
In [ ]:
pdata = ad[ptype, 'particle_sync_spec_150.0MHz']
null = plt.hist(np.log10(pdata), range=(-60, -35), bins='auto')
In [ ]:
gamc_dtau = (ad[ptype, 'particle_dens']/ad[ptype, 'particle_den0'])**(1./3.)/ ad[ptype, 'particle_dtau']
B2 = (ad[ptype, 'particle_magx']**2 + ad[ptype, 'particle_magy']*2 + ad[ptype, 'particle_magz']**2)/2
In [ ]:
plt.scatter(np.log10(gamc_dtau), np.log10(B2), s=1, linewidth=0)
plt.xlim(1,6)
In [ ]:
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(gamc_dtau), c=np.log10(ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1']), vmin=-1, vmax=1, s=1, cmap='seismic', linewidth=0, alpha=0.5)
plt.ylim(1,7)
plt.xlim(-29, -24)
plt.ylabel(r'$\log \gamma_c$')
plt.xlabel(r'$\log \rho_1$')
cb = plt.colorbar()
cb.set_label(r'$\rho / \rho_1$')
In [ ]:
import matplotlib
plt.scatter(np.log10(ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1']), np.log10(gamc_dtau), \
c=np.log10(ad[ptype, 'particle_den1']), vmin=-28, vmax=-25, s=1, linewidth=0, alpha=0.5)
plt.ylim(1,7)
plt.xlim(-3, 2)
plt.ylabel(r'$\log \gamma_c$')
plt.xlabel(r'$\log \rho / \rho_1$')
cb = plt.colorbar()
cb.set_label(r'$\rho_1$')
In [ ]:
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(gamc_dtau*ad[ptype, 'particle_den1']**(1/3)), \
s=1, linewidth=0, alpha=0.5)
plt.ylim(-7,-2)
plt.twinx()
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(gamc_dtau), \
s=1, linewidth=0, alpha=0.5, c='r')
plt.ylim(1,6)
plt.xlim(-28, -24)
#plt.ylabel(r'$\log \gamma_c$')
#plt.xlabel(r'$\log \rho / \rho_1$')
#cb = plt.colorbar()
#cb.set_label(r'$\rho_1$')
In [ ]:
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(gamc_dtau*ad[ptype, 'particle_den1']**(1/3)), \
s=1, linewidth=0, alpha=0.5)
plt.ylim(-7,-2)
plt.twinx()
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(gamc_dtau), \
s=1, linewidth=0, alpha=0.5, c='r')
plt.ylim(1,6)
plt.xlim(-28, -24)
#plt.ylabel(r'$\log \gamma_c$')
#plt.xlabel(r'$\log \rho / \rho_1$')
#cb = plt.colorbar()
#cb.set_label(r'$\rho_1$')
In [ ]:
import matplotlib
plt.scatter(np.log10(ad[ptype, 'particle_den1']),
np.log10(ad[ptype, 'particle_dens']/ad[ptype, 'particle_den1']),
vmin=-28, vmax=-25, s=1, linewidth=0, alpha=0.5)
plt.xlim(-28,-24)
plt.ylim(-3, 2)
plt.xlabel(r'$\log \rho_1$')
plt.ylabel(r'$\log \rho / \rho_1$')
In [ ]:
import matplotlib
Myr = yt.units.Myr.in_units('s').v
plt.scatter(ad[ptype, 'particle_tadd']/Myr,
np.log10(ad[ptype, 'particle_den1']),
s=1, linewidth=0, alpha=0.5)
plt.xlim(0, 10)
plt.ylim(-28,-24)
plt.xlabel('time (Myr)')
plt.ylabel(r'$\log \rho_1$')
In [ ]:
null = plt.hist(np.log10(ad[ptype, 'particle_den1']), bins='auto')
In [ ]:
e = 4.803E-10 #esu
me = 9.109E-28 #g
c = 2.998E10 #cm/s
nuc = 3*gamc_dtau**2*e*np.sqrt(B2)/(4*np.pi*me*c)
plt.scatter(np.log10(ad[ptype, 'particle_den1']), np.log10(nuc), c=np.log10(np.sqrt(B2)), vmin=-4, vmax=-3, s=1, linewidth=0)
plt.colorbar()
plt.ylim(6,15)
In [ ]:
plt.scatter(np.log10(ad['gas', 'density']), np.log10(ad['gas', 'magnetic_field_strength']), s=1, linewidth=0)
In [ ]:
plt.scatter(np.log10(ad['gas', 'density']), np.log10(ad['gas', 'magnetic_field_strength']), s=1, linewidth=0)
plt.ylim(-8,-3)
In [ ]: