In [ ]:

    

import matplotlib.pyplot as plt
import numpy as np


    

In [ ]:

    

r_temp = 60
rr = np.logspace(np.log10(8), np.log10(400), 17)
Ts4 = 6.4*(1+(rr/r_temp)**4)/(2.1+(rr/r_temp)**4)
Ts3 = 6.4*(1+(rr/r_temp)**3)/(2.1+(rr/r_temp)**3)
Ts2 = 6.4*(1.0+(rr/r_temp)**2)/(2.1+(rr/r_temp)**2)
Ts = 6.4*(1+(rr/100)**3)/(2.1+(rr/100)**3)


    

In [ ]:

    

fig=plt.figure(figsize=(6,5))
plt.plot(rr, Ts, ls='-',  label='r^3 100')
plt.plot(rr, Ts4, ls='-',  label='r^4')
plt.plot(rr, Ts3, ls='--', label='r^3')
plt.plot(rr, Ts2, ls=':',  label='r^2')
plt.semilogx()
plt.xlim(8,400)
plt.legend(loc=2)
plt.axvline(40)
plt.axvline(100)


    

In [ ]:

    

gasConst=8.31447E+07
sim_mu = 0.61
sim_densityBeta = 0.53
sim_rCore = 8.1E23
sim_rCoreT = 1.85E23
sim_Tcore = 3.48E7
sim_Tout = 7.42E7

rr = np.linspace(0,100,100)*3.086E21

grav = gasConst/sim_mu*rr\
           *(-3.*sim_densityBeta/(1.+rr*rr/sim_rCore**2)/sim_rCore**2\
           *sim_Tout*(1.0+(rr/sim_rCoreT)**3)\
           /(sim_Tout/sim_Tcore+(rr/sim_rCoreT)**3)\
           +3.0*rr*sim_Tout*(sim_Tout/sim_Tcore-1.0)*(sim_rCoreT)**3\
           /(sim_Tout/sim_Tcore*(sim_rCoreT)**3+rr**3)**2)
plt.plot(rr/3.086E21, grav)


    

In [ ]: