In [ ]:

    

%matplotlib inline
import yt
import logging
logging.getLogger('yt').setLevel(logging.ERROR)
import numpy as np
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt

from tools import setup_cl


    

In [ ]:

    

print(matplotlib.rcParams['figure.dpi'])
dirs = ['/home/ychen/data/0only_0330_h0_nojiggle/',\
        '/home/ychen/data/0only_0518_hydro_nojiggle/',\
        '/home/ychen/data/0only_0314_h1_nojiggle/',\
        '/home/ychen/data/0only_0525_hinf_nojiggle/',\
        ]

colors, labels = setup_cl(dirs)

tt = np.linspace(0, 5, 10)

for dir in dirs:
    data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
    label = labels[dir]
    plt.fill_between(data[:,1], data[:,2], -data[:,3],\
                     facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.7, label=label)

    #plt.fill_between(data[:,1], data[:,4], -data[:,5],\
    #                 facecolor=colors[dir], lw=1, linestyle=':', edgecolor=colors[dir], alpha=0.3)
plt.plot(tt, 10*tt, ':', c='orange', label='10 kpc / Myr')
plt.plot(tt, 8*tt, '--', c='cyan', label='8 kpc / Myr')
plt.legend(loc=2)

handles, labels = plt.axes().get_legend_handles_labels()
# reverse the order
plt.legend(handles[::-1], labels[::-1])
plt.grid(ls='--', alpha=0.5)
plt.ylabel('z (kpc)')
plt.xlabel('t (Myr)')
plt.xlim(0,5)
plt.ylim(0,60)
plt.savefig('LobeSize_nojiggle.pdf')
#plt.show()
#plt.semilogx()


    

In [ ]:

    




    

In [ ]:

    

def cal_vel(locations, times):
    loc_c = (locations[:-1] + locations[1:])/2
    vel = (locations[:-1] - locations[1:])/(times[:-1]-times[1:])
    return loc_c, vel


    

In [ ]:

    

for zEdge in [zEdgep, zEdgem]:
    zEdge_c, vEdge = cal_vel(zEdge, time)
    plt.plot(np.abs(zEdge_c), np.abs(vEdge))
plt.xlabel('z (kpc)')
plt.ylabel('velocity (kpc/Myr)')


    

In [ ]:

    

sim_Tcore		= 3.48E7
sim_Tout		= 7.42E7
sim_rhoCore		= 9.6E-26		# g/cm3 = N particles/cm3 * 1.67E-24 = 1.91 * N electrons/cm3 * 1.67E-24
sim_mu			= 0.61			# mean molecular weight
sim_windVel     = 0.0
sim_gammaICM	= 1.66666
sim_bzAmbient	= 0.0 #1.8623E-6
sim_densityProfile = "betacore"
sim_densityBeta = 0.53			# for the beta-model density profile
sim_rCore		= 8.1E22		# 1.26 arcmin *20.81 kpc/arcmin
sim_rCoreT		= 1.85E23		# 60kpc

kpc = 3.086e+21
kB = 1.38E-16
R = 8.3144725E7

distance = np.linspace(0.1*kpc, 60*kpc, 150)

densityBG = sim_rhoCore*(1.0 + (distance/sim_rCore)**2)**(-1.5*sim_densityBeta)
tempBG = sim_Tout*(1.0+(distance/sim_rCoreT)**3)/(sim_Tout/sim_Tcore+(distance/sim_rCoreT)**3)
pressureBG = densityBG*tempBG*R

Myr = 3.155760e+13

sim_velJet = 3.0E9
rhoJet = 1.27E-26


    

In [ ]:

    




    

In [ ]:

    

eta = rhoJet*2/densityBG
velHead = np.sqrt(eta)/(1+np.sqrt(eta))*sim_velJet
plt.plot(distance/kpc, velHead/kpc*Myr, ls='--', c='grey', label=r'$2\rho_j$')

eta = rhoJet/densityBG
velHead = np.sqrt(eta)/(1+np.sqrt(eta))*sim_velJet
plt.plot(distance/kpc, velHead/kpc*Myr, ls='-', c='grey', label=r'$\rho_j$')

eta = rhoJet*0.5/densityBG
velHead = np.sqrt(eta)/(1+np.sqrt(eta))*sim_velJet
plt.plot(distance/kpc, velHead/kpc*Myr, ls='-.', c='grey', label=r'$\frac{1}{2}\rho_j$')

eta = rhoJet/4/densityBG
velHead = np.sqrt(eta)/(1+np.sqrt(eta))*sim_velJet
plt.plot(distance/kpc, velHead/kpc*Myr, ls=':', c='grey', label=r'$\frac{1}{4}\rho_j$')


for dir in dirs:
    data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
    fnumber = data[:,0]
    time = data[:,1]
    zEdgep = data[:,2]
    zEdgem = data[:,3]
    zCenterp = data[:,4]
    zCenterm = data[:,5]
    rEdgep = data[:,6]
    rEdgem = data[:,7]

    zEdge_c, vEdge = cal_vel(zEdgem, time)
    plt.plot(np.abs(zEdge_c), np.abs(vEdge), color=colors[dir], label=labels[dir])

    #zEdge_c, vEdge = cal_vel(zEdgem, time)
    #plt.plot(np.abs(zEdge_c), np.abs(vEdge), color=colors[dir])

#plt.xlim(0, 60)
plt.xlabel('distance (kpc)')
plt.ylim(0, 16)
plt.ylabel('Jet head velocity (kpc/Myr)')
plt.legend(loc=4)


    

In [ ]: