In [ ]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 150


    

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
mH = 1.6737352238051868e-24 #g

distance = np.linspace(0.1*kpc, 150*kpc, 150)
e = 1.602176562e-09

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/sim_mu
entropyBG = pressureBG*((densityBG/mH*sim_mu)**(-sim_gammaICM))/e


    

In [ ]:

    

keV = 1.1604525E7
rho = 1E-26
T = 1
P = 1E-10
entropy = 1000

plt.plot(distance/kpc, densityBG/rho, label='density / $(10^{-26} g/cm^3)$')
plt.plot(distance/kpc, tempBG/keV/T, ls='--', label='temperature / $keV$')
plt.plot(distance/kpc, pressureBG/P, ls=':', label='pressure / $(10^{-10} dyn / cm^2)$')
plt.plot(distance/kpc, entropyBG/entropy*10, ls=':', label='entropy / 1000 $keV*cm^2$')
plt.grid(ls=':')
plt.xlabel('radius (kpc)')
plt.xlim(0,150)
plt.legend(loc=5)
plt.title("ICM Profile")
plt.semilogy()
plt.savefig("ICM_profile.pdf")


    

In [ ]:

    

density_derivative = (densityBG[1:]-densityBG[:-1])/(distance[1:]-distance[:-1])
pressure_derivative = (pressureBG[1:]-pressureBG[:-1])/(distance[1:]-distance[:-1])

yr_in_sec = 31557600

plt.plot(distance[1:]/kpc, 1/np.sqrt(pressure_derivative*density_derivative/densityBG[1:]**2)/yr_in_sec/1E6)
#plt.axvline(10)
plt.ylabel('Brunt Vaisala period (Myr)')
plt.xlabel('radius (kpc)')
plt.semilogy()


    

In [ ]:

    

plt.plot(distance/kpc, tempBG)
plt.ylabel('Temperature (K)')
plt.xlabel('radius (kpc)')
ax1 = plt.gca()
ma, mi = ax1.get_ylim()
ax2 = plt.twinx()
keV = 1.1604525E7
ax2.set_ylim(ma/keV, mi/keV)
ax2.set_ylabel('(keV)')


    

In [ ]:

    

plt.plot(distance/kpc, pressureBG)
plt.ylabel('pressure ($dyn/cm^2$)')
plt.xlabel('radius (kpc)')
plt.semilogy()


    

In [ ]:

    

Myr = 3.155760e+13

sim_velJet = 3.0E9
rhoJet = 5E-27

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

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

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

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


    

In [ ]:

    

distance = 10*kpc

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

print(densityBG, pressureBG)