stellarinteriors.ipynb

Kernel: Python 2 (SageMath)

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: $$\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:

$$\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 = \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: $$\frac{3M_o}{4\pi r_o^3} = -3\rho_c|\frac{1}{\xi}\frac{d\theta}{d\xi}|$$

$$\rho_c = -\frac{M_o\xi}{4\pi r_o^3}\frac{d\xi}{d\theta}$$

Density: $$\rho = \rho_c\theta^n$$

Log Density: $$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.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:30: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/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\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\rho^{\frac{n+1}{n}} = K\rho_c^{\frac{n+1}{n}}\theta^{n+1}$$

K: $$\alpha = \sqrt{\frac{(n+1)K}{4\pi G\rho_c^{\frac{(n-1)}{n}}}}$$ $$K = \frac{\alpha^{2}4\pi G\rho_c^{\frac{(n-1)}{n}}}{(n+1)}$$

Combining these 2 equations: $$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.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:15: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/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 = \frac{kT\rho}{m_H\mu K} = K\rho_c^{\frac{(n+1)}{n}}\theta^{n+1}$$

Rearrange and substitute with density equation.

$$T = \frac{m_H\mu K \rho_c^{\frac{1}{n}}\theta}{k}$$

Input the K equation:

$$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.6825000000000001

#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.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:10: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/local/lib/python2.7/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in log10
/ext/sage/sage-8.6_1804/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()

#plot each variable against radius for n=3
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_3, log_rho_3, 'c', label = 'n=3')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_3, M_3, 'c', label = 'n=3')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_3, logP_3, 'c', label = 'n=3')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_3, logT_3, 'c', label = 'n=3')
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=4
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_4, log_rho_4, 'm', label = 'n=4')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_4, M_4, 'm', label = 'n=4')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_4, logP_4, 'm', label = 'n=4')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_4, logT_4, 'm', label = 'n=4')
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=5
f, (ax1, ax2, ax3, ax4) = plt.subplots(4, sharex=True, figsize = (9,9))
ax1.plot(r_5, log_rho_5, 'y', label = 'n=5')
ax1.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
ax1.set_ylabel('log(rho) (kgm^-3)')
ax1.set_ylim(-1,)
ax2.plot(r_5, M_5, 'y', label = 'n=5')
ax2.plot(r_ssm, m_ssm, 'k--', label = 'SSM')
ax2.set_ylabel('M (M_o)')
ax3.plot(r_5, logP_5, 'y', label = 'n=5')
ax3.plot(r_ssm, log_P_ssm, 'k--', label ='SSM')
ax3.set_ylabel('log(P) (Nm^-2)')
ax3.set_ylim(5,)
ax4.plot(r_5, logT_5, 'y', label = 'n=5')
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()

From the above figures it can be seen that the n=3 polytrope is the best fit to the standard solar model, more analysis of the data can be found in the report.

#zoom in the plot for density for n=2 & 3 polytropes, to show clear relationship
fig, axis = plt.subplots(figsize = (7,4))
axis.plot(r_2, log_rho_2, label = 'n=2', color='g')
axis.plot(r_3, log_rho_3, label = 'n=3', color='c')
axis.plot(r_ssm, log_rho_ssm, 'k--', label = 'SSM')
axis.set_ylim(0, 4)
axis.set_xlim(0.6,)
axis.set_ylabel('log(rho) (kgm^-3)')
axis.set_xlabel('r/R_o')
axis.legend()
plt.show()

#zoom in the plot for pressure for n=2 & 3 polytropes, to show clear relationship
fig, axis = plt.subplots(figsize = (7,4))
axis.plot(r_2, logP_2, label = 'n=2', color='g')
axis.plot(r_3, logP_3, label = 'n=3', color='c')
axis.plot(r_ssm, log_P_ssm, 'k--', label = 'SSM')
axis.set_ylim()
axis.set_xlim(0.6,)
axis.set_ylabel('log(P) (kgm^-2)')
axis.set_xlabel('r/R_o')
axis.legend()
plt.show()

#find values of each variable at the centre and surface of the Sun using slicing for n=0
print'central density n=0:', rho_0[0]
print'surface density n=0:', rho_0[-1]
print'central mass n=0:', M_0[0]
print'surface mass n=0:', M_0[-1]
print'central pressure n=0:', P_0[0]
print'surface pressure n=0:', P_0[-1]
print'central temperature n=0:', T_0[0]
print'surface temperature n=0:', T_0[-1]

central density n=0: 1409.6224013661824
surface density n=0: 1409.6224013661824
central mass n=0: 2.085196890921245e-10
surface mass n=0: 0.9999999999999996
central pressure n=0: 134551543657213.67
surface pressure n=0: -2219201009.021052
central temperature n=0: 7896785.156978242
surface temperature n=0: -130.24416600551606

#find values of each variable at the centre and surface of the Sun using slicing for n=1
print'central density n=1:', rho_1[0]
print'surface density n=1:', rho_1[-1]
print'central mass n=1:', M_1[0]
print'surface mass n=1:', M_1[-1]
print'central pressure n=1:', P_1[0]
print'surface pressure n=1:', P_1[-1]
print'central temperature n=1:', T_1[0]
print'surface temperature n=1:', T_1[-1]

central density n=1: 4640.393479871425
surface density n=1: -0.9391604491021224
central mass n=1: 3.2496091098874176e-10
surface mass n=1: 1.0
central pressure n=1: 442844127970057.06
surface pressure n=1: 18139314.23118434
central temperature n=1: 7895152.321344089
surface temperature n=1: -1597.884927648124

#find values of each variable at the centre and surface of the Sun using slicing for n=2
print'central density n=2:', rho_2[0]
print'surface density n=2:', rho_2[-1]
print'central mass n=2:', M_2[0]
print'surface mass n=2:', M_2[-1]
print'central pressure n=2:', P_2[0]
print'surface pressure n=2:', P_2[-1]
print'central temperature n=2:', T_2[0]
print'surface temperature n=2:', T_2[-1]

central density n=2: 16092.908044829437
surface density n=2: 6.952243848535457e-05
central mass n=2: 4.2339940396405507e-10
surface mass n=2: 1.0000000000000002
central pressure n=2: 1848752684800063.0
surface pressure n=2: -524.9458634614224
central temperature n=2: 9504050.138644136
surface temperature n=2: -624.6747506957851

#find values of each variable at the centre and surface of the Sun using slicing for n=3
print'central density n=3:', rho_3[0]
print'surface density n=3:', rho_3[-1]
print'central mass n=3:', M_3[0]
print'surface mass n=3:', M_3[-1]
print'central pressure n=3:', P_3[0]
print'surface pressure n=3:', P_3[-1]
print'central temperature n=3:', T_3[0]
print'surface temperature n=3:', T_3[-1]

central density n=3: 76561.0741050157
surface density n=3: -9.593438873378746e-10
central mass n=3: 5.057594412168205e-10
surface mass n=3: 1.0000000000000002
central pressure n=3: 1.2490142139297048e+16
surface pressure n=3: 0.003635150874981288
central temperature n=3: 13496581.881126745
surface temperature n=3: -313.48196386374104

#find values of each variable at the centre and surface of the Sun using slicing for n=4
print'central density n=4:', rho_4[0]
print'surface density n=4:', rho_4[-1]
print'central mass n=4:', M_4[0]
print'surface mass n=4:', M_4[-1]
print'central pressure n=4:', P_4[0]
print'surface pressure n=4:', P_4[-1]
print'central temperature n=4:', T_4[0]
print'surface temperature n=4:', T_4[-1]

central density n=4: 881912.540230316
surface density n=4: 3.95630100226574e-19
central mass n=4: 5.678643458752823e-10
surface mass n=4: 1.0000000000000002
central pressure n=4: 2.808113010443602e+17
surface pressure n=4: -1.0309659968824032e-13
central temperature n=4: 26342288.05137631
surface temperature n=4: -21.558544987223293

#find values of each variable at the centre and surface of the Sun using slicing for n=5
print'central density n=5:', rho_5[0]
print'surface density n=5:', rho_5[-1]
print'central mass n=5:', M_5[0]
print'surface mass n=5:', M_5[-1]
print'central pressure n=5:', P_5[0]
print'surface pressure n=5:', P_5[-1]
print'central temperature n=5:', T_5[0]
print'surface temperature n=5:', T_5[-1]

central density n=5: 2195902.9241341255
surface density n=5: 10.582294117181593
central mass n=5: 5.956738137340104e-10
surface mass n=5: 0.9999999999999999
central pressure n=5: 8.153165227918121e+17
surface pressure n=5: 339536347569.56964
central temperature n=5: 30716910.857173055
surface temperature n=5: 2654427.545914755

#find values of each variable at the centre and surface of the Sun using slicing for the SSM
print'central density SSM:', rho_ssm[0]
print'surface density SSM:', rho_ssm[-1]
print'central mass SSM:', m_ssm[0]
print'surface mass SSM:', m_ssm[-1]
print'central pressure SSM:', P_ssm[0]
print'surface pressure SSM:', P_ssm[-1]
print'central temperature SSM:', T_ssm[0]
print'surface temperature SSM:', T_ssm[-1]

central density SSM: 151000.0
surface density SSM: 0.9460000000000001
central mass SSM: 4e-07
surface mass SSM: 0.9999897
central pressure SSM: 2.34e+16
surface pressure SSM: 970000000.0
central temperature SSM: 15500000.0
surface temperature SSM: 81200.0