In [ ]:
%matplotlib inline
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 150
plt.rcParams['font.family'] = 'stixgeneral'
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.08567758128e+21
kB = 1.38064852E-16
gasConst = 8.3144598E7
mH = 1.6737352238051868e-24 #g
keV_erg = 1.602176562e-09 #erg
rmin, rmax = 0.1,1000
distance = np.linspace(rmin*kpc, rmax*kpc, rmax-rmin+1)
r = distance
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*gasConst/sim_mu
entropyBG = pressureBG*((densityBG/mH*sim_mu)**(-sim_gammaICM))/e
In [ ]:
from scipy.special import hyp2f1
Msun = 2E33 #g
Mgas = sim_rhoCore/Msun/3*4*np.pi*r[-1]**3*hyp2f1(1.5, 1.5*sim_densityBeta, 2.5, -(r[-1]/sim_rCore)**2)
print('%e' % Mgas)
In [ ]:
keV = 1.1604525E7 # K
# normalization
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(rmin,rmax)
plt.legend(ncol=2)
#plt.title("ICM Profile")
plt.semilogy()
plt.semilogx()
In [ ]:
grav = gasConst/sim_mu*r \
*(-3.*sim_densityBeta/(1.+r*r/sim_rCore**2)/sim_rCore**2 \
*sim_Tout*(1.0+(r/sim_rCoreT)**3)\
/(sim_Tout/sim_Tcore+(r/sim_rCoreT)**3)\
+3.0*r*sim_Tout*(sim_Tout/sim_Tcore-1.0)*(sim_rCoreT)**3\
/(sim_Tout/sim_Tcore*(sim_rCoreT)**3+r**3)**2)
plt.plot(r/kpc, grav, label='g')
plt.legend()
#plt.semilogx()
In [ ]:
i0 = 0
print('initial radius = %.1f kpc' % (r[i0]/kpc))
M = 1
rho_ad = densityBG[i0]*(pressureBG/pressureBG[i0])**(1/sim_gammaICM)
plt.plot(r/kpc, rho_ad, label=r'$\rho_{ad}$')
plt.plot(r/kpc, densityBG, label=r'$\rho$')
plt.plot(r/kpc, rho_ad-densityBG, ':', label=r'$\Delta \rho$')
plt.axvline(r[i0]/kpc, ls=':', color='k')
plt.legend()
plt.grid()
In [ ]:
V_ad = M/densityBG[i0]*(pressureBG/pressureBG[i0])**(-1/sim_gammaICM)
Fg = (rho_ad-densityBG)*V_ad*grav
plt.plot(r/kpc, -Fg, label=r'$|F_g|$')
plt.axvline(r[i0]/kpc, ls=':', color='k')
plt.legend()
In [ ]:
integrate.simps(-Fg[i0:], r[i0:])
In [ ]:
T_ad = tempBG[i0]*(pressureBG/pressureBG[i0])/(rho_ad/densityBG[i0])
plt.plot(r/kpc, T_ad, label=r'$T_{ad}$')
plt.plot(r/kpc, tempBG, label=r'$T$')
plt.plot(r/kpc, tempBG-T_ad, label=r'$\Delta T$')
plt.axvline(r[i0]/kpc, ls=':', color='k')
plt.ylim(0, 1.3E8)
plt.legend()
In [ ]:
W = integrate.cumtrapz(-Fg[i0:], r[i0:])
W
In [ ]:
Eth = 1/(sim_gammaICM-1)*M/sim_mu*gasConst*(tempBG-T_ad)
plt.plot(r/kpc, Eth, label='thermal energy')
plt.plot(r[i0+1:]/kpc, W, label='work')
plt.axvline(r[i0]/kpc, ls=':', color='k')
plt.legend()
plt.grid()
In [ ]:
xi = (Eth[i0+1:]-W)/W
print(xi[-1])
plt.plot(r[i0+1:]/kpc, xi)
plt.semilogy()
plt.axhline(1, ls=':', color='k')
plt.grid()
In [ ]:
ls = {10: ['solid', 2],
20: ['dotted', 2],
40: ['dashed', 1],
65: ['solid', 1],
100: ['dotted', 1]}
colors = {10: 'C0', 20: 'C1', 40: 'C2', 65: 'C3', 100: 'C4'}
radii = colors.keys()
plt.figure(figsize=(4,2.5))
for r0 in radii:
i0 = r0-rmin
print('initial radius = %.1f kpc' % (r[i0]/kpc))
V_ad = M/densityBG[i0]*(pressureBG/pressureBG[i0])**(-1/sim_gammaICM)
rho_ad = densityBG[i0]*(pressureBG/pressureBG[i0])**(1/sim_gammaICM)
Fg = (rho_ad-densityBG)*V_ad*grav
W = integrate.cumtrapz(-Fg[i0:], r[i0:])
T_ad = tempBG[i0]*(pressureBG/pressureBG[i0])/(rho_ad/densityBG[i0])
Eth = 1/(sim_gammaICM-1)*M/sim_mu*gasConst*(tempBG[i0+1:]-T_ad[i0+1:])
xi = (Eth-W)/W
plt.plot(r[i0+1:]/kpc, xi, label=r'$R_0$ = %.1f kpc' % (r[i0]/kpc),
color=colors[r0], ls=ls[r0][0], lw=ls[r0][1])
plt.legend(fontsize=8)
plt.semilogy()
#plt.semilogx()
plt.axhline(1, ls=(0, (1, 1)), lw=1, color='k')
plt.grid(alpha=0.4, ls='--')
plt.ylabel(r'$\Delta\xi_{\rm{max}} \equiv \frac{E_{th}-W}{W}$')
plt.xlabel(r'R [kpc]')
plt.ylim(0.1, 10)
plt.xlim(0, 400)
plt.savefig('gain_efficiency_perseus.pdf', bbox_inches='tight')
In [ ]:
ls = {10: ['solid', 2],
20: ['dotted', 2],
40: ['dashed', 1],
65: ['solid', 1],
100: ['dotted', 1]}
colors = {10: 'C0', 20: 'C1', 40: 'C2', 65: 'C3', 100: 'C4'}
radii = colors.keys()
plt.figure(figsize=(4,2.5))
for r0 in radii:
i0 = r0-rmin
print('initial radius = %.1f kpc' % (r[i0]/kpc))
V_ad = M/densityBG[i0]*(pressureBG/pressureBG[i0])**(-1/sim_gammaICM)
rho_ad = densityBG[i0]*(pressureBG/pressureBG[i0])**(1/sim_gammaICM)
Fg = (rho_ad-densityBG)*V_ad*grav
W = integrate.cumtrapz(-Fg[i0:], r[i0:])
T_ad = tempBG[i0]*(pressureBG/pressureBG[i0])/(rho_ad/densityBG[i0])
Eth = 1/(sim_gammaICM-1)*M/sim_mu*gasConst*(tempBG[i0+1:]-T_ad[i0+1:])
#dE = (Eth-W)/M/(-grav[i0])/r[i0]
dE = (Eth-W)/M*sim_mu*mH/keV_erg
plt.plot(r[i0+1:]/kpc, dE, label=r'$R_0$ = %.1f kpc' % (r[i0]/kpc),
color=colors[r0], ls=ls[r0][0], lw=ls[r0][1])
plt.legend(fontsize=8, ncol=2, frameon=False)
#plt.semilogy()
#plt.semilogx()
plt.axhline(0, ls='-', lw=1, color='k', alpha=0.7)
plt.grid(alpha=0.4, ls='--')
#plt.ylabel(r'$\Delta \epsilon_{\rm{max}} \equiv \frac{E_{th}-W}{M g_0 R_0}$')
plt.ylabel(r'$\Delta \epsilon_{\rm{max}} \equiv (E_{th}-W)\ \frac{\mu m_p}{M}$ [keV/particle]')
plt.xlabel(r'radius R [kpc]')
plt.ylim(-4, 4)
plt.xlim(0, 400)
plt.savefig('energy_gain_perseus.pdf', bbox_inches='tight')
In [ ]:
r0_ind = np.arange(2, 300)
plt.figure(figsize=(4,3))
dR_max = np.zeros(len(r0_ind))
dE_max = np.zeros(len(r0_ind))
for i, r0 in enumerate(r0_ind):
i0 = r0-rmin
V_ad = M/densityBG[i0]*(pressureBG/pressureBG[i0])**(-1/sim_gammaICM)
rho_ad = densityBG[i0]*(pressureBG/pressureBG[i0])**(1/sim_gammaICM)
Fg = (rho_ad-densityBG)*V_ad*grav
W = integrate.cumtrapz(-Fg[i0:], r[i0:])
T_ad = tempBG[i0]*(pressureBG/pressureBG[i0])/(rho_ad/densityBG[i0])
Eth = 1/(sim_gammaICM-1)*M/sim_mu*gasConst*(tempBG[i0+1:]-T_ad[i0+1:])
dE = (Eth-W)/M*sim_mu*mH/keV_erg
imax = np.argmax(dE)
dR_max[i] = r[imax] - r0
dE_max[i] = dE[imax]
c1 = 'navy'
plt.plot(r[r0_ind]/kpc, dE_max, c=c1, label=r'$\Delta \epsilon_{\rm{max,peak}}$')
plt.legend(frameon=False, loc=3)
plt.ylabel(r'$\Delta \epsilon_{\rm{max,peak}}$ [keV/particle]', color=c1)
plt.tick_params(axis='y', labelcolor=c1)
plt.ylim(0, 4)
plt.xlabel(r'initial radius $R_0$ [kpc]')
plt.twinx()
c2 = 'orangered'
plt.plot(r[r0_ind]/kpc, dR_max/kpc, c=c2, ls='--', label=r'$\Delta R_{\rm{peak}}$')
plt.legend(frameon=False, loc=4)
plt.ylim(0, 250)
plt.tick_params(axis='y', labelcolor=c2)
plt.ylabel(r'$\Delta R_{\rm{peak}} \equiv R_{\rm{peak}} - R_0$ [kpc]', color=c2)
plt.xlim(0, 200)
plt.savefig('energy_gain_peak_perseus.pdf', bbox_inches='tight')
In [ ]:
ls = {10: ['solid', 2],
20: ['dotted', 2],
40: ['dashed', 1],
65: ['solid', 1],
100: ['dotted', 1]}
colors = {10: 'C0', 20: 'C1', 40: 'C2', 65: 'C3', 100: 'C4'}
radii = colors.keys()
plt.figure(figsize=(4,2.5))
for r0 in radii:
i0 = r0-rmin
print('initial radius = %.1f kpc' % (r[i0]/kpc))
V_ad = M/densityBG[i0]*(pressureBG/pressureBG[i0])**(-1/sim_gammaICM)
rho_ad = densityBG[i0]*(pressureBG/pressureBG[i0])**(1/sim_gammaICM)
Fg = (rho_ad-densityBG)*V_ad*grav
W = integrate.cumtrapz(-Fg[i0:], r[i0:])
T_ad = tempBG[i0]*(pressureBG/pressureBG[i0])/(rho_ad/densityBG[i0])
Eth = 1/(sim_gammaICM-1)*M/sim_mu*gasConst*(tempBG[i0+1:]-T_ad[i0+1:])
xi = (Eth)/W
plt.plot(r[i0+1:]/kpc, xi, label=r'$R_0$ = %.1f kpc' % (r[i0]/kpc),
color=colors[r0], ls=ls[r0][0], lw=ls[r0][1])
plt.legend(fontsize=8, frameon=False)
plt.semilogy()
#plt.semilogx()
plt.axhline(1, ls=(0, (1, 1)), lw=1, color='k')
plt.grid(alpha=0.4, ls='--')
plt.ylabel(r'$\xi_{\rm{max}} \equiv \frac{E_{th}}{W}$')
plt.xlabel(r'radius R [kpc]')
plt.ylim(0.5, 20)
plt.xlim(0, 400)
plt.savefig('energy_efficiency_perseus.pdf', bbox_inches='tight')
In [ ]:
(gasConst*tempBG/sim_mu) / (grav*r/(3*sim_densityBeta-d))
In [ ]: