Modelling Stellar Interiors: Jupyter Notebook¶
import math
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
Defining the Dimensionless Lane-Emden equation as: $$\frac{1}{\xi^2} \frac{d}{d\xi} \big(\xi^2 \frac{d\theta}{d\xi} \big) =-\theta^n$$ Where using the Chain rule expands the equation to give: $$\frac{d^2\theta}{d\xi^2}=-\frac{2}{\xi} \frac{d\theta}{d\xi} -\theta^n$$ As this can't be solved analytically, we need to consider what happens at the $\theta_{i+1}$ step, and how it compares when at position $i$, leading the equation: $$\theta_{i+1}=\theta_{i}+ \big[\Delta\xi \big(\frac{d\theta}{d\xi}_{i+1} \big) \big]$$ Where $\Delta \xi$ is a known step value, as well as $\frac{d\theta}{d\xi}$ also being known.
At position $i+1$: $$\big ( \frac{d\theta}{d\xi} \big)_{i+1} =\Delta\xi \big[ \frac{2}{\xi_{i}} \frac{d\theta}{d\xi_{i}}+\theta^n \big] $$
In which all quantities on the Right Hand Side (RHS) are known, and hence can be calculated using numerical integration.
Writing the function to solve the Lane Emden Equation¶
#Need to define a function in order to carry out the iteration over the total extent of the graph.
def function(n):
theta=1 #initial y displacement of the graph
xi=0.0000001 #starts at the y intercept (or as close to the intercept in order to not create a divide by 0 error in the function)
delta_xi=0.0001 #step size- small enough to be accurate, but not too small to break the computer
dtheta_dxi=0 #intial gradient of the graph
theta_values=[] #empty lists for theta and xi
xi_values=[]
while theta >= 0 and xi <= 25: # xi <=25 due to the n=5 going off to infinity
dtheta_dxi=dtheta_dxi-delta_xi*(((2/xi)*dtheta_dxi)+theta**n) #taken from equation
xi=xi+delta_xi #value of xi_{i+1}, is the next value of xi for the function to act on
theta= theta+dtheta_dxi*delta_xi #value of theta_{i+1}, taken from equation
theta_values.append(theta)
xi_values.append(xi)
return (xi_values, theta_values)
xi, theta=function(0) #calling the function for the required values of n (n=0, 1, 2, 3, 4, 5)
xi1, theta1=function(1) #this block was mainly used to check that the function was outputting correct (or seemingly correct) values.
xi2, theta2=function(2)
xi3, theta3=function(3)
xi4, theta4=function(4)
xi5, theta5=function(5)
fig, axis = plt.subplots(figsize=(18,9))
axis.plot(xi, theta, label='n=0') #plotting all the data for the required values of n, with plots of theta against xi
axis.plot(xi1, theta1, label='n=1')
axis.plot(xi2, theta2, label='n=2')
axis.plot(xi3, theta3, label='n=3')
axis.plot(xi4, theta4, label='n=4')
axis.plot(xi5, theta5, label='n=5')
axis.set_xlabel('xi')
axis.set_ylabel('theta')
axis.axis('tight')
axis.legend()
plt.show
Theta intercepts of the Lane Emden graph¶
#Finding the values of xi at theta=0
linear_interpolation=interp1d(theta, xi)
#calculating the theta intercepts through linear interpolation for each value of n (except n=5)
Xint0=interp1d(theta,xi)(0)
Xint1=interp1d(theta1,xi1)(0)
Xint2=interp1d(theta2,xi2)(0)
Xint3=interp1d(theta3,xi3)(0)
Xint4=interp1d(theta4,xi4)(0)
print('The value of xi, at theta=0, when n=0 is: {:.5}'.format(Xint0)) #printing the values to 3 decimal places
print("The value of xi, at theta=0, when n=1 is: {:.5}".format(Xint1))
print("The value of xi, at theta=0, when n=2 is: {:.5}".format(Xint2))
print("The value of xi, at theta=0, when n=3 is: {:.5}".format(Xint3))
print("The value of xi, at theta=0, when n=4 is: {:.6}".format(Xint4))
The theta intercept when n=5 is $\infty$, hence it cannot be calculated through linear interpolation like the other intercepts.
Comparison to the Sun¶
SSM_data=np.loadtxt('../PHY213 Coursework/ssm.txt', skiprows=1)
r_SSM=SSM_data[:,1] #Selecting the r/R_Sun column
rho_SSM=SSM_data[:,3] #Selecting the density column
M_SSM=SSM_data[:,0] #Selecting the Mass column in the data file
P_SSM=SSM_data[:,4] #Selecting the Pressure column in the data file
T_SSM=SSM_data[:,2] #Selecting the Temperature column in the data file
In order to test which polytropic index is the most accurate when compared to the Sun, the dimensionless radius and density must be converted back into their real quantities.
Starting with the equation of mass conservation: $$\frac{dM}{dr} = 4\pi r^2\rho$$ Re-arranging to calculate what the mass of the star is in integral form: $$ M= \int_{0}^{r} 4\pi r^2\rho dr$$ As the value of r can be rewritten as $\alpha=\frac{r}{\xi}$ , so at the surface $R=\alpha\xi_R$, therefore: $$M= 4\pi\alpha^3\rho_C\int_{0}^{\xi} d\xi$$ Therefore: $$M=-4\pi\alpha^3\rho_C \frac{d\theta}{d\xi}$$ Where $\frac{d\theta}{d\xi}$ is known, and $\rho_C$ and $\alpha$ can be calculated.
Using the average density of the sun: $$\rho_{sun}=\frac {3M_{sun}}{4 \pi R^3_{sun}} = 3\rho_C \big|\frac{1}{\xi} \frac{d\theta}{d\xi} \big|_{\xi=\xi_R}$$ $$\rho_C=-\frac {\rho_{sun} \xi_R}{3 \big|\frac{d\theta}{d\xi_R}\big|}$$ When re-arranged for $\rho_C$
With this the central density can be calculated as all the values are known. The value of $\xi_R$ is equal to the x-intercept of the xi, theta graph for each value of n.
As $\rho_C$ is able to be calculated, the value of $rho$ can be calculated using $\rho=\rho_C\theta^n$, and hence a graph of log$\rho$ against radius (in units of solar radii) can be plotted.
Constants and calculating r/R_sun¶
#n=0 has a constant density and hence a constant mass distribution throughout the star, so the graph of density vs radius is a straight line
#Solar units
M_sun=1.99e30 #solar mass
R_sun=6.957e8 #solar radius
rho_sun=(3*M_sun)/(4*np.pi*R_sun**3) #average density of the sun
G=6.67e-11 #Universal Gravitational constant
mH=1.67e-27 #mass of the hydrogen atom
k=1.38e-23 #Boltzmann Constant
mu=9.6
#defining r and alpha (important for all the graphs)
alpha0=R_sun/Xint0
alpha1=R_sun/Xint1
alpha2=R_sun/Xint2
alpha3=R_sun/Xint3
alpha4=R_sun/Xint4
r0=(alpha0*xi)/R_sun
r1=(alpha1*xi1)/R_sun
r2=(alpha2*xi2)/R_sun
r3=(alpha3*xi3)/R_sun
r4=(alpha4*xi4)/R_sun
Mass and Density¶
#defining the central density rho_C
def func(n):
theta=1
xi=0.0000001
delta_xi=0.0001
dtheta_dxi=0
dtheta_dxi_values=[]
while theta >=0 and xi <= 25:
dtheta_dxi=dtheta_dxi-delta_xi*(((2/xi)*dtheta_dxi)+theta**n) #taken from equation
xi=xi+delta_xi #value of xi_{i+1}, is the next value of xi for the function to act on
theta= theta+dtheta_dxi*delta_xi #value of theta_{i+1}, taken from equation
dtheta_dxi_values.append(dtheta_dxi)
return dtheta_dxi_values
rho_C=-(rho_sun*Xint0)/(3*func(0)[-1]) #Due to the x intercepts of each value on n is the value of xi_R
rho_C1=-(rho_sun*Xint1)/(3*func(1)[-1])
rho_C2=-(rho_sun*Xint2)/(3*func(2)[-1])
rho_C3=-(rho_sun*Xint3)/(3*func(3)[-1])
rho_C4=-(rho_sun*Xint4)/(3*func(4)[-1])
#rho_C5=-(rho_sun*Xint5)/(3*func(5)[-1]) no intercept for n=5 as it intercepts at infinity
#Calculations for density
rho=rho_C*np.power(theta, 0) #n=0 hence theta^n=1
rho1=rho_C1*np.power(theta1, 1)
rho2=rho_C2*np.power(theta2, 2)
rho3=rho_C3*np.power(theta3, 3)
rho4=rho_C4*np.power(theta4, 4)
#Graphs require log_10 density to be plotted:
logrho=np.log10(rho)
logrho1=np.log10(rho1)
logrho2=np.log10(rho2)
logrho3=np.log10(rho3)
logrho4=np.log10(rho4)
#Mass calculations
M0=(-4*np.pi*np.power(alpha0, 3)*rho_C*np.power(xi, 2)*func(0))/M_sun
M1=(-4*np.pi*np.power(alpha1, 3)*rho_C1*np.power(xi1, 2)*func(1))/M_sun
M2=(-4*np.pi*np.power(alpha2, 3)*rho_C2*np.power(xi2, 2)*func(2))/M_sun
M3=(-4*np.pi*np.power(alpha3, 3)*rho_C3*np.power(xi3, 2)*func(3))/M_sun
M4=(-4*np.pi*np.power(alpha4, 3)*rho_C4*np.power(xi4, 2)*func(4))/M_sun
Pressure and Temperature¶
Using the Polytropic equation of state: $$P=K\rho^\frac{n+1}{n} = K\rho_C^\frac{n+1}{n} \theta^{n+1}$$ Where K is a known constant from the definition of $\alpha$: $$\alpha^2=\frac{(n+1)K}{4\pi G\rho_C^\frac{n-1}{n}} \therefore K=\frac{4\pi G \rho_C^\frac{n-1}{n}}{\alpha^2 (n+1)}$$ Hence a value of K can be calculated for each value of n, which can then be used to calculate values of the internal Pressure (P) of the star, which was graphed against radius in solar radii
#Calculations for Pressure
#First have to define the equation of K- which depends on every value of n
#Define 4piG as a constant as it doesn't vary- to make the coding easier basically
H=4*np.pi*G
P0=(H*rho_C**2*np.power(theta, 1)*np.power(alpha0, 2))/(0+1)
P1=(H*rho_C1**2*np.power(theta1, 2)*np.power(alpha1, 2))/(1+1)
P2=(H*rho_C2**2*np.power(theta2, 3)*np.power(alpha2, 2))/(2+1)
P3=(H*rho_C3**2*np.power(theta3, 4)*np.power(alpha3, 2))/(3+1)
P4=(H*rho_C4**2*np.power(theta4, 5)*np.power(alpha4, 2))/(4+1)
#Calculations for Temperature
#Have to use the previously defined values of K to calculate temperature, hence can only work for n=1, 2, 3, 4
T0=np.log10((mH*mu*H*rho_C*np.power(alpha0, 2)*np.power(theta, 1))/((0+1)*k))
T1=np.log10((mH*mu*H*rho_C1*np.power(alpha1, 2)*np.power(theta1, 1))/((1+1)*k))
T2=np.log10((mH*mu*H*rho_C2*np.power(alpha2, 2)*np.power(theta2, 1))/((2+1)*k))
T3=np.log10((mH*mu*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T4=np.log10((mH*mu*H*rho_C4*np.power(alpha4, 2)*np.power(theta4, 1))/((4+1)*k))
SSM Data¶
Plotting the graphs¶
fig, axis = plt.subplots(figsize=(18,9))
axis.plot(r0, logrho, label='n=0')
axis.plot(r1, logrho1, label='n=1')
axis.plot(r2, logrho2, label='n=2')
axis.plot(r3, logrho3, label='n=3')
axis.plot(r4, logrho4, label='n=4')
axis.plot(r_SSM, np.log10(rho_SSM*1000), 'k--', label='SSM')
axis.set_xlabel('Radius (r/R_sun)')
axis.set_ylabel('log density')
axis.axis('tight')
axis.legend()
plt.ylim(-8,8) #setting limits on the y axis to improve the aesthetics of the graph
plt.show()
fig,axis=plt.subplots(figsize=(18,9))
axis.plot(r0, M0, label='n=0')
axis.plot(r1, M1, label='n=1')
axis.plot(r2, M2, label='n=2')
axis.plot(r3, M3, label='n=3')
axis.plot(r4, M4, label='n=4')
axis.plot(r_SSM, M_SSM, 'k--', label='SSM')
axis.set_xlabel('Radius (r/R_sun)')
axis.set_ylabel('Mass (M/M_sun)')
axis.axis('tight')
axis.legend()
plt.show()
fig,axis=plt.subplots(figsize=(18,9))
axis.plot(r0, np.log10(P0), label='n=0')
axis.plot(r1, np.log10(P1), label='n=1')
axis.plot(r2, np.log10(P2), label='n=2')
axis.plot(r3, np.log10(P3), label='n=3')
axis.plot(r4, np.log10(P4), label='n=4')
axis.plot(r_SSM, np.log10(P_SSM), 'k--', label='SSM')
axis.set_xlabel('Radius (r/R_sun)')
axis.set_ylabel('Log Pressure (Pa)')
axis.axis('tight')
axis.legend()
plt.ylim(0, 20)
plt.show()
fig,axis=plt.subplots(figsize=(18,9))
axis.plot(r0, T0, label='n=0')
axis.plot(r1, T1, label='n=1')
axis.plot(r2, T2, label='n=2')
axis.plot(r3, T3, label='n=3')
axis.plot(r4, T4, label='n=4')
axis.plot(r_SSM, np.log10(T_SSM), 'k--', label='SSM')
axis.set_xlabel('Radius (r/R_sun)')
axis.set_ylabel('Log Temperature (K)')
axis.axis('tight')
axis.legend()
plt.ylim(5, 10)
plt.show()
