of course, $\rho $ and $g$ are dependent of radius and the pressure $P$ is a function of density $\rho$
We insert this into the above equation of the hydrostatic equilibrium.
Other useful relations are:
We can write $d\rho / \rho$ as a logarithmic derivative:
$$\frac{d\rho}{\rho} = -\frac{GM(r)}{R_sT}\frac{dr}{r^2} \\ d\ln \rho = - \frac{GM(r)}{R_sT}\frac{dr}{r^2}$$And some more sorting of the terms in the equation:
This is the differential equation for the isothermal sphere. Now to the solution wwhich has to be obtained numerically.
The derivative $d/dr$ changes to $\frac{d}{dx}\sqrt{\frac{4\pi G\rho_c}{R_s T}}$ and eliminates the prefactor $\frac{4\pi G}{R_s T}$ in the equation. We then have the dimensionless differential equation:
$$y'' - \frac{y'^2}{y} + 2\frac{y'}{x} + y^2 = 0$$To get this into a system of two first order differential equations, we do the following:
$$ \frac{dy}{dx}= y' \\ \frac{dy'}{dx}=\frac{y'^2}{y} - 2\frac{y'}{x} - y^2 $$We can use this to integrate the density. SciPy provides a function odeint which takes a vector of derivatives as input and returns the value.
In [22]:
%matplotlib inline
plt.rcParams['figure.figsize']=(9.0,6.0)
plt.rcParams['font.size']=16
plt.rcParams['font.family']='serif'
In [23]:
import numpy as np
import scipy.constants as C
from scipy.integrate import odeint
import matplotlib.pyplot as plt
N=10000 # number of points on R
T0=200. #arbitraty temperature
rho_c = 1e-10
c=np.sqrt((4*C.pi*C.G*rho_c)/(C.R*T0))
#derivative function which is needed for the integration
def deriv(y,r):
return np.array([y[1], y[1]**2/y[0] - 2*y[1]/r - y[0]**2 ])
yinit = np.array([1., 0.]) # initial value vector
radius = np.linspace(0.1,100.,N) # radius r=0. cannot be used, as it is divided by r
y = odeint(deriv, yinit, radius) # calling the integration function
fig1= plt.figure(1)
sub1 = fig1.add_subplot(111)
sub1.plot(radius, y[:,0])
sub1.set_xlabel("radius")
sub1.set_ylabel("density $\\rho / \\rho_c$")
sub1.set_yscale('log')
sub1.set_xscale('log')
sub1.set_ylim(1e-2,1.5)
sub1.set_xlim(1e-1,1e2)
sub1.set_title("Logarithmic plot of the Sphere")
Out[23]:
In [24]:
fig2 = plt.figure(2)
sub2 = fig2.add_subplot(111)
sub2.plot(radius, y[:,0])
sub2.set_xlim(0,20)
sub2.set_title("Linear plot")
sub2.set_xlabel("radius")
sub2.set_ylabel("density $\\rho/\\rho_c$")
Out[24]:
In [24]: