In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 100
import numpy as np
from scipy.special import gamma

In [ ]:
kpc = 3.0857E21
nc = 1E-2
beta = 0.7
rc = 10
rr = np.logspace(-1, 4, 100)

metallicity = 0.5
# X-ray 0.1~2.4keV metal cooling rate = 2.95E-23 erg cm^3 s^{-1}
#             primordial cooling rate = 7.86E-24 erg cm^3 s^{-1}
Lambda = 2.95E-23*metallicity + 7.86E-24


S = 2*nc**2*Lambda*np.sqrt(np.pi)/2*gamma(3*beta-0.5)/gamma(3*beta)*rc*kpc \
    *(1+(rr/rc)**2)**(-3*beta+0.5)
plt.plot(rr, S)
plt.xlabel('r (kpc)')
plt.ylabel(r'S (erg cm$^{-2}$ s$^{-1}$)')
plt.semilogx()
plt.semilogy()

In [ ]:
def Lx(rr, rc, beta, nc, metallicity=0.5):
    Lambda = 2.95E-23*metallicity + 7.86E-24
    # X-ray 0.1~2.4keV metal cooling rate = 2.95E-23 erg cm^3 s^{-1}
    #             primordial cooling rate = 7.86E-24 erg cm^3 s^{-1}
    L = 2*np.pi*nc**2*(rc*kpc)**3*Lambda*np.sqrt(np.pi)/2*gamma(3*beta-0.5)/gamma(3*beta)
    if beta == 0.5:
        return L*np.log(1+(rr/rc)**2)
    else:
        return L*1/(1.5-3*beta)*((1+(rr/rc)**2)**(-3*beta+1.5)-1)
        

for beta in [0.3, 0.4, 0.5, 0.6, 0.7]:
    nc = 0.04
    rc = 26
    L = Lx(rr, rc, beta, nc)
    plt.plot(rr, L, label=r'$\beta=%.1f$' % beta)

plt.legend()
plt.semilogx()
plt.semilogy()
plt.xlabel('r (kpc)')
plt.ylabel(r'L (erg s$^{-1}$)')
plt.grid()

In [ ]:
Lx(100, 100, 0.7, 1E-2)

In [ ]: