star_evollution_analysis.py

# Stellar Evolution Model With Neural Networks - Yan Barros

#### Reference: http://spiff.rit.edu/classes/phys370/lectures/statstar/statstar.html

#### The repository was used as the main base for generating data, however, i had to choose initial conditions for stars, based on Milky Way patterns. This is stored in a csv file which I used to feed the package tool and run a loop to generate my data for the project.

# Step 1 - Importing packages
"""

import os
import csv
import numpy as np
import pandas as pd
import seaborn as sns
import tensorflow as tf
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
from statsmodels.graphics.tsaplots import plot_acf

"""# Step 2 - Adaptating the main functions of the StatStar package"""

#  This program is a adaptation of a Python version of the code STATSTAR written
#  in FORTRAN and described in Carroll and Ostlie.  Program
#  has been modified and converted to Python in Fall 2010
#  by Prof. D. Schiminovich at Columbia University
#  for use in C3101 Stellar Structure and Evolution
#
#  Notes from FORTRAN code:
#
#  This program is listed in Appendix H of both "An Introduction
#  to Modern Astrophysics," Bradley W. Carroll and Dale A. Ostlie,
#  Addison-Wesley Publishing Company, copyright 1996, and "An
#  Introduction to Modern Stellar Astrophysics," Dale A. Ostlie
#  and Bradley W. Carroll, Addison-Wesley Publishing Company,
#  copyright 1996.
#
#  This program will calculate a static stellar model using the
#  equations developed in the text.  The user is expected to supply the
#  star's mass, luminosity, effective temperature, and composition
#  (X and Z).  If the choices for these quantities are not consistent
#  with the central boundary conditions, an error message will be
#  generated and the user will then need to supply a different set of
#  initial values.

#  SUBROUTINES

#  Subroutine STARTMDL computes values of M_r, L_r, P, and T, near the
#  surface of the star using expansions of the stellar structure
#  equations (M_r and L_r are assumed to be constant).
#
def STARTMDL(deltar, X, Z, mu, Rs, r_i, M_ri, L_ri, tog_bf,irc,cst):
#
#     X = hydrogen mass fraction
#     Z = metal mass fraction
#     mu = mean molecular weight
#     Rs = radius of star
#     r_i = radius of shell
#     M_ri = mass enclosed within shell at this radius
#          (assumed constant at start of model)
#     L_ri = enclosed luminosity, assumed constant
#
      r = r_i + deltar
      M_rip1 = M_ri
      L_rip1 = L_ri
#
#  This is the radiative approximation (neglect radiation pressure
#  and electron scattering opacity)# see Prialnik Eq. 5.1, 5.3 and Sec. 3.7 or C&O Eqs. (H.1), (H.2), (9.19),
#  and (9.20).
#
      if (irc == 0):
          T_ip1 = cst.G*M_rip1*mu*cst.m_H/(4.25e0*cst.k_B)*(1.0e0/r - 1.0e0/Rs)
          A_bf = 4.34e25*Z*(1.0e0 + X)/tog_bf
          A_ff = 3.68e22*cst.g_ff*(1.0e0 - Z)*(1.0e0 + X)
          Afac = A_bf + A_ff
          P_ip1 = np.sqrt((1.0e0/4.25e0)*(16.0e0/3.0e0*np.pi*cst.a*cst.c)*(cst.G*M_rip1/L_rip1)*(cst.k_B/(Afac*mu*cst.m_H)))*T_ip1**4.25e0

#
#  This is the convective approximation# see Prialnik Sec 6.5, 6.6 or C&O Eqs. (H.3) and (10.75).
#
      else:
          T_ip1 = cst.G*M_rip1*mu*cst.m_H/cst.k_B*(1.0e0/r - 1.0e0/Rs)/cst.gamrat
          P_ip1 = cst.kPad*T_ip1**cst.gamrat

      return r,P_ip1, M_rip1, L_rip1, T_ip1

#
#  Subroutine EOS calculates the values of density, opacity, the
#  guillotine-to-gaunt factor ratio, and the energy generation rate at
#  the radius r.
#
def   EOS(X, Z, XCNO, mu, P, T,izone,cst):

#
#  Solve for density from the ideal gas law (remove radiation pressure)#
#  see Prialnik Eq. 5.5 or C&O Eq. (10.26).
#
      if ((T < 0.0e0) or (P < 0.0e0)):
          print(' Something is a little wrong here.')
          print(' You are asking me to deal with either a negative temperature')
          print(' or a negative pressure.  I am sorry but that is not in my')
          print(' contract! You will have to try again with different')
          print(' initial conditions.')
          print(' In case it helps, I detected the problem in zone ',izone)
          print(' with the following conditions:')
          print('T = ',T,' K')
          print('P = ',P,' dynes/cm**2')
          return (0.0, 0.0, 0.0,0.0,1)


      Prad = cst.a*T**4/3.0e0
      Pgas = P - Prad
      rho=(mu*cst.m_H/cst.k_B)*(Pgas/T)
      if (rho < 0.0e0):
          print(' I am sorry, but a negative density was detected.')
          print(' my equation-of-state routine is a bit baffled by this new')
          print(' physical system you have created.  The radiation pressure')
          print(' is probably too great, implying that the star is unstable.')
          print(' Please try something a little less radical next time.')
          print(' In case it helps, I detected the problem in zone ')
          print(' with the following conditions:')
          print('T       = {0:12.5E} K'.format(T))
          print('P_total = {0:12.5E} dynes/cm**2'.format(P))
          print('P_rad   = {0:12.5E} dynes/cm**2'.format(Prad))
          print('P_gas   = {0:12.5E} dynes/cm**2'.format(Pgas))
          print('rho     = {0:12.5E} g/cm**3'.format(rho))
          return (rho, 0.0 , 0.0 ,0.0,1)

#
#  Calculate opacity, including the guillotine-to-gaunt factor ratio#
#  see Novotny (1973), p. 469. k_bf, k_ff, and k_e are the bound-free,
#  free-free, and electron scattering opacities, given by Prialnik Sec 3.7 or C&O Eqs. (9.19),
#  (9.20), and (9.21), respectively.
#

      tog_bf = 2.82e0*(rho*(1.0e0 + X))**0.2e0
      k_bf = 4.34e25/tog_bf*Z*(1.0e0 + X)*rho/T**3.5e0
      k_ff = 3.68e22*cst.g_ff*(1.0e0 - Z)*(1.0e0 + X)*rho/T**3.5e0
      k_e = 0.2e0*(1.0e0 + X)
      kappa = k_bf + k_ff + k_e
#
#  Compute energy generation by the pp chain and the CNO cycle.  These
#  are calculated using Eqs. (10.49) and (10.53) Prialnik Eq. 4.6, 4.7 + screening factor, which come from
#  Fowler, Caughlan, and Zimmerman (1975). The screening factor for the
#  pp chain is calculated as fpp# see Clayton (1968), p. 359ff.
#
#     Calculate 1/3 and 2/3 for convenience
      oneo3=0.333333333e0
      twoo3=0.666666667e0

      T6 = T*1.0e-06
      fx = 0.133e0*X*np.sqrt((3.0e0 + X)*rho)/T6**1.5e0
      fpp = 1.0e0 + fx*X
      psipp = 1.0e0 + 1.412e8*(1.0e0/X - 1.0e0)*np.exp(-49.98*T6**((-1.0)*oneo3))
      Cpp = 1.0e0 + 0.0123e0*T6**oneo3 + 0.0109e0*T6**twoo3 + 0.000938e0*T6
      epspp = 2.38e6*rho*X*X*fpp*psipp*Cpp*T6**(-twoo3)*np.exp(-33.80e0*T6**(-oneo3))
      CCNO = 1.0e0 + 0.0027e0*T6**oneo3 - 0.00778e0*T6**twoo3- 0.000149e0*T6
      epsCNO = 8.67e27*rho*X*XCNO*CCNO*T6**(-twoo3)*np.exp(-152.28e0*T6**(-oneo3))
      epslon = epspp + epsCNO

      return (rho, kappa, epslon,tog_bf,0)


#
#  The following four function subprograms calculate the gradients of
#  pressure, mass, luminosity, and temperature at r.

def dPdr(r, M_r, rho,cst):
      return -cst.G*rho*M_r/r**2

def dMdr(r, rho, cst):
      return (4.0e0*np.pi*rho*r**2)

def dLdr(r, rho, epslon,cst):
      return (4.0e0*np.pi*rho*epslon*r**2)

def dTdr(r, M_r, L_r, T, rho, kappa, mu, irc,cst):
      if (irc == 0):
          return (-(3.0e0/(16.0e0*np.pi*cst.a*cst.c))*kappa*rho/T**3*L_r/r**2)
#  This is the adiabatic convective temperature gradient (Prialnik Eq. 6.29 or C&O Eq. 10.81).
      else:
          return (-1.0e0/cst.gamrat*cst.G*M_r/r**2*mu*cst.m_H/cst.k_B)

#
# Runge-kutta algorithm
#
def RUNGE(f_im1, dfdr, r_im1, deltar, irc, X, Z, XCNO,mu, izone,cst):

      f_temp=np.zeros(4)
      f_i=np.zeros(4)

      dr12 = deltar/2.0e0
      dr16 = deltar/6.0e0
      r12  = r_im1 + dr12
      r_i  = r_im1 + deltar
#
#  Calculate intermediate derivatives from the four fundamental stellar
#  structure equations found in Subroutine FUNDEQ.
#
      for i in range(0,4):
          f_temp[i] = f_im1[i] + dr12*dfdr[i]

      df1, ierr = FUNDEQ(r12, f_temp, irc, X, Z, XCNO, mu, izone,cst)
      if (ierr != 0):
          return f_i,ierr
#

      for i in range(0,4):
          f_temp[i] = f_im1[i] + dr12*df1[i]

      df2,ierr = FUNDEQ(r12, f_temp, irc, X, Z, XCNO, mu, izone,cst)
      if (ierr != 0):
          return f_i,ierr

#
      for i in range(0,4):
          f_temp[i] = f_im1[i] + deltar*df2[i]

      df3,ierr=FUNDEQ(r_i, f_temp, irc, X, Z, XCNO, mu, izone,cst)
      if (ierr != 0):
          return f_i,ierr
#
#  Calculate the variables at the next shell (i + 1).
#
      for i in range(0,4):
          f_i[i] = f_im1[i] + dr16*(dfdr[i] + 2.0e0*df1[i] + 2.0e0*df2[i] + df3[i])

      return f_i,0


#
#  This subroutine returns the required derivatives for RUNGE, the
#  Runge-Kutta integration routine.
#
#      Subroutine FUNDEQ(r, f, irc, X, Z, XCNO, mu, izone,cst)

def FUNDEQ(r,f,irc,X,Z,XCNO,mu,izone,cst):

      dfdr=np.zeros(4)
      P   = f[0]
      M_r = f[1]
      L_r = f[2]
      T   = f[3]
      rho,kappa,epslon,tog_bf,ierr = EOS(X, Z, XCNO, mu, P, T, izone,cst)
      dfdr[0] = dPdr(r, M_r, rho,cst)
      dfdr[1] = dMdr(r, rho,cst)
      dfdr[2] = dLdr(r, rho, epslon,cst)
      dfdr[3] = dTdr(r, M_r, L_r, T, rho, kappa, mu, irc,cst)
      return (dfdr,ierr)



def StatStar(Msolar,Lsolar,Te,X,Z,iname):

#
#  Main program for calculating stellar structure
#
#  Variables, run-time parameters and settings:
#
#  deltar = radius integration step
#  idrflg = set size flag
#         = 0 (initial surface step size of Rs/1000.)
#         = 1 (standard step size of Rs/100.)
#         = 2 (core step size of Rs/5000.)
#
#  Nstart = number of steps for which starting equations are to be used
#           (the outermost zone is assumed to be radiative)
#  Nstop = maximum number of allowed zones in the star
#  Igoof = final model condition flag
#        = -1 (number of zones exceeded; also the initial value)
#        =  0 (good model)
#        =  1 (core density was extreme)
#        =  2 (core luminosity was extreme)
#        =  3 (extrapolated core temperature is too low)
#        =  4 (mass became negative before center was reached)
#        =  5 (luminosity became negative before center was reached)
#  X, Y, Z = mass fractions of hydrogen, helium, and metals
#  T0, P0 = surface temperature and pressure (T0 = P0 = 0 is assumed)
#  Ms, Ls, Rs = mass, luminosity, and radius of the star (cgs units)
#

      nsh=999   # maximum number of shells

      # initialize r, P, M_r, L_r, T, rho, kappa, epslon, dlPdlT to
      # an array with nsh elements

      r=np.zeros(nsh,float)
      P=np.zeros(nsh,float)
      M_r=np.zeros(nsh,float)
      L_r=np.zeros(nsh,float)
      T=np.zeros(nsh,float)
      rho=np.zeros(nsh,float)
      kappa=np.zeros(nsh,float)
      epslon=np.zeros(nsh,float)
      tog_bf=np.zeros(nsh,float)
      dlPdlT=np.zeros(nsh,float)

      # initialize other variables

      deltar=0.0
      XCNO=0.0
      mu=0.0
      Ms=0.0
      Ls=0.0
      Rs=0.0
      T0=0.0
      P0=0.0
      Pcore=0.0
      Tcore=0.0
      rhocor=0.0
      epscor=0.0
      rhomax=0.0
      Rsolar=0.0

      Y=1.0-(X+Z)


      #     tog_bf = bound-free opacity constant
      #             (ratio of guillotine to gaunt factors)
      tog_bf0=0.01


      # initialize variables used for 4 structure equation derivatives

      f_im1=np.zeros(4,float)
      dfdr=np.zeros(4,float)
      f_i=np.zeros(4,float)


      # Various run time parameters

      Nstart=10
      Nstop=999
      Igoof=-1
      ierr=0
      P0=0.0
      T0=0.0
      dlPlim=99.9
      debug=0

      #
      #  Assign values to constants (cgs units)
      #     Rsun = radius of the Sun
      #     Msun = mass of the Sun
      #     Lsun = luminosity of the Sun
      #     sigma = Stefan-Boltzmann constant
      #     c = speed of light in vacuum
      #     a = 4*sigma/c (radiation pressure constant)
      #     G = universal gravitational constant
      #     K_B = Boltzmann constant
      #     m_H = mass of hydrogen atom
      #     gamma = 5/3 (adiabatic gamma for a monatomic gas)
      #     gamrat = gamma/(gamma-1)
      #     kPad = P/T**(gamma/(gamma-1)) (adiabatic gas law constant)
      #     g_ff = free-free opacity gaunt factor (assumed to be unity)
      #

      Rsun=6.9599e10
      Msun=1.989e33
      Lsun=3.826e33

      # For convenience we will store most constants in a single variable

      class Constants:
           pass
      cst = Constants()

      cst.sigma=5.67051e-5
      cst.c=2.99792458e10
      cst.a=7.56591e-15
      cst.G=6.67259e-8
      cst.k_B=1.380658e-16
      cst.m_H=1.673534e-24
      cst.gamma= 5.0e0/3
      cst.g_ff= 1.0e0

#
#  Select the mass fraction CNO to be 50% of Z.
#
      XCNO = Z/2.0e0
#
#  Calculate the mass, luminosity, and radius of the star.
#  The radius is calculated from Prialnik Eq. 1.3 (C&O Eq. 3.17).
#
      Ms = Msolar*Msun
      Ls = Lsolar*Lsun
      Rs = np.sqrt(Ls/(4.e0*np.pi*cst.sigma))/Te**2
      Rsolar = Rs/Rsun
#
#  Begin with a very small step size since surface conditions vary
#  rapidly.
#
      deltar = -Rs/1000.0e0
      idrflg = 0
#
#  Calculate mean molecular weight mu assuming complete ionization
#  (see Prialnik Eq. 3.18, 3.25 or C&O Eq. 10.21).
#
      mu = 1.0e0/(2.0*X + 0.75*Y + 0.5*Z)
#
#  Calculate the delimiter between adiabatic convection and radiation
#  (see Prialnik Eq. 6.28 or C&O Eq. 10.87).
#
      cst.gamrat = cst.gamma/(cst.gamma - 1.0e0)
#
#  Initialize values of r, P, M_r, L_r, T, rho, kappa, and epslon at
#  the surface.  The outermost zone is assumed to be zone 1.  The zone
#  number increases toward the center.
#

      initsh=0
      r[initsh]   = Rs
      M_r[initsh] = Ms
      L_r[initsh] = Ls
      T[initsh]   = T0
      P[initsh]   = P0
      tog_bf[initsh]   = tog_bf0


      if (P0 <= 0.0) or (T0 <= 0.0):
          rho[initsh]    = 0.0
          kappa[initsh]  = 0.0
          epslon[initsh] = 0.0
          tog_bf[initsh] = 0.01
      else:
          rho[initsh],kappa[initsh],epslon[initsh],tog_bf[initsh],ierr=EOS(X, Z, XCNO, mu, P[initsh], T[initsh], 0 ,cst)
          if ierr != 0:
              print ("we're stopping now")
              istop=0

#
#  Apply approximate surface solutions to begin the integration,
#  assuming radiation transport in the outermost zone (the do 20 loop).
#  irc = 0 for radiation, irc = 1 for convection.
#  Assume arbitrary initial values for kPad, and dlPdlT.
#  dlPdlT = dlnP/dlnT (see Prialnik Eq. 6.28 or C&O Eq. 10.87)
#
      cst.kPad = 0.3e0
      irc = 0
      dlPdlT[initsh] = 4.25e0
      for i in range(0,Nstart):
          ip1 = i + 1
          r[ip1],P[ip1],M_r[ip1],L_r[ip1],T[ip1]=STARTMDL(deltar, X, Z, mu, Rs, r[i], M_r[i], L_r[i], tog_bf[i], irc,cst)
          rho[ip1],kappa[ip1],epslon[ip1],tog_bf[ip1],ierr=EOS(X, Z, XCNO, mu, P[ip1], T[ip1], ip1 ,cst)

          if ierr != 0:
              print('Values from the previous zone are:')
              print('r/Rs      = {0:12.5E}'.format(r[i]/Rs))
              print('rho       = {0:12.5E}  g/cm**3'.format(rho[i]))
              print('M_r/Ms    = {0:12.5E}'.format(M_r[i]/Ms))
              print('kappa     = {0:12.5E}  cm**2/g'.format(kappa[i]))
              print('T         = {0:12.5E}  K'.format(T[i]))
              print('epsilon   = {0:12.5E}  ergs/g/s'.format(epslon[i]))
              print('P         = {0:12.5E}  dynes/cm**2'.format(P[i]))
              print('L_r/Ls    = {0:12.5E}'.format(L_r[i]/Ls))
              break
#
#  Determine whether convection will be operating in the next zone by
#  calculating dlnP/dlnT numerically between zones i and i+1 [ip1].
#  Update the adiabatic gas constant if necessary.
#
          if (i > initsh):
              dlPdlT[ip1] = np.log(P[ip1]/P[i])/np.log(T[ip1]/T[i])
          else:
              dlPdlT[ip1] = dlPdlT[i]

          if (dlPdlT[ip1] < cst.gamrat):
              irc = 1
          else:
              irc = 0
              cst.kPad = P[ip1]/T[ip1]**cst.gamrat

#
#  Test to see whether the surface assumption of constant mass is still
#  valid.
#
          deltaM = deltar*dMdr(r[ip1], rho[ip1],cst)
          M_r[ip1] = M_r[i] + deltaM
          if (np.abs(deltaM) > (0.001e0*Ms)):
              if (ip1 > 1):
                  ip1 = ip1 - 1
                  print(' The variation in mass has become larger than 0.001*Ms')
                  print(' leaving the approximation loop before Nstart was reached')
                  break
#
#
#  =============  MAIN LOOP ==========================
#
#  This is the main integration loop.  The assumptions of constant
#  interior mass and luminosity are no longer applied.
#
#  Initialize the Runge-Kutta routine by specifying zone i-1 quantities
#  and their derivatives.  Note that the pressure, mass, luminosity,
#  and temperature are stored in the memory locations f_im1(1),
#  f_im1(2), f_im1(3), and f_im1(4), respectively.  The derivatives of
#  those quantities with respect to radius are stored in dfdr(1),
#  dfdr(2), dfdr(3), and dfdr(4).  Finally, the resulting values for
#  P, M_r, L_r, and T are returned from the Runge-Kutta routine in
#  f_i(1), f_i(2), f_i(3), and f_i(4), respectively.
#
#  The stellar structure equations
#           dPdr (Prialnik 5.1 or C&O Eq. 10.7)
#           dMdr (Prialnik 5.2 or C&O Eq. 10.8)
#           dLdr (Prialnik 5.4 or C&O Eq. 10.45) and
#           dTdr (Prialnik 5.3 or C&O Eq. 10.61 or Eq. 10.81)
#           are calculated in function calls, defined previously in the code.
#

      Nsrtp1 = ip1 + 1

      if (ierr != 0):    # exit if we've arrived at this point after an error in initialization
          Nstop=Nsrtp1-1
          istop=Nstop

      for i in range(Nsrtp1,Nstop):
          im1 = i - 1
          f_im1[0] = P[im1]
          f_im1[1] = M_r[im1]
          f_im1[2] = L_r[im1]
          f_im1[3] = T[im1]
          dfdr[0]  = dPdr(r[im1], M_r[im1], rho[im1],cst)
          dfdr[1]  = dMdr(r[im1], rho[im1],cst)
          dfdr[2]  = dLdr(r[im1], rho[im1], epslon[im1],cst)
          dfdr[3]  = dTdr(r[im1], M_r[im1], L_r[im1], T[im1], rho[im1],kappa[im1], mu, irc,cst)
          f_i,ierr=RUNGE(f_im1, dfdr, r[im1], deltar, irc, X, Z, XCNO, mu, i,cst)

          if (ierr != 0):
              print(' The problem occurred in the Runge-Kutta routine')
              print(' Values from the previous zone are:')
              print('r/Rs    = {0:12.5e}'.format(r[im1]/Rs))
              print('rho     = {0:12.5e} g/cm**3'.format(rho[im1]))
              print('M_r/Ms  = {0:12.5e}'.format(M_r[im1]/Ms))
              print('kappa   = {0:12.5e} cm**2/g'.format(kappa[im1]))
              print('T       = {0:12.5e} K'.format(T[im1]))
              print('epsilon = {0:12.5e} ergs/g/s'.format(epslon[im1]))
              print('P       = {0:12.5e} dynes/cm**2'.format(P[im1]))
              print('L_r/Ls  = {0:12.5e}'.format(L_r[im1]/Ls))
              break
#
#  Update stellar parameters for the next zone, including adding
#  dr to the old radius (note that dr <  0 since the integration is
#  inward).
#
          r[i]   = r[im1] + deltar
          P[i]   = f_i[0]
          M_r[i] = f_i[1]
          L_r[i] = f_i[2]
          T[i]   = f_i[3]
#
#  Calculate the density, opacity, and energy generation rate for
#  this zone.
#
          rho[i],kappa[i],epslon[i],tog_bf[i],ierr=EOS(X, Z, XCNO, mu, P[i], T[i], i, cst)

          if (ierr != 0):
              print(' Values from the previous zone are:')
              print('r/Rs    = {0:12.5e}'.format(r[im1]/Rs))
              print('rho     = {0:12.5e} g/cm**3'.format(rho[im1]))
              print('M_r/Ms  = {0:12.5e}'.format(M_r[im1]/Ms))
              print('kappa   = {0:12.5e} cm**2/g'.format(kappa[im1]))
              print('T       = {0:12.5e} K'.format(T[im1]))
              print('epsilon = {0:12.5e} ergs/g/s'.format(epslon[im1]))
              print('P       = {0:12.5e} dynes/cm**2'.format(P[im1]))
              print('L_r/Ls  = {0:12.5e}'.format(L_r[im1]/Ls))
              istop = i
              break

          if (debug == 1): print (i,r[i],M_r[i],L_r[i],T[i],P[i],rho[i],kappa[i],epslon[i],tog_bf[i])



#
#  Determine whether convection will be operating in the next zone by
#  calculating dlnP/dlnT and comparing it to gamma/(gamma-1)
#  (see Prialnik Eq. 6.28 or C&O Eq. 10.87).  Set the convection flag appropriately.
#
          dlPdlT[i] = np.log(P[i]/P[im1])/np.log(T[i]/T[im1])
          if (dlPdlT[i] < cst.gamrat):
              irc = 1
          else:
              irc = 0

#
#  Check to see whether the center has been reached.  If so, set Igoof and
#  estimate the central conditions rhocor, epscor, Pcore, and Tcore.
#  The central density is estimated to be the average density of the
#  remaining central ball, the central pressure is determined by
#  using Taylor expansion at center (Prialnik - Exercise. 5.1; CO Eq. H.4)
#   and the central value for the energy
#  generation rate is calculated to be the remaining interior
#  luminosity divided by the mass of the central ball.  Finally, the
#  central temperature is computed by applying the ideal gas law
#  (where radiation pressure is neglected).
#
          if ((r[i] <= np.abs(deltar)) and ((L_r[i] >= (0.1e0*Ls)) or (M_r[i] >= (0.01e0*Ms)))):
              #   Hit center before mass/luminosity depleted
              Igoof = 6
          elif (L_r[i] <= 0.0e0):
              #   Obtained negative central luminosity
              Igoof = 5
              # rho: Taylor expansion at center
              rhocor = M_r[i]/(4.0e0/3.0e0*np.pi*r[i]**3)
              if (M_r[i] != 0.0e0):
                  # energy generation rate
                  epscor = L_r[i]/M_r[i]
              else:
                  epscor = 0.0e0
              # P:  Taylor expansion at center
              Pcore = P[i] + 2.0e0/3.0e0*np.pi*cst.G*rhocor**2*r[i]**2
              # Assume ideal gas
              Tcore = Pcore*mu*cst.m_H/(rhocor*cst.k_B)
          elif (M_r[i] <= 0.0e0):
              Igoof  = 4  #  Model has a hole in center (negative density!)
              Rhocor = 0.0e0
              epscor = 0.0e0
              Pcore  = 0.0e0
              Tcore  = 0.0e0
          elif ((r[i] < (0.02e0*Rs)) and ((M_r[i] < (0.01e0*Ms)) and ((L_r[i] < 0.1e0*Ls)))):
              #  if we've reached <2% star's radius,
              #<1% mass enclosed and <10% luminosity then....
              # rho: Taylor expansion at center
              rhocor = M_r[i]/(4./3.*np.pi*r[i]**3)
              # set maximum reasonable core mass
              rhomax = 10.0e0*(rho[i]/rho[im1])*rho[i]
              epscor = L_r[i]/M_r[i]
              # P: Taylor expansion at center
              Pcore  = P[i] + 2.0e0/3.0e0*np.pi*cst.G*rhocor**2*r[i]**2
              # Assume ideal gas
              Tcore  = Pcore*mu*cst.m_H/(rhocor*cst.k_B)
            # In general, these should all produce values
              # that rise towards center (but not too high)
              if ((rhocor < rho[i]) or (rhocor > rhomax)):
                  # rho is off either large or small
                  Igoof = 1
              elif (epscor < epslon[i]):
                  # energy generation rate a bit off (low)
                  Igoof = 2
              elif (Tcore < T[i]):
                  # Temperature a bit off (low)
                  Igoof = 3
              else:
                  # number of allowed shells has been exceeded
                  Igoof = 0

          if (Igoof != -1):
              istop = i
              break
#
#  Is it time to change the step size?
#
          if ((idrflg == 0) and (M_r[i] < (0.99e0*Ms))):
               deltar = (-1.0)*Rs/100.0e0
               idrflg = 1

          if ((idrflg == 1) and (deltar >= (0.5*r[i]))):
               deltar = (-1.0)*Rs/5000.0e0
               idrflg = 2

          istop = i

#
#  Generate warning messages for the central conditions.
#
      rhocor = M_r[istop]/(4.0e0/3.0e0*np.pi*r[istop]**3)
      epscor = L_r[istop]/M_r[istop]
      Pcore  = P[istop] + 2.0e0/3.0e0*np.pi*cst.G*rhocor**2*r[istop]**2
      Tcore  = Pcore*mu*cst.m_H/(rhocor*cst.k_B)

      if  (Igoof != 0):
          if (Igoof == -1):
              print('Sorry to be the bearer of bad news, but...')
              print('       Your model has some problems')
              print('The number of allowed shells has been exceeded')

          if (Igoof == 1):
              print('It looks like you are getting close,')
              print('however, there are still a few minor errors')
              print('The core density seems a bit off,')
              print(' density should increase smoothly toward the center.')
              print(' The density of the last zone calculated was rho = ',rho[istop],' gm/cm**3')
              print (rhocor,rhomax)
          if (rhocor > 1e10):
              print('It looks like you will need a degenerate')
              print(' neutron gas and general relativity')
              print(' to solve this core.  Who do you think I am, Einstein?')

          if (Igoof == 2):
              print('It looks like you are getting close,')
              print('however, there are still a few minor errors')
              print('The core epsilon seems a bit off,')
              print(' epsilon should vary smoothly near the center.')
              print(' The value calculated for the last zone was eps =',epslon[istop],' ergs/g/s')

          if (Igoof == 3):
              print('It looks like you are getting close,')
              print('however, there are still a few minor errors')
              print(' Your extrapolated central temperature is too low')
              print(' a little more fine tuning ought to do it.')
              print(' The value calculated for the last zone was T = ',T[istop],' K')

          if (Igoof == 4):
              print('Sorry to be the bearer of bad news, but...')
              print('       Your model has some problems')
              print('You created a star with a hole in the center!')

          if (Igoof == 5):
              print('Sorry to be the bearer of bad news, but...')
              print('       Your model has some problems')
              print('This star has a negative central luminosity!')

          if (Igoof == 6):
              print('Sorry to be the bearer of bad news, but...')
              print('       Your model has some problems')
              print('You hit the center before the mass and/or ')
              print('luminosity were depleted!')
      else:
          print('CONGRATULATIONS, I THINK YOU FOUND IT!')
          print('However, be sure to look at your model carefully.')

#
#  Print the central conditions.  If necessary, set limits for the
#  central radius, mass, and luminosity if necessary, to avoid format
#  field overflows.
#

      Rcrat = r[istop]/Rs
      if (Rcrat < -9.999e0): Rcrat = -9.999e0
      Mcrat = M_r[istop]/Ms
      if (Mcrat < -9.999e0): Mcrat = -9.999e0
      Lcrat = L_r[istop]/Ls
      if (Lcrat < -9.999e0): Lcrat = -9.999e0

      f=open(f'starmodl_py{iname}.dat','w')

      f.write('A Homogeneous Main-Sequence Model\n')
      f.write(' The surface conditions are:        The central conditions are:\n')
      f.write(' Mtot = {0:13.6E} Msun          Mc/Mtot     = {1:12.5E}\n'.format(Msolar,Mcrat))
      f.write(' Rtot = {0:13.6E} Rsun          Rc/Rtot     = {1:12.5E}\n'.format(Rsolar,Rcrat))
      f.write(' Ltot = {0:13.6E} Lsun          Lc/Ltot     = {1:12.5E}\n'.format(Lsolar,Lcrat))
      f.write(' Teff = {0:13.6E} K             Density     = {1:12.5E}\n'.format(Te,rhocor))
      f.write(' X    = {0:13.6E}               Temperature = {1:12.5E}\n'.format(X,Tcore))
      f.write(' Y    = {0:13.6E}               Pressure    = {1:12.5E} dynes/cm**2\n'.format(Y,Pcore))
      f.write(' Z    = {0:13.6E}               epsilon     = {1:12.5E} ergs/s/g\n'.format(Z,epscor))
      f.write('                                    dlnP/dlnT   = {0:12.5E}\n'.format(dlPdlT[istop]))

      f.write('Notes:\n')
      f.write(' (1) Mass is listed as Qm = 1.0 - M_r/Mtot, where Mtot = {0:13.6}\n'.format(Msun))
      f.write(' (2) Convective zones are indicated by c, radiative zones by r\n')
      f.write(' (3) dlnP/dlnT may be limited to +99.9 or -99.9# if so it is\n')
      f.write(' labeled by *\n')



#
#  Print data from the center of the star outward, labeling convective
#   or radiative zones by c or r, respectively.  If abs(dlnP/dlnT)
#  exceeds 99.9, set a print warning flag (*) and set the output limit
#  to +99.9 or -99.9 as appropriate to avoid format field overflows.
#
      f.write('   r        Qm       L_r       T        P        rho      kap      eps     dlPdlT\n')

      for ic in range(0,istop+1):
          i = istop - ic
          Qm = 1.0e0 - M_r[i]/Ms    # Total mass fraction down to radius


          if (dlPdlT[i] < cst.gamrat):
              rcf = 'c'
          else:
              rcf = 'r'
          if (np.abs(dlPdlT[i]) > dlPlim):
              dlPdlT[i] = np.copysign(dlPlim,dlPdlT[i])
              clim = '*'
          else:
              clim = ' '
          s='{0:7.2E} {1:7.2E} {2:7.2E} {3:7.2E} {4:7.2E} {5:7.2E} {6:7.2E} {7:6.2E}{8:1s}{9:1s} {10:5.1f}\n'.format(r[i], Qm, L_r[i], T[i], P[i], rho[i], kappa[i],epslon[i], clim, rcf, dlPdlT[i])
          f.write(s)

#     Output to screen
      print
      print('***** The integration has been completed *****')
      print('      The model has been stored in starmodl_py.dat')
      print
      return Igoof,ierr,istop, f'starmodl_py{iname}.dat'

"""# Step 2 - Load the Initial Conditions dataset to extract data with the StatStar Model."""

df=pd.read_csv('stars_data_gen.csv').values
models=[]
for i in range(0,len(df)):
  Msolar=df[i][2]
  Lsolar=df[i][3]
  Te=df[i][4]
  X=df[i][5]
  Z=df[i][6]
  Y = 1.e0 - X - Z
  if Y < 0:
    continue

  Igoof,ierr,istop,model=StatStar(Msolar,Lsolar,Te,X,Z,i)
  models.append([model,df[i][0],df[i][1]])

"""# Checking the .dat Files"""

path='/content'
files=os.listdir(path)
files

"""# Step 3 - Converting the outputs from .dat file's to .csv files."""

i=1
for file in files[1:]:
  try:
    datContent = [i.strip().split() for i in open('/content/'+file).readlines()]

    with open(f"/content/starmodl_{i}.csv", "w") as f:
        writer = csv.writer(f)
        writer.writerows(datContent)
  except:
    pass
  print(f'Iteration {i} of {len(files)-1}.')
  i+=1

"""# Step 4 - A Loop to Generate dataframes from .csv files and separating the important variables for the model."""

files=os.listdir(path)
for file in files:
  if file[-3:]=='csv':
    try:
      df=pd.read_csv(f'/content/{file}',skiprows=15)
      radius=df['r']
      temperature=df['T']
      density=df['rho']
      mass=df['Qm']
    except:
      pass

"""# Step 5 - Using the first Dataset as an example"""

df=pd.read_csv('/content/starmodl_1.csv',skiprows=15)
radius=df['r'].values
temperature=df['T'].values
density=df['rho'].values
mass=df['Qm'].values
star_age=df.reset_index()['index'].values

"""# Step 6 - Visualizing Data"""

variables=[radius,temperature,density,mass]
labels=['Radius of Star/Radius of Sun','Temperature of Star', 'Density of Star', 'Mass of Star/Mass of Sun']
sns.set(style="whitegrid")
i=0
for value in variables:
  plt.plot(star_age,value)
  # Add x and y-axis labels
  plt.xlabel("Star_Age")
  plt.ylabel(labels[i])
  plt.show()
  i+=1

"""# Step 7 - Normalizing Data"""

radius=radius-np.mean(radius) #exp
temperature=temperature-np.mean(temperature) #exp
density=density-np.mean(density) #exp
mass=mass-np.mean(mass) #exp

"""# Step 8 - Plotting Radius Function by Exponential fitting"""

plt.plot(radius)
plt.plot(np.exp(radius))

"""# Step 9 - Visualizing the Difference between the square of delta maximum values and the square of delta minimum values of the radius function."""

diff=((max(np.exp(radius))-max(radius))**2-(min(np.exp(radius))-min(radius))**2)
print(diff)

plt.plot(radius, label="Real Value for Radius/Sun Radius")
plt.plot(np.exp(radius)-diff, label="Predicted by exponential model")
sns.set(style="whitegrid")
plt.xlabel("Steps")
plt.legend()
plt.title("Comparison of Real and Predicted Values for Radius/Sun Radius")
plt.show()

"""# Step 10 - Creating and visualizing the Autocorrelation Function"""

def autocor(data):
    return acf(data, nlags=len(data) - 1)

fig=plt.figure(figsize=(20,10))
x1 = autocor(radius)
shifted_output = np.concatenate((np.zeros(3), x1[:-3]))
# Subtract the corresponding samples
auto_correlation = x1 - shifted_output
plt.plot(auto_correlation)
plt.title("Autocorrelation of values (Should be a Delta Function)")

fig, ax = plt.subplots(figsize=(8,4))
plot_acf(radius, ax=ax, lags=20)
ax.set_xlim(0.20)
ax.set_ylim(-0.2,1.0)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Sample Autocorrelation Function')
plt.grid(True)
plt.show()

fig=plt.figure(figsize=(20,10))
x1 = autocor(radius-(np.exp(radius)-diff))
shifted_output = np.concatenate((np.zeros(3), x1[:-3]))
auto_correlation = x1 - shifted_output
plt.title("Autocorrelation of Residual values (Should be a Delta Function)")
plt.plot(auto_correlation)

fig, ax = plt.subplots(figsize=(8,4))
plot_acf(radius-(np.exp(radius)-diff), ax=ax, lags=20)
ax.set_xlim(0.20)
ax.set_ylim(-0.2,1.0)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Sample Autocorrelation Function')
plt.grid(True)
plt.show()

"""# Step 11 - Visualizing Real and Model Functions"""

x1 = autocor(radius)
x2 = autocor(radius-(np.exp(radius)-diff))

plt.plot(x1[:200])
plt.plot(x2[:200])

"""# Step 12 - Developing a linear regression Algorithm for testing"""

X=radius
y_list=[]
x_list=[]
for i in range(3, len(X)):
  try:
    y_list.append([X[i-1],X[i-2],X[i-3]])
    x_list.append([X[i]])
  except Exception as e:
    print(e)
    pass
x=np.array(y_list)
y=np.array(x_list)

"""# Step 13 - Developing a linear Neural Network with 1 layer"""

def build_nn(input_shape):
    model = tf.keras.Sequential()
    model.add(tf.keras.layers.Dense(1, input_shape=input_shape, activation='linear'))

    return model

input_shape = (x.shape[1],)

model = build_nn(input_shape)
model.compile(optimizer='adam', loss='mse')

"""#Step 14 - Training the model with 10000 steps"""

history=model.fit(x, y, epochs=4000, batch_size=32)

"""# Step 15 - Visualizing the loss function"""

fig,ax=plt.subplots(figsize=(14,12))
plt.plot(history.history['loss'])
plt.xlabel('Epochs')
plt.ylabel('Loss Function')
plt.title('Loss Function per Epoch')
plt.grid(True)
plt.show()

"""# Step 16 - Applying Fourier Transform on Loss Function for valuating results"""

fig,ax=plt.subplots(figsize=(14,12))
plt.plot(np.fft.fft(history.history['loss'])[:4000])
plt.xlabel('Epochs')
plt.ylabel('Loss Function')
plt.title('Loss Function per Epoch')
plt.grid(True)
plt.show()

"""# Step 17 - Analyzing Autocorrelation"""

fig, ax = plt.subplots(figsize=(8,4))
plot_acf(np.fft.fft(history.history['loss'])[:4000], ax=ax, lags=20)
ax.set_ylim(-0.2,1.5)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Sample Autocorrelation Function')
plt.grid(True)
plt.show()

"""# Step 18 - Predicting Data and Visualizating Data"""

yHat=model.predict(x)

fig=plt.figure(figsize=(20,10))
plt.plot(y, label='Real output')
plt.scatter(range(len(yHat)), yHat, label='Predicted output', s=7.5)
sns.set(style="whitegrid")
plt.xlabel("Steps")
plt.legend()
plt.title("Comparison of Real and Predicted Values for Radius/Sun Radius using a 1 Layer Neural Network")
plt.show()

"""# Step 19 - Autocorrelation"""

fig=plt.figure(figsize=(20,10))
x1 = autocor(y)
shifted_output = np.concatenate((np.zeros(3), x1[:-3]))
auto_correlation = x1 - shifted_output
plt.plot(auto_correlation)

fig, ax = plt.subplots(figsize=(8,4))
plot_acf(auto_correlation, ax=ax, lags=20)
ax.set_xlim(0.20)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Sample Autocorrelation Function')
plt.grid(True)
plt.show()

fig=plt.figure(figsize=(20,10))
x1 = autocor(y-yHat)
shifted_output = np.concatenate((np.zeros(3), x1[:-3]))
auto_correlation = x1 - shifted_output
plt.plot(auto_correlation)

fig, ax = plt.subplots(figsize=(8,4))
plot_acf(auto_correlation, ax=ax, lags=20)
ax.set_xlim(0.20)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Sample Autocorrelation Function for Residual Data')
plt.grid(True)
plt.show()

"""# Step 20 - Saving Model"""

model.save('radius_star.h5')