#include <stdio.h>
#include <math.h>                  /* The necessary includes */
#include <float.h>
#include <conio.h>

#define PI 3.14159
#define G 6.67259e-8
#define SIGMA 7.56471e-15
#define H 1.67e-24                 /* The definitions of the constants */
#define K 1.38e-16
#define C 3.0e10


/*
 *
 * S T E L L A R   I N T E R I O R
 *
 * James Marshall
 *
 *
 *
 * Calculates the mass, luminosity, pressure, temperature, and density
 * for the solar interior using the appropriate equations.  To do so, it
 * also calculates the opacity and energy generated.  In order to make the
 * representation more acurrate, the fractions of H, He, and metals are
 * changed from surface to core values once the core region is reached.
 * The stellar mass must be kept within 3 solar masses or else the
 * approximating equations fail to hold true.  The user also has the option
 * of including the CNO cycle with the PP chain.
 *
 *
 *
 * These are the possible return values and when they occur:
 *
 * (-1) -- The output file could not be opened.
 * (0) --- The program ended successfully (reached stellar interior).
 * (1) --- The mass and/or luminosity hit zero before the radius did.
 * (3) --- The user entered a mass greater than 3 solar masses and then
 *         quit the program at the warning message.
 * (any other value) -- Program was quit for a reason other than the above.
 *
 */



double FinalCalc (double delta, double dvalue, double outside)
{
  double newvalue;
					     /* Q(r) = -((delta r * dQ/dr) */
  newvalue = -((dvalue * delta) - outside);  /*          - Q(r + delta r)) */
  return (newvalue);
}



double Mass (double mass, double rho, double r, double delta)
{
  double dmass, newmass;
						   /* dM/dr = rho 4 PI r^2 */
  dmass = (rho * 4 * PI * r * r);
  newmass = FinalCalc (delta, dmass, mass);
  return (newmass);
}



double Luminosity (double luminosity, double rho, double epsilon,
		   double r, double delta)
{
  double dlum, newlum;
					       /* dL/dr = eps rho 4 PI r^2 */
  dlum = (epsilon * rho * 4 * PI * r * r);
  newlum = FinalCalc (delta, dlum, luminosity);
  return (newlum);
}



double Pressure (double pressure, double mass, double rho, double r,
		 double delta)
{
  double dpres, newpres;
					       /* dP/dr = -(G M rho / r^2) */
  dpres = -((G * mass * rho) / (r * r));
  newpres = FinalCalc (delta, dpres, pressure);
  return (newpres);
}



double Temperature (double temp, double kappa, double rho, double lum,
		    double r, double delta, double gamma, double mass,
		    double pres)
{
  double dtemp1, dtemp2, newtemp1, newtemp2, finaltemp;

  dtemp1 = -((3 * kappa * rho * lum) /
	     (16 * PI * SIGMA * C * temp * temp * temp * r * r));
  newtemp1 = FinalCalc (delta, dtemp1, temp);
  dtemp2 = ((1 - (1/gamma)) * (temp / pres) * (-(G * mass * rho) / (r * r)));
  newtemp2 = FinalCalc (delta, dtemp2, temp);
  if (newtemp1 < newtemp2)
    finaltemp = newtemp1;                        /* newtemp1 and newtemp2  */
  else finaltemp = newtemp2;                     /* dTrad/dr and dTconv/dr */
  return (finaltemp);
}



double Density (double pres, double temp, double X, double Y, double Z)
{
  double newrho;

  newrho = ((pres * H) / (K * temp * ((2 * X) + (3 * Y / 4) + (Z / 2))));
  return (newrho);
}



double Kappa (double X, double Z, double rho, double temp)
{
  double newkappa, tempcalc;

  tempcalc = pow (sqrt (temp), 7.0);
  newkappa = ((4.34e25 * Z * (1 + X) * rho) /  (3.162 * tempcalc));
  return (newkappa);
}



double Epsilon (double rho, double X, double Z, double temp, int CNO)
{
  double eps_pp, eps_cno, newepsilon, exponent, powerterm;

  exponent = exp (-33.8 * pow ((1.0e6 / temp), 0.33333));
  powerterm = pow ((1.0e6 / temp), 0.66666);
  eps_pp = 2.5e6 * rho * X * X * powerterm * exponent;
  if (CNO == 1)
  {
    exponent = exp (-152.3 * pow ((1.0e6 / temp), 0.33333));
    powerterm = pow ((1.0e6 / temp), 0.66666);
    eps_cno = 9.5e28 * rho * X * (Z / 3) * powerterm * exponent;
    newepsilon = (eps_pp + eps_cno);
  }
  else
    newepsilon = eps_pp;
  return (newepsilon);
}



main (void)
{
  /* Set up variables */

  char filename [51];
  int count, steps, valuesOK, choice, dummy2, CNO, proceed;
  double mass, lum, pres, temp, rho, kappa, epsilon, radius, gamma,
	 newmass, newlum, newpres, newtemp, newrho, newkappa, newepsilon,
	 newradius, X, Y, Z, Xcore, Ycore, Zcore, delta, dummy;
  FILE *output;



  /* Clear out the filename string so it contains only spaces */

  for (count = 0; count <= 50; count++)
    filename [count] = ' ';



  /* Print out general information about program */

  printf ("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
  printf ("This program will allow you to make a model of a star.\n");
  printf ("The data obtained will be output to a file which can be used\n");
  printf ("to create graphs or for use in other calculations.  The model\n");
  printf ("is made by starting with initial values, which I will ask you\n");
  printf ("to input, and how many steps you want to take.  The steps are\n");
  printf ("the total number of small shells of the star at which I make\n");
  printf ("calculations as I move towards the center.  The star's mass\n");
  printf ("must be within 3 solar masses or else the approximations used\n");
  printf ("in the calculations no longer hold true.\n\n");
  printf ("<Press any key to continue.>");
  proceed = getch ();
  printf ("\n\n");                               /* Space between sections */


  /* Open file for output data of solar values at each step */

  printf ("Now I need the name of the file you want the data saved to.\n");
  printf ("Remember that the file name can be a maximum size of \n");
  printf ("up to 8 characters, a period, and up to 3 more characters.\n");
  printf ("A maximum length file name will look like ????????.???\n");
  printf ("but please enter the full path and not just a file name.\n");
  printf ("You can use up to 4 directories with names up to 8 characters\n");
  printf ("each before you must enter the name of the file.  If you do\n");
  printf ("not enter the file name in this way, either the data will be\n");
  printf ("saved somewhere else or the file won't be opened and the\n");
  printf ("program will be ended with an error message.\n");
  printf ("For example, a valid path is a:\\sun&star\\sunmodel.dat\n\n");
  printf ("<Press any key to enter path and file name.>");
  proceed = getch();
  printf ("\n\n");                                         /* Leave a space */
  printf ("Path and file name:  ");
  scanf ("%s51", filename);
  output = fopen (filename, "w");
  if (output == NULL)
  {
    printf ("Unable to open file %s\n", filename);
    printf ("The program is being terminated.\n");
    return (1);
  }
  printf ("\n\n");                               /* Space between sections */


  /* Print instructions for entering the necessary stellar data */

  printf ("Now I need the initial values (at the star's surface)\n");
  printf ("for the following variables.  I will also give you the proper\n");
  printf ("values to use for our own sun.  If you want to use the values\n");
  printf ("for the sun, you can just enter a zero instead of the number.\n");
  printf ("Otherwise, you must type in the new value.\n\n");
  printf ("<Press any key to begin entering the data.>");
  proceed = getch();
  printf ("\n\n");                                        /* Leave a space */



  /* Set up initial values for the sun */

  mass = 1.99e33;
  lum = 3.90e33;
  pres = 2.47e5;
  temp = 5780.0;
  rho = 3.00e-7;
  radius = 6.96e10;
  X = 0.710;
  Y = 0.270;
  Z = 0.020;
  Xcore = 0.330;
  Ycore = 0.630;
  Zcore = 0.040;
  gamma = 1.666;
  steps = 100;
  CNO = 0;


  /* Get all of the necessary information */

  printf ("Mass in g (%.2le, max. %.2le):  ", mass, (3 * mass));
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    mass = dummy;
  if (mass > 5.97e33)
  {
    printf ("The approximating calculations will not hold true for this mass.\n");
    printf ("Do you still want to continue the program? (0=N, 1=Y):  ");
    scanf ("%d", &choice);
    if (choice == 0)     /* Mass entered too big for accurate calculations */
      return (3);
  }
  printf ("Luminosity in erg/s (%.2le):  ", lum);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    lum = dummy;
  printf ("Pressure in dyne/cm^2 (%.2le):  ", pres);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    pres = dummy;
  printf ("Temperature in K (%.2le):  ", temp);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    temp = dummy;
  printf ("Density in g/cm^3 (%.2le):  ", rho);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    rho = dummy;
  printf ("Radius in cm (%.2le):  ", radius);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    radius = dummy;
  printf ("surface fraction of H (%.3lf):  ", X);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    X = dummy;
  printf ("surface fraction of He (%.3lf):  ", Y);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    Y = dummy;
  printf ("surface fraction of metals (%.3lf):  ", Z);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    Z = dummy;
  printf ("core fraction of H (%.3lf):  ", Xcore);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    Xcore = dummy;
  printf ("core fraction of He (%.3lf):  ", Ycore);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    Ycore = dummy;
  printf ("core fraction of metals (%.3lf):  ", Zcore);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    Zcore = dummy;
  printf ("Ratio of specific heats (%.3lf):  ", gamma);
  scanf ("%lf", &dummy);
  if (dummy != 0.0)
    gamma = dummy;
  printf ("Number of steps (%d):  ", steps);
  scanf ("%d", &dummy2);
  if (dummy2 != 0)
    steps = dummy2;
  printf ("Should I use the CNO cycle as well as the PP chain? (0=N, 1=Y):  ");
  scanf ("%d", &CNO);
  printf ("Are all of these values correct? (1=Y, 0=N):  ");
  scanf ("%d", &valuesOK);

  while (valuesOK != 1)
  {
    printf ("\nChoose by number the incorrect value:\n");
    printf (" 0. No correction.\n");
    printf (" 1. Mass\n");
    printf (" 2. Luminsoity\n");
    printf (" 3. Pressure\n");
    printf (" 4. Temperature\n");
    printf (" 5. Density\n");
    printf (" 6. Radius\n");
    printf (" 7. surface fraction of H\n");
    printf (" 8. surface fraction of He\n");
    printf (" 9. surface fraction of metals\n");
    printf ("10. core fraction of H\n");
    printf ("11. core fraction of He\n");
    printf ("12. core fraction of metals\n");
    printf ("13. Ratio of specific heats\n");
    printf ("14. Number of steps\n");
    printf ("15. Use CNO cycle\n");
    printf ("Which value is incorrect? : ");
    scanf ("%d", &choice);

    switch (choice)
    {
      case 0:  break;
      case 1:  printf ("Mass in g (1.99e+033, 5.97e+033):  ");
	       scanf ("%lf", &mass);
	       if (mass > 5.97e33)
	       {
		 printf ("The approximating calculations will not hold true for this mass.\n");
		 printf ("Do you still want to continue the program? (0=N, 1=Y):  ");
		 scanf ("%d", &choice);
		 if (choice == 0)     /* Mass entered too big for */
		 return (3);           /* accurate calculations */
	       }
	       break;
      case 2:  printf ("Luminosity in erg/s (3.90e+033):  ");
	       scanf ("%lf", &lum);
	       break;
      case 3:  printf ("Pressure in dyne/cm^2 (2.47e+005):  ");
	       scanf ("%lf", &pres);
	       break;
      case 4:  printf ("Temperature in K (5.78e+003):  ");
	       scanf ("%lf", &temp);
	       break;
      case 5:  printf ("Density in g/cm^3 (3.00e-007):  ");
	       scanf ("%lf", &rho);
	       break;
      case 6:  printf ("Radius in cm (6.96e+010):  ");
	       scanf ("%lf", &radius);
	       break;
      case 7:  printf ("surface fraction of H (0.710):  ");
	       scanf ("%lf", &X);
	       break;
      case 8:  printf ("surface fraction of He (0.270):  ");
	       scanf ("%lf", &Y);
	       break;
      case 9:  printf ("surface fraction of metals (0.020):  ");
	       scanf ("%lf", &Z);
	       break;
     case 10:  printf ("core fraction of H (0.330):  ");
	       scanf ("%lf", &Xcore);
	       break;
     case 11:  printf ("core fraction of He (0.630):  ");
	       scanf ("%lf", &Ycore);
	       break;
     case 12:  printf ("core fraction of metals (0.040):  ");
	       scanf ("%lf", &Zcore);
	       break;
     case 13:  printf ("Ratio of specific heats (1.666):  ");
	       scanf ("%lf", &gamma);
	       break;
     case 14:  printf ("Number of steps (100):  ");
	       scanf ("%d", &steps);
	       break;
     case 15:  printf ("Should I use the CNO cycle as well as the PP chain? (1=Y, 0=N):  ");
	       scanf ("%d", &CNO);
	       break;
     default:  printf ("That is an invalid choice.  Try again.\n");
	       break;
    }
    printf ("Have all of the wrong values been corrected? (1=Y, 0=N):  ");
    scanf ("%d", &valuesOK);
  }



  /* Set change in radius and kappa and epsilon for surface */

   delta = (radius / steps);
   kappa  = Kappa (X, Z, rho, temp);
   epsilon = Epsilon (rho, X, Z, temp, CNO);



  /* Send headings and initial values to file */

  printf ("\nStep  Mass      Lum.      Pres.     Temp.     Density   Kappa     Epsilon\n");
  printf ("       (g)     (erg/s)  (dyne/cm^2)  (K)      (g/cm^3)  (1/g)     (erg/g)\n\n");
  fprintf (output,"Radius\t\t Mass\t\t Luminosity\t Pressure\t Temperature\t Density\t Kappa\t\t Epsilon\n");
  fprintf (output," (cm) \t\t (g) \t\t  (erg/s)\t  (dyne/cm^2)\t   (K)    \t (g/cm^3)\t(1/g)\t\t (erg/g)\n\n");
  printf ("%d    %9.2le %9.2le %9.2le %9.2le %9.2le %9.2le %9.2le\n",
	   0, mass, lum, pres, temp, rho, kappa, epsilon);
  fprintf (output, "%le\t %le\t %le\t %le\t %le\t %le\t %le\t %le\n",
	   radius, mass, lum, pres, temp, rho, kappa, epsilon);



  /* Loop for calculations */

  for (count = 1; count <= steps; count++)
  {
    /* Step radius with counter */

    newradius = (radius - delta);



    /* Make the calculations for new values */

    newmass = Mass (mass, rho, newradius, delta);
    newlum = Luminosity (lum, rho, epsilon, newradius, delta);
    newpres = Pressure (pres, mass, rho, newradius, delta);
    newtemp = Temperature (temp, kappa, rho, lum, newradius, delta,
			   gamma, mass, pres);
    newrho = Density (pres, temp, X, Y, Z);
    newkappa  = Kappa (X, Z, rho, temp);
    newepsilon = Epsilon  (rho, X, Z, temp, CNO);



    /* Send results to screen and to file */

    printf ("%d    %9.2le %9.2le %9.2le %9.2le %9.2le %9.2le %9.2le\n",
	     count, newmass, newlum, newpres, newtemp, newrho, newkappa,
	     newepsilon);
    fprintf (output, "%le\t %le\t %le\t %le\t %le\t %le\t %le\t %le\n",
	      newradius, newmass, newlum, newpres, newtemp, newrho, newkappa,
	      newepsilon);



    /* Set initial values for next step to new values of last step */

    mass = newmass;
    lum = newlum;
    pres = newpres;
    temp = newtemp;
    rho = newrho;
    kappa = newkappa;
    epsilon =  newepsilon;
    radius = newradius;



    /* Check to see if the core region of the star has been reached */
    /* (if epsilon > 0.0) and if so, set X, Y, and Z to core values */

    if (epsilon >= 0.0)
    {
      X = Xcore;
      Y = Ycore;
      Z = Zcore;
    }



    /* Check to see if the mass or luminosity has become less than zero */
    /* before the radius is zero.  If so, end the program and close the */
    /* output file.  Otherwise, tell user all data is in the file and   */
    /* close the file. */

     if (((mass < 0.0) || (lum < 0.0)) && (radius != 0.0))
    {
      printf ("\n<Press any key to continue.>");
      proceed = getch();                            /* Pause before ending */
      printf ("\n");

      if ((mass < 0.0) && (lum < 0.0))
	printf ("\nThe mass and luminosity have reached zero before I could\n");
      else if (lum < 0.0)
	     printf ("\nThe luminosity has reached zero before I could\n");
	   else
	     printf ("\nThe mass has reached zero before I could\n");
      printf ("reach the center of the star.  The program is being ended\n");
      printf ("and the data file (with incomplete data) closed.\n");
      fclose (output);
      return (-1);
    }

    else if (radius == 0.0)
	 {
	   printf ("\nThe stellar center has been reached!\n");
	   printf ("A complete set of data is now held in the file.\n");
	   printf ("The file is being closed and the program completed.\n");
	   fclose (output);
	   return (0);
	  }



    /* Pause after 20 new values are calculated.  This lets the user */
    /* take a look at the data before it all rushes past on the screen. */
    /* The first time through only calculates 19 new values in order to */
    /* keep the headings and units on the screen at the same time. */

    if (count <= 20)
    {
      if ((fmod( (double) count, 19.0)) == 0.0)
      {
	printf ("\n<Press any key to continue calculations.>");
	proceed = getch();
	printf ("\n\n");
      }
    }
    else if ((fmod( (double) count, 20.0)) == 0.0)
	   {
	     printf ("\n<Press any key to continue calculations.>");
	     proceed = getch();
	     printf ("\n\n");
	   }



  }                                                  /* End of for loop */



}                                                      /* End of main */