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


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

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


/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:30: RuntimeWarning: invalid value encountered in log10 /projects/sage/sage-7.3/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)

/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:15: 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
#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.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)

/projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:10: RuntimeWarning: invalid value encountered in log10 /projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in log10 /projects/sage/sage-7.3/local/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in log10 /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
#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)
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)
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)
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)
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)
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)
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.62240137 surface density n=0: 1409.62240137 central mass n=0: 2.08519689092e-10 surface mass n=0: 1.0 central pressure n=0: 1.34551543657e+14 surface pressure n=0: -2219201009.02 central temperature n=0: 7896785.15698 surface temperature n=0: -130.244166006
#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.39347987 surface density n=1: -0.939160449102 central mass n=1: 3.24960910989e-10 surface mass n=1: 1.0 central pressure n=1: 4.4284412797e+14 surface pressure n=1: 18139314.2312 central temperature n=1: 7895152.32134 surface temperature n=1: -1597.88492765
#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.9080448 surface density n=2: 6.95224384854e-05 central mass n=2: 4.23399403964e-10 surface mass n=2: 1.0 central pressure n=2: 1.8487526848e+15 surface pressure n=2: -524.945863461 central temperature n=2: 9504050.13864 surface temperature n=2: -624.674750696
#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.074105 surface density n=3: -9.59343887338e-10 central mass n=3: 5.05759441217e-10 surface mass n=3: 1.0 central pressure n=3: 1.24901421393e+16 surface pressure n=3: 0.00363515087498 central temperature n=3: 13496581.8811 surface temperature n=3: -313.481963864
#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.54023 surface density n=4: 3.95630100227e-19 central mass n=4: 5.67864345875e-10 surface mass n=4: 1.0 central pressure n=4: 2.80811301044e+17 surface pressure n=4: -1.03096599688e-13 central temperature n=4: 26342288.0514 surface temperature n=4: -21.5585449872
#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.92413 surface density n=5: 10.5822941172 central mass n=5: 5.95673813734e-10 surface mass n=5: 1.0 central pressure n=5: 8.15316522792e+17 surface pressure n=5: 339536347570.0 central temperature n=5: 30716910.8572 surface temperature n=5: 2654427.54591
#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.946 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