SharedPHY213 Coursework / stellarinteriors.ipynbOpen in CoCalc
Jupyter notebook PHY213 Coursework/stellarinteriors.ipynb

Modelling Stellar Interiors: Jupyter Notebook

#import appropriate packages
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Part 1

Lane emden equation: d2θdξ2=[2ξdθdξ]θn\frac{d^2\theta}{d\xi^2} = -[\frac{2}{\xi}\frac{d\theta}{d\xi}] - \theta^n

#set up function which can solve the lame emden equation for any value of 'n'.
def solve(n):
    """Solves the Lame-Emdem equation for any value of n.
    
    args: n
    """
    #create empty lists for theta and xi values to store them into
    theta_values = []
    xi_values = []
    
    #define values of xi, delta xi, the gradient and theta
    xi = 0.00001
    d_xi = 0.001
    dtheta = 0
    theta = 1
    
    #create a variable to store xi
    xi_now = xi
    
    #use while loop to find theta and xi values until theta=0
    while (theta >= 0) and (xi_now < 20):
        
        #increase xi value by small amount
        xi_now = xi_now + d_xi
    
        #calculate values after small increase in xi
        dtheta_next = dtheta - (((2/xi_now)*dtheta)+theta**n)*d_xi
        theta_next = theta + dtheta_next*d_xi

        #update the old values to be the new ones
        dtheta = dtheta_next
        theta  = theta_next
        
        #store these values in list
        theta_values.append(theta)
        xi_values.append(xi_now)
    
    #convert lists to arrays to make it easier to deal with
    xi_values = np.array(xi_values)
    theta_values = np.array(theta_values)
    return (xi_values, theta_values)

#call the function to find the theata and xi values for each n
xi_0, theta_0 = solve(0)
xi_1, theta_1 = solve(1)
xi_2, theta_2 = solve(2)
xi_3, theta_3 = solve(3)
xi_4, theta_4 = solve(4)
xi_5, theta_5 = solve(5)
#plot the values of xi vs theta for each n value
fig, axis = plt.subplots(figsize = (9,5))
axis.plot(xi_0, theta_0, label = 'n=0')
axis.plot(xi_1, theta_1, label = 'n=1')
axis.plot(xi_2, theta_2, label = 'n=2')
axis.plot(xi_3, theta_3, label = 'n=3')
axis.plot(xi_4, theta_4, label = 'n=4')
axis.plot(xi_5, theta_5, label = 'n=5')
#set limits on the axes
axis.set_ylim(0)
axis.set_xlim(0, 20)
#give title and axes labels
axis.set_title('Numerical Solutions to the Lame Emden Equation')
axis.set_ylabel('Dimensionless density')
axis.set_xlabel('Dimensionless radius')
#add a legend
axis.legend()
#show
plt.show()
#find the intercept of each plot, by using slicing to find the last known value of xi which is the intercept (theta =0)
int_0 = xi_0[-1]
int_1 = xi_1[-1]
int_2 = xi_2[-1]
int_3 = xi_3[-1]
int_4 = xi_4[-1]
int_5 = xi_5[-1]

print 'Values of xi at theta = 0 for various n are as follows'
print 'n=0:', '%.3f' %int_0
print 'n=1:', '%.3f' %int_1
print 'n=2:', '%.3f' %int_2
print 'n=3:', '%.3f' %int_3
print 'n=4:', '%.3f' %int_4
print 'n=5: unknown, as it never crosses the x-axis'
Values of xi at theta = 0 for various n are as follows n=0: 2.448 n=1: 3.141 n=2: 4.353 n=3: 6.900 n=4: 14.993 n=5: unknown, as it never crosses the x-axis

Part 2

Standard Solar Model

#Standard solar model can be found using the data from the file ssm.txt
#import data
ssm_data = np.loadtxt("../PHY213 Coursework/ssm.txt", skiprows=1)
#find radius from the SSM
r_ssm = ssm_data[:,1]
#find mass from the SSM
m_ssm = ssm_data[:,0]
#Find temperature from the SSM
T_ssm = ssm_data[:,2]
log_T_ssm = np.log10(T_ssm)
#Find density from the SSM and convert into the aprropriate units (kgcm^-3)
rho_ssm = ssm_data[:,3]
rho_ssm = rho_ssm*1000
log_rho_ssm = np.log10(rho_ssm)
#Find pressure from the SSM and convert into the appropriate units  (kgcm^-2)
P_ssm = ssm_data[:,4]
P_ssm = P_ssm/10
log_P_ssm = np.log10(P_ssm)
#Find the value of mu from the SSM
mu_ssm = ssm_data[:,12]

Alpha:

α=RξR \alpha = \frac{R}{\xi_R}

#radius of the sun (known)
R = 6.957e8

#find the value of alpha for each n
alpha_0 = R/int_0
alpha_1 = R/int_1
alpha_2 = R/int_2
alpha_3 = R/int_3
alpha_4 = R/int_4
alpha_5 = R/int_5

Radius:

r=αξ r = \alpha\xi

#find array of radius values for each n

r_0 = (xi_0*alpha_0)/R
r_1 = (xi_1*alpha_1)/R
r_2 = (xi_2*alpha_2)/R
r_3 = (xi_3*alpha_3)/R
r_4 = (xi_4*alpha_4)/R
r_5 = (xi_5*alpha_5)/R

#check
print(r_0)
[ 4.12580014e-04 8.21075077e-04 1.22957014e-03 ..., 9.99183010e-01 9.99591505e-01 1.00000000e+00]

Density

Central Density: 3Mo4πro3=3ρc1ξdθdξ \frac{3M_o}{4\pi r_o^3} = -3\rho_c|\frac{1}{\xi}\frac{d\theta}{d\xi}|

ρc=Moξ4πro3dξdθ \rho_c = -\frac{M_o\xi}{4\pi r_o^3}\frac{d\xi}{d\theta}

Density: ρ=ρcθn\rho = \rho_c\theta^n

Log Density: log(ρ)=log(ρcθn)log(\rho) = log(\rho_c\theta^n)

#set up function to find the gradient for any value of n.
def gradient(n):
    """Finds the gradient (dtheta/dxi) for any value of n.
    
    args: n
    """
    #create empty lists for theta, xi and dtheta values to store them into
    theta_values = []
    xi_values = []
    dtheta_values = []
    
    #define values of xi, delta xi, the gradient and theta
    xi = 0.00001
    d_xi = 0.001
    dtheta = 0
    theta = 1
    
    #create a variable to store xi
    xi_now = xi
    
    #use while loop to find theta and xi values until theta=0
    while (theta >= 0) and (xi_now < 20):
        
        #increase xi value by small amount
        xi_now = xi_now + d_xi
    
        #calculate values after small increase in xi
        dtheta_next = dtheta - (((2/xi_now)*dtheta)+theta**n)*d_xi
        theta_next = theta + dtheta_next*d_xi

        #update the old values to be the new ones
        dtheta = dtheta_next
        theta  = theta_next
        
        #store these values in list
        theta_values.append(theta)
        xi_values.append(xi_now)
        dtheta_values.append(dtheta)
    
    #convert lists to arrays to make it easier to deal with
    xi_values = np.array(xi_values)
    theta_values = np.array(theta_values)
    dtheta_values = np.array(dtheta_values)
    return (dtheta_values)

#calculate gradient by calling this function for each value of n.
dtheta_0 = gradient(0)
dtheta_1 = gradient(1)
dtheta_2 = gradient(2)
dtheta_3 = gradient(3)
dtheta_4 = gradient(4)
dtheta_5 = gradient(5)
#Mass of the sun (known)
M_sun = 1.989e30

#Find the gradient (dtheta/dxi) at the surface of the sun for each n
dtheta_R_0 = dtheta_0[-1]
dtheta_R_1 = dtheta_1[-1]
dtheta_R_2 = dtheta_2[-1]
dtheta_R_3 = dtheta_3[-1]
dtheta_R_4 = dtheta_4[-1]
dtheta_R_5 = dtheta_5[-1]

#find the central density of the sun for each n using the equation stated above
rho_c_0 = -M_sun*int_0/(4*np.pi*R**3*dtheta_R_0)
rho_c_1 = -M_sun*int_1/(4*np.pi*R**3*dtheta_R_1)
rho_c_2 = -M_sun*int_2/(4*np.pi*R**3*dtheta_R_2)
rho_c_3 = -M_sun*int_3/(4*np.pi*R**3*dtheta_R_3)
rho_c_4 = -M_sun*int_4/(4*np.pi*R**3*dtheta_R_4)
rho_c_5 = -M_sun*int_5/(4*np.pi*R**3*dtheta_R_5)

#find array of density values for each value of n using the equation stated above
rho_0 = rho_c_0*theta_0**0
rho_1 = rho_c_1*theta_1**1
rho_2 = rho_c_2*theta_2**2
rho_3 = rho_c_3*theta_3**3
rho_4 = rho_c_4*theta_4**4
rho_5 = rho_c_5*theta_5**5

#log these density values to make the plot of the results more readable
log_rho_0 = np.log10(rho_0)
log_rho_1 = np.log10(rho_1)
log_rho_2 = np.log10(rho_2)
log_rho_3 = np.log10(rho_3)
log_rho_4 = np.log10(rho_4)
log_rho_5 = np.log10(rho_5)

/ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:30: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:32: RuntimeWarning: invalid value encountered in log10
#plot the log of density against radius for each n
fig, axis = plt.subplots(figsize = (9,5))
axis.plot(r_0, log_rho_0, label = 'n=0')
axis.plot(r_1, log_rho_1, label = 'n=1')
axis.plot(r_2, log_rho_2, label = 'n=2')
axis.plot(r_3, log_rho_3, label = 'n=3')
axis.plot(r_4, log_rho_4, label = 'n=4')
axis.plot(r_5, log_rho_5, label = 'n=5')
axis.set_ylim(0, 7)
axis.set_ylabel('log(rho) (kgm^-3)')
axis.set_xlabel('r/R_o')
axis.legend()
plt.show()

Mass

Mass: M=4πα3ρcξ2dθdξ M = -4\pi\alpha^3\rho_c\xi^2\frac{d\theta}{d\xi}

#Apply the formula for mass to the varaibles found
M_0 = -4*np.pi*alpha_0**3*rho_c_0*xi_0**2*dtheta_0/M_sun
M_1 = -4*np.pi*alpha_1**3*rho_c_1*xi_1**2*dtheta_1/M_sun
M_2 = -4*np.pi*alpha_2**3*rho_c_2*xi_2**2*dtheta_2/M_sun
M_3 = -4*np.pi*alpha_3**3*rho_c_3*xi_3**2*dtheta_3/M_sun
M_4 = -4*np.pi*alpha_4**3*rho_c_4*xi_4**2*dtheta_4/M_sun
M_5 = -4*np.pi*alpha_5**3*rho_c_5*xi_5**2*dtheta_5/M_sun
#plot the mass values for each n against radius
fig, axis = plt.subplots(figsize = (9,5))
axis.plot(r_0, M_0, label = 'n=0')
axis.plot(r_1, M_1, label = 'n=1')
axis.plot(r_2, M_2, label = 'n=2')
axis.plot(r_3, M_3, label = 'n=3')
axis.plot(r_4, M_4, label = 'n=4')
axis.plot(r_5, M_5, label = 'n=5')
axis.set_ylim(0,1.0)
axis.set_ylabel('M/M_o')
axis.set_xlabel('r/R_o')
axis.legend(loc = 'lower right')
plt.show()

Pressure

Pressure: P=Kρn+1n=Kρcn+1nθn+1 P = K\rho^{\frac{n+1}{n}} = K\rho_c^{\frac{n+1}{n}}\theta^{n+1}

K: α=(n+1)K4πGρc(n1)n \alpha = \sqrt{\frac{(n+1)K}{4\pi G\rho_c^{\frac{(n-1)}{n}}}} K=α24πGρc(n1)n(n+1) K = \frac{\alpha^{2}4\pi G\rho_c^{\frac{(n-1)}{n}}}{(n+1)}

Combining these 2 equations: P=α24πGρc2θn+1(n+1) P = \frac{\alpha^{2}4\pi G\rho_c^{2}\theta^{n+1}}{(n+1)}

#Value of the gravitational constant
G = 6.672e-11

#Input all variables into the combined equation for pressure for each value of n
P_0 = (alpha_0**2*4*np.pi*G*rho_c_0**2*theta_0**(0+1))/(0+1)
P_1 = (alpha_1**2*4*np.pi*G*rho_c_1**2*theta_1**(1+1))/(1+1)
P_2 = (alpha_2**2*4*np.pi*G*rho_c_2**2*theta_2**(2+1))/(2+1)
P_3 = (alpha_3**2*4*np.pi*G*rho_c_3**2*theta_3**(3+1))/(3+1)
P_4 = (alpha_4**2*4*np.pi*G*rho_c_4**2*theta_4**(4+1))/(4+1)
P_5 = (alpha_5**2*4*np.pi*G*rho_c_5**2*theta_5**(5+1))/(5+1)

#log these pressure values to make a better plot
logP_0 = np.log10(P_0)
logP_1 = np.log10(P_1)
logP_2 = np.log10(P_2)
logP_3 = np.log10(P_3)
logP_4 = np.log10(P_4)
logP_5 = np.log10(P_5)
/ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:15: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in log10
#plot the log of the pressure for each n value against radius
fig, axis = plt.subplots(figsize = (9,5))
axis.plot(r_0, logP_0, label = 'n=0')
axis.plot(r_1, logP_1, label = 'n=1')
axis.plot(r_2, logP_2, label = 'n=2')
axis.plot(r_3, logP_3, label = 'n=3')
axis.plot(r_4, logP_4, label = 'n=4')
axis.plot(r_5, logP_5, label = 'n=5')
axis.set_ylim(0)
axis.set_ylabel('logP (Nm^-2)')
axis.set_xlabel('r/R_o')
axis.legend(loc = 'lower left')
plt.show()

Temperature

Equate the equation of state of an ideal gas with the polytropic equation of state:

P=kTρmHμK=Kρc(n+1)nθn+1 P = \frac{kT\rho}{m_H\mu K} = K\rho_c^{\frac{(n+1)}{n}}\theta^{n+1}

Rearrange and substitute with density equation.

T=mHμKρc1nθk T = \frac{m_H\mu K \rho_c^{\frac{1}{n}}\theta}{k}

Input the K equation:

T=mHμα24πGρcθk(n+1) T = \frac{m_H\mu\alpha^{2}4\pi G\rho_c\theta}{k(n+1)}

#mass of hydrogen atom
m_H = 1.674e-27

#value of k (Boltzmann constant)
k = 1.381e-23
#to find the value of mu, plot it against each radius value to give a visual representation
fig, axis = plt.subplots(figsize = (6,3))
axis.plot(r_ssm, mu_ssm, 'k')
axis.set_ylabel('mu')
axis.set_xlabel('r/R_o')
plt.show()
#from this we can see the median mu value would be acceptable for most values of n (though it could have limitations towards the centre of the Sun)
print(np.median(mu_ssm))
mu = np.median(mu_ssm)
0.6825
#Input variables into the T equation stated above
T_0 = (m_H*mu*alpha_0**2*4*np.pi*G*rho_c_0*theta_0)/(k*(0+1))
T_1 = (m_H*mu*(alpha_1**2)*4*np.pi*G*rho_c_1*theta_1)/(k*(1+1))
T_2 = (m_H*mu*(alpha_2**2)*4*np.pi*G*rho_c_2*theta_2)/(k*(2+1))
T_3 = (m_H*mu*(alpha_3**2)*4*np.pi*G*rho_c_3*theta_3)/(k*(3+1))
T_4 = (m_H*mu*(alpha_4**2)*4*np.pi*G*rho_c_4*theta_4)/(k*(4+1))
T_5 = (m_H*mu*(alpha_5**2)*4*np.pi*G*rho_c_5*theta_5)/(k*(5+1))

#log values of T to make for a more readable plot
logT_0 = np.log10(T_0)
logT_1 = np.log10(T_1)
logT_2 = np.log10(T_2)
logT_3 = np.log10(T_3)
logT_4 = np.log10(T_4)
logT_5 = np.log10(T_5)
/ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:10: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in log10 /ext/sage/sage-8.2_1604/local/lib/python2.7/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in log10
#plot log of temperature against radius
fig, axis = plt.subplots(figsize = (9,5))
axis.plot(r_0, logT_0, label = 'n=0')
axis.plot(r_1, logT_1, label = 'n=1')
axis.plot(r_2, logT_2, label = 'n=2')
axis.plot(r_3, logT_3, label = 'n=3')
axis.plot(r_4, logT_4, label = 'n=4')
axis.plot(r_5, logT_5, label = 'n=5')
axis.set_ylabel('logT (K)')
axis.set_xlabel('r/R_o')
axis.legend(loc = 'lower left')
plt.show()

Plotting

#plot each variable against radius for each n
#plot for n=0
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_0, log_rho_0, label = 'n=0')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_0, M_0, label = 'n=0')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_0, logP_0, label = 'n=0')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_0, logT_0, label = 'n=0')
ax4.plot(r_ssm, log_T_ssm, 'k--', label = 'SSM')
ax4.set_ylabel('logT (K)')
ax4.set_xlabel('r/R_o')
ax4.set_ylim(4,)

ax1.legend(framealpha = 0.1, borderpad = 0.1)
f.subplots_adjust(hspace=0.1)
plt.show()
#plot each variable against radius for n=1
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_1, log_rho_1, 'g', label = 'n=1')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_1, M_1, 'g', label = 'n=1')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_1, logP_1, 'g', label = 'n=1')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_1, logT_1, 'g', label = 'n=1')
ax4.plot(r_ssm, log_T_ssm, 'k--', label = 'SSM')
ax4.set_ylabel('logT (K)')
ax4.set_xlabel('r/R_o')
ax4.set_ylim(4,)

ax1.legend(framealpha = 0.1, borderpad = 0.1)
f.subplots_adjust(hspace=0.1)
plt.show()
#plot the each variable against radius for n=2
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_2, log_rho_2, 'r', label = 'n=2')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_2, M_2, 'r', label = 'n=2')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_2, logP_2, 'r', label = 'n=2')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_2, logT_2, 'r', label = 'n=2')
ax4.plot(r_ssm, log_T_ssm, 'k--', label = 'SSM')
ax4.set_ylabel('logT (K)')
ax4.set_xlabel('r/R_o')
ax4.set_ylim(4,)

ax1.legend(framealpha = 0.1, borderpad = 0.1)
f.subplots_adjust(hspace=0.1)
plt.show()