import numpy as np
import matplotlib.pyplot as plt

c=3e10
x=0.73
y=0.25
z=0.02
a=7.565e-15
gamma=5/3
h=1.67e-24
k=1.3806e-16
G=6.67e-8
pi=np.pi

rho = np.zeros([101])
rho[0] = 3e-7

M = np.zeros([101])
M[0]=1.99e33

R = 6.96e10

L = np.zeros([101])
L[0] = 3.82e33

T = np.zeros([101])
T[0] = 6000.

pressure = np.zeros([101])
pressure[0] = rho[0]*((2*x)+((3/4)*y)+(z/2))*k*T[0]/h

shellthick = R/100
r = np.zeros([101])
r[0] = R


for o in range(1,100):
    r[o] = R - (o*shellthick)
    dM = 4 * pi* (r[o]**2) * rho[o-1] * shellthick
    dpressure = ((-1*rho[o-1]) * G * M[o-1] * shellthick) / (r[o]**2)
    epsilon = 2.5e6 * rho[o-1] * (0.73**2) * ((100000/T[o-1])**(2/3)) * np.exp(-33.8*((100000/T[o-1])**(1/3)))
    dL = 4 * pi * (r[o]**2) * rho[o-1] * shellthick * epsilon
    kappa = 4.34e25*1/(3.16*z)*(z*((1+x)))*(rho[o-1])/((T[o-1])**3.5)
    tcon = (1-(1/gamma)) * (T[o-1] / pressure[o-1]) * dpressure
    trad = (-3 * kappa * rho[o-1] *L[o-1] * shellthick) / (16 * pi * a * c * ((T[o-1])**3) * (r[o]**2))
    L[o] = L[o-1] - dL
    M[o] = M[o-1] - dM
    pressure[o]=pressure[o-1] - dpressure

    if tcon > trad:
        T[o] = T[o-1] + trad
    else:
        T[o] = T[o-1] + tcon

    rho[o] = rho[o-1] + (pressure[o]*h)/((k*T[o])*((2*x)+((3/4)*y)+(z/2)))