Jupyter notebook PHY213 Coursework/stellarinteriors.ipynb

Frederick Ashman - PHY213: Stellar Structure & Evolution

PHY213 Coursework/stellarinteriors.ipynb

⁷⁰ views

Kernel: Python 2 (SageMath)

Modelling Stellar Interiors: Jupyter Notebook

In [1]:

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

In [3]:

#Need to define a function in order to carry out the iteration over the total extent of the graph.
def polytropes(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=polytropes(0) #calling the function for the required values of n (n=0, 1, 2, 3, 4, 5)
xi1, theta1=polytropes(1) #this block was mainly used to check that the function was outputting correct (or seemingly correct) values.
xi2, theta2=polytropes(2)
xi3, theta3=polytropes(3)
xi4, theta4=polytropes(4)
xi5, theta5=polytropes(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$', fontsize=20)
axis.set_ylabel('$\\theta$', fontsize=20)
axis.axis('tight')
axis.legend()
plt.show

Out[3]:

<function matplotlib.pyplot.show>

Theta intercepts of the Lane Emden graph

In [4]:

#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))

Out[4]:

The value of xi, at theta=0, when n=0 is:   2.449
The value of xi, at theta=0, when n=1 is:   3.141
The value of xi, at theta=0, when n=2 is:   4.352
The value of xi, at theta=0, when n=3 is:   6.897
The value of xi, at theta=0, when n=4 is:   14.973

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

In [5]:

#In order to compare the models to the known observational data from the Sun, the Standard Solar Model data must be imported:
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])/10 #Selecting the Pressure column in the data file
T_SSM=SSM_data[:,2] #Selecting the Temperature column in the data file
mu_SSM=SSM_data[:,-1]


#calculating central values to compare to the polytropes
rho_SSM_C=rho_SSM[0]
P_SSM_C=P_SSM[0]
T_SSM_C=T_SSM[0]
print(rho_SSM_C)
print(P_SSM_C)
print(T_SSM_C)

Out[5]:

151.0
2.34e+16
15500000.0

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

In [8]:

#Useful constants to define
M_sun=1.99e30 #solar mass (kg)
R_sun=6.957e8 #solar radius (m)
rho_sun=(3*M_sun)/(4*np.pi*R_sun**3) #average density of the sun (kgm^-3)
G=6.67e-11 #Universal Gravitational constant (m^3 kg^-1 s^-2)
mH=1.67e-27 #mass of the hydrogen atom (kg)
k=1.38e-23 #Boltzmann Constant (m^2 kgs^-2 K^-1)
mu=0.62 #is a good approximation for the majority of the Star, however does deviate away from the true value of mu towards the central regions of the star (unitless)

#defining r (essential for all the graphs)
#owing that alpha=R_sun/Xint and r=alpha*xi/R_sun, simplifying r to r=R_sun*xi/Xint*R_sun, therefore r=xi/Xint
r0=xi/Xint0
r1=xi1/Xint1
r2=xi2/Xint2
r3=xi3/Xint3
r4=xi4/Xint4

alpha0=R_sun/Xint0
alpha1=R_sun/Xint1
alpha2=R_sun/Xint2
alpha3=R_sun/Xint3
alpha4=R_sun/Xint4

Mass and Density

In [9]:

#defining the central density rho_C
def dtheta_dxi(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*dtheta_dxi(0)[-1]) #Due to the x intercepts of each value on n is the value of xi_R
rho_C1=-(rho_sun*Xint1)/(3*dtheta_dxi(1)[-1])
rho_C2=-(rho_sun*Xint2)/(3*dtheta_dxi(2)[-1])
rho_C3=-(rho_sun*Xint3)/(3*dtheta_dxi(3)[-1])
rho_C4=-(rho_sun*Xint4)/(3*dtheta_dxi(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)*dtheta_dxi(0))/M_sun
M1=(-4*np.pi*np.power(alpha1, 3)*rho_C1*np.power(xi1, 2)*dtheta_dxi(1))/M_sun
M2=(-4*np.pi*np.power(alpha2, 3)*rho_C2*np.power(xi2, 2)*dtheta_dxi(2))/M_sun
M3=(-4*np.pi*np.power(alpha3, 3)*rho_C3*np.power(xi3, 2)*dtheta_dxi(3))/M_sun
M4=(-4*np.pi*np.power(alpha4, 3)*rho_C4*np.power(xi4, 2)*dtheta_dxi(4))/M_sun

#calculating the central value of density to compare to the SSM data
print(rho_C3)

Out[9]:

/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:34: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:36: RuntimeWarning: invalid value encountered in log10

76463.0584822

Plotting for Density and Mass

In [11]:

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_{\odot}$)', fontsize=20)
axis.set_ylabel('$log_{10} \\rho $', fontsize=20)
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_{\odot}$)' ,fontsize=20)
axis.set_ylabel('Mass ($M/M_{\odot}$)' ,fontsize=20)
axis.axis('tight')
axis.legend()
plt.show()

Out[11]:

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

In [13]:

#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
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))

#calculating values for the central values to compare to the SSM
P_C3= P3[0]
T_C3= T3[0]
print(P_C3)
print(T_C3)

Out[13]:

1.24645375663e+16
7.08745359253

/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:15: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:16: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in log10

Graphs for Pressure and Temperature

In [14]:

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_{\odot}$)', fontsize=20)
axis.set_ylabel('$log_{10} P$', fontsize=20)
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_{\odot}$)', fontsize=20)
axis.set_ylabel('$log_{10} T$', fontsize=20)
axis.axis('tight')
axis.legend()
plt.ylim(5, 8)
plt.show()

Out[14]:

/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in log10
  from ipykernel import kernelapp as app
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in log10

In [15]:

#Calculating and graphing the temperature graphs using the value of mu taken from the SSM data, for just the n=3 index
#mu_SSM values are taken at intervals of 0.1R_0

mu_0=0.961
mu_1=0.729
mu_2=0.675
mu_3=0.651
mu_4=0.640
mu_5=0.635
mu_6=0.632
mu_7=0.630
mu_8=0.629
mu_9=0.627
mu_10=0.613

#calculating for just the n=3 polytrope - as for the other models this was the closest to the SSM data.

T3_0=np.log10((mH*mu_0*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_1=np.log10((mH*mu_1*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_2=np.log10((mH*mu_2*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_3=np.log10((mH*mu_3*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_4=np.log10((mH*mu_4*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_5=np.log10((mH*mu_5*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_6=np.log10((mH*mu_6*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_7=np.log10((mH*mu_7*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_8=np.log10((mH*mu_8*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_9=np.log10((mH*mu_9*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))
T3_10=np.log10((mH*mu_10*H*rho_C3*np.power(alpha3, 2)*np.power(theta3, 1))/((3+1)*k))

fig,axis=plt.subplots(figsize=(18,9))
axis.plot(r3, T3_0, label='mu=0.961')
axis.plot(r3, T3_1, label='mu=0.729')
axis.plot(r3, T3_2, label='mu=0.675')
axis.plot(r3, T3_3, label='mu=0.651')
axis.plot(r3, T3_4, label='mu=0.640')
axis.plot(r3, T3_5, label='mu=0.635')
axis.plot(r3, T3_6, label='mu=0.632')
axis.plot(r3, T3_7, label='mu=0.630')
axis.plot(r3, T3_8, label='mu=0.629')
axis.plot(r3, T3_9, label='mu=0.627')
axis.plot(r3, T3_10, label='mu=0.613')
axis.plot(r_SSM, np.log10(T_SSM), 'k--', label='SSM')
axis.set_xlabel('Radius $(r/R_{\odot})$', fontsize=20)
axis.set_ylabel('$log_{10}(T)$ $(K)$', fontsize=20)
axis.axis('tight')
axis.legend()
plt.ylim(5, 8)
plt.show()

Out[15]:

/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:20: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:21: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:22: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:23: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:24: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:25: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:26: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:27: RuntimeWarning: invalid value encountered in log10
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:28: RuntimeWarning: invalid value encountered in log10

In [21]:

#calculating the min, max, and most comparable values of central temperature to quantitively compare to the SSM data
#min value occurs at mu=0.613, max value at mu=0.961
T_C_min=T3_10[0]
T_C_max=T3_0[0]
T_C_comparable=T3_1[0]

print(T_C_min, T_C_max, T_C_comparable)

Out[21]:

(7.0825223775478729, 7.2777852906980032, 7.1577894313474326)

In [0]: