#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <float.h>

#define PI 3.14159265359                               /* Definition of pi */

struct HoldDate {double JD; double month; double day; double year;
		 double hour; double min; double sec; double last;
		 double lha;};
		    /* Define the structure that will hold dates and times */

/*
 *
 * T I M E   O F   M I N I M A
 * James Marshall
 *
 * This program will calculate the times when a specific phase occurs for
 * variable stars given the values for a time of minium, the period of
 * variation, the first date from which to begin listing times, the final
 * date to list times, and the phase you want the times for.
 * Among the calculations involved here are JD to Calendar date, Calendar
 * date to JD, hms UT to dec UT, dec UT to hms UT, decimal day to whole day,
 * time UT, Local time to UT, and UT to Local time.  A structure is used to
 * hold the dates and times of phase, done one at a time, and written to a
 * file, one line at a time.  This allows an unlimited number of times of
 * minima to be found and written; using an array restricts the amount of
 * times that can be stored.
 *
 */



void Instructions (void)
{
  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 calculate the times of a phase\n");
  printf ("for any variable star.  In order to make the calculations, I\n");
  printf ("need you to input a previously calculated time of minimum,\n");
  printf ("the period of variation, the date you wish to start finding\n");
  printf ("times from, the date to stop finding times at, and the phase\n");
  printf ("to find times at.  I will also need the right ascencsion (RA)\n");
  printf ("of the star being observed, the longitude you are observing at,\n");
  printf ("and the conversion factor to change from UT into Local time.\n");
  printf ("I will prompt you for the information that I need and then\n");
  printf ("print out a list of the times of phase.\n\n");
  printf ("Please enter the filename to output data to, including the path.\n");
  printf ("For example, c:\\temp\\test.dat -- but please do not use\n");
  printf ("more than 51 characters.\n\n");
}



double sign (double d_number)
{
  double d_signof;                                      /* Set up variable */

  if (d_number >= 0.0)                   /* Calculate the sign of d_number */
    d_signof = 1.0;                                   /* (+ or - ; 0 is +) */
    else d_signof = -1.0;

  return (d_signof);                       /* Send the result back to main */
}



void TimeTooHigh (double *p_hour, double *p_month, double *p_day,
		  double *p_year)
{
  /* Start long loop to do checks if time must change */

  (*p_hour) -= 24.0;            /* Convert to proper hour of one day after */
  (*p_day) += 1.0;                             /* Convert to one day after */


  if ((((*p_month) == 1.0) || ((*p_month) == 3.0) || ((*p_month) == 5.0) ||
       ((*p_month) == 7.0) || ((*p_month) == 8.0) || ((*p_month) == 10.0) ||
       ((*p_month) == 12.0)) && ((*p_day) > 31.0))
    {
      (*p_day) = 1.0;               /* See if 31 day month reached 32 days */
      (*p_month) += 1.0;                             /* and correct for it */
    }

  if ((((*p_month) == 4.0) || ((*p_month) == 6.0) || ((*p_month) == 9.0) ||
       ((*p_month) == 11.0)) && ((*p_day) > 30.0))
    {
      (*p_day) = 1.0;               /* See if 30 day month reached 31 days */
      (*p_month) += 1.0;                             /* and correct for it */
    }

  if (((*p_month) == 2.0) && ((*p_day) >= 28.0))
    {
      (*p_day) = 1.0;               /* See if 28 day month reached 29 days */
      (*p_month) += 1.0;                             /* and correct for it */
    }

  if ((*p_month) > 12.0)
    {
      (*p_month) = 1.0;                   /* See if month is past December */
      (*p_year) += 1.0;                              /* and correct for it */
    }

  if ((floor((*p_year)/4.0) == ((*p_year)/4.0)) && ((*p_year) != 2000.0))
    {                                              /* If it is a leap year */
      if (((*p_month) == 3.0) && ((*p_day) == 1.0))
	{                                /* and date was changed up to 3-1 */
	  *p_month = 2.0;
	  *p_day = 29.0;                        /* Correct it back to 2-29 */
	}
    }
}



void TimeTooLow (double *p_hour, double *p_month, double *p_day,
		 double *p_year)
{
  /* Start long loop to do checks if time must change */

  (*p_hour) += 24.0;           /* Convert to proper hour of one day before */
  (*p_day) -= 1.0;                            /* Convert to one day before */


  if ((((*p_month) == 2.0) || ((*p_month) == 4.0) || ((*p_month) == 6.0) ||
       ((*p_month) == 8.0) || ((*p_month) == 9.0) || ((*p_month) == 11.0) ||
       ((*p_month) == 1.0)) && ((*p_day) < 1.0))
    {
      (*p_day) = 31.0;              /* See if month went into 31 day month */
      (*p_month) -= 1.0;                             /* and correct for it */
    }

  if ((((*p_month) == 5.0) || ((*p_month) == 7.0) || ((*p_month) == 10.0) ||
       ((*p_month) == 12.0)) && ((*p_day) < 1.0))
    {
      (*p_day) = 30.0;              /* See if month went into 30 day month */
      (*p_month) -= 1.0;                             /* and correct for it */
    }

  if (((*p_month) == 3.0) && ((*p_day) < 1.0))
    {
      (*p_day) = 28.0;              /* See if month went into 28 day month */
      (*p_month) -= 1.0;                             /* and correct for it */
    }

  if ((*p_month) < 1.0)
    {
      (*p_month) = 12.0;                 /* See if month is before January */
      (*p_year) -= 1.0;                              /* and correct for it */
    }

  if ((floor((*p_year)/4.0) == ((*p_year)/4.0)) && ((*p_year) != 2000.0))
    {                                              /* If it is a leap year */
      if (((*p_month) == 2.0) && ((*p_day) == 28.0)) /* & date was changed */
	*p_day = 29.0;              /* down to 2-28, correct it up to 2-29 */
    }

}



void Local_to_UT (double *p_hour, double *p_month, double *p_day,
		  double *p_year, double d_correct)

{
  (*p_hour) -= d_correct;
  if (*p_hour >= 24.0)                        /* If gone more than one day */
    TimeTooHigh (p_hour, p_month, p_day, p_year);            /* Correct it */
  else
    if ((*p_hour) <= 0.0)                     /* If gone less than one day */
      TimeTooLow (p_hour, p_month, p_day, p_year);           /* Correct it */
}



void UT_to_Local (double *p_hour, double *p_month, double *p_day,
		  double *p_year, double d_correct)

{
  (*p_hour) += d_correct;
  if (*p_hour <= 0.0)                         /* If gone less than one day */
    TimeTooLow (p_hour, p_month, p_day, p_year);             /* Correct it */
  else
    if (*p_hour >= 24.0)                      /* If gone more than one day */
      TimeTooHigh (p_hour, p_month, p_day, p_year);          /* Correct it */
}



double FindDayOfYear (double d_month, double d_day, double d_year)
{
  double d_N, d_term1, d_term2;          /* Set up variable */

  d_term1 = floor((275.0 * d_month) / 9);
  if (((d_year/4)==(floor(d_year/4))) && ((d_year/400)==(floor(d_year/400))))
    d_term2 = floor((d_month + 9) / 12);
  d_term2 = 2 * floor((d_month + 9) / 12);
  d_N = d_term1 - d_term2 + d_day - 30.0;

  return (d_N);
}



void GetName (char *sname)
{
  printf ("Please enter the star's name (max. 20 chars.):  ");   /* Prompt */
  scanf ("%s20", sname);                                       /* Get name */
}



void PrintHeadings (double d_phase, char *sname)
{
  int i_proceed;

  printf ("\n<Press any key to begin calculations.>");
  i_proceed = getch();                               /* Pause before start */
  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");    /* CLS */
  printf ("\nSTAR NAME = %s\t\t\t\tPHASE = %lf\n", sname, d_phase);
  printf ("\nJulian Date\t   Date\t\tLocal time\tSidereal\tHour Angle");
  printf ("\n          \tmm-dd-yyyy\t hh:mm:ss\t hours  \t  hours   \n");
						      /* Print col. labels */
}



void WriteToFile (struct HoldDate date_time, FILE *p_file)
{

  fprintf (p_file, "%lf\t%2.0lf-%2.0lf-%4.0lf\t%2.0lf:%2.0lf:%2.0lf\t%lf\t%lf\n",
	   date_time.JD, date_time.month, date_time.day, date_time.year,
	   date_time.hour, date_time.min, date_time.sec, date_time.last,
	   date_time.lha);
}



void StructDate (FILE *p_file, double d_JD, double d_month, double d_day,
		 double d_year, double d_hour, double d_min, double d_sec,
		 double d_last, double d_lha)
{
  struct HoldDate date_time;

  date_time.JD = d_JD;
  date_time.month = floor(d_month);
  date_time.day = floor(d_day);
  date_time.year = floor(d_year);
  date_time.hour = floor(d_hour);
  date_time.min = floor(d_min);
  date_time.sec = floor(d_sec);
  date_time.last = d_last;
  date_time.lha = d_lha;

  printf ("\n%lf\t%2.0lf-%2.0lf-%4.0lf\t%2.0lf:%2.0lf:%2.0lf\t%lf\t%lf",
	  date_time.JD, date_time.month, date_time.day, date_time.year,
	  date_time.hour, date_time.min, date_time.sec, date_time.last,
	  date_time.lha);                               /* Check structure */

  WriteToFile (date_time, p_file);

}



double hmsUT_to_decUT (double d_hour, double d_min, double d_sec)
{
  double d_UT;                                         /* Set up variables */

  d_UT = (d_hour + (d_min / 60.0) + (d_sec / 3600.0));   /* Convert to dec */

  return (d_UT);                                   /* Send UT back to main */
}



double CalculateLAST (double d_N, double d_UT, double d_lambda)
{
  double d_gmst, d_omega, d_E, d_gast, d_last;         /* Set up variables */

  d_gmst = 6.6265313 + 0.0657098242 * d_N + 1.00273791 * d_UT;
						    /* Formula to get GMST */
  d_omega = 318.5111 - 0.0529538 * (d_N + d_UT / 24);
						   /* Formula to get omega */
  d_E = -0.00029 * sin (d_omega * PI / 180);
			   /* Formula to get E (omega is converted to rads */
  d_gast = d_gmst + d_E;                            /* Formula to get GAST */
  d_last = d_gast + d_lambda / 15;                  /* Fromula to get LAST */
  if (d_last >= 24.0)
    d_last -= 24.0;
  if (d_last < 0.0)                     /* Keep LAST in range 0 - 24 hours */
    d_last += 24.0;

  return (d_last);                                       /* Send back LAST */
}



double CalculateLHA (double d_last, double d_RA)
{
  double d_lha;                                         /* Set up variable */

  d_lha = d_last - d_RA;                    /* Formula to get LHA in hours */
  if (d_lha >= 24.0)
    d_lha -= 24.0;
  if (d_lha < 0.0)                       /* Keep LHA in range 0 - 24 hours */
    d_lha += 24.0;
  if (d_lha > 12.0)                    /* Put LHA in range -12 -> 12 hours */
    d_lha -= 24.0;

  return (d_lha);                                         /* Send back LHA */
}



double Date_to_JD (double d_month, double d_day, double d_year, double d_UT)
{
  /* Set up the variables */

  double d_JD, d_t1, d_t2, d_t3, d_s1, d_q1, d_q2, d_q3, d_q4;



  /* Calculate the JD  - broken down in small portions of the equation */

  d_t1 = floor ((d_month + 9.0) / 12.0);                /* Truncate part 1 */
  d_t2 = floor ((7.0 * (d_year + d_t1)) / 4.0);         /* Truncate part 2 */
  d_t3 = floor ((275.0 * d_month) / 9.0);               /* Truncate part 3 */
  d_s1 = sign (100.0 * d_year + d_month - 190002.5);   /* Get sign of part */
  d_q1 = ((367.0 * d_year) - d_t2);    /* Find first 2 terms of equation   */
  d_q2 = (d_t3 + d_day);               /* Find second 2 terms of equation  */
  d_q3 = (1721013.5 + (d_UT / 24));    /* Find third 2 terms of equation   */
  d_q4 = ((-0.5 * d_s1) + 0.5);        /* Find fourth 2 terms of equation  */
  d_JD = (d_q1 + d_q2 + d_q3 + d_q4);  /* Sum the 4 quarters to get the JD */

  return (d_JD);                               /* Send the JD back to main */
}



void JD_to_Date (FILE *p_file, double d_JD, double d_phase, double d_correct,
		 double d_RA, double d_lambda)
{
  /* Set up the variables */

  double d_Z, d_F, d_A, d_alpha, d_B, d_C, d_D, d_E, d_day, d_month, d_year,
	 d_hour, d_min, d_sec, d_JD2, *p_hour, *p_month, *p_day, *p_year,
	 d_last, d_lha, d_N, d_UT;



  /* Set up the pointers */

  p_hour = &d_hour;
  p_month = &d_month;
  p_day = &d_day;
  p_year = &d_year;



  /* Calculate the calendar date - done in small parts of the formula */

  d_JD2 = d_JD + 0.5;
  d_Z = floor (d_JD2);                         /* Integer part of JD + 0.5 */
  d_F = (d_JD2 - d_Z);                      /* Fractional part of JD + 0.5 */
  if (d_Z < 2299161.0)
    d_A = d_Z;                                        /* Term 1 in formula */
    else
    {
      d_alpha = floor ((d_Z - 1867216.25) / 36524.25);   /* Part of Opt. 2 */
      d_A = (d_Z + 1 + d_alpha - floor (d_alpha / 4));  /* Term 1 - Opt. 2 */
    }
  d_B = (d_A + 1524);                                 /* Term 2 in formula */
  d_C = floor ((d_B - 122.1) / 365.25);               /* Term 3 in formula */
  d_D = floor (365.25 * d_C);                         /* Term 4 in formula */
  d_E = floor ((d_B - d_D) / 30.6001);                /* Term 5 in formula */
  if (d_E < 13.5)
    d_month = (d_E - 1);                                 /* Find the month */
    else
      d_month = (d_E - 13);                     /* Find the month - Opt. 2 */
  if (d_month > 2.5)
    d_year = (d_C - 4716);                                /* Find the year */
    else
      d_year = (d_C - 4715);                     /* Find the year - Opt. 2 */
  d_day = (d_B - d_D - floor (30.6001 * d_E) + d_F);      /* Find dec. day */



  /* Convert decimal day to whole day and hms UT */

  d_hour = (24 * (d_day - floor (d_day)));
  d_min = (60 * (d_hour - floor (d_hour)));
  d_sec = (60 * (d_min - floor (d_min)));
  d_day = floor (d_day);



  /* Section to find day of year, local apparent sidereal time, and */
  /* local hour angle. */

  d_N = FindDayOfYear (d_month, d_day, d_year);
  d_UT = d_hour + d_min / 60 + d_sec / 3600;
  d_last = CalculateLAST (d_N, d_UT, d_lambda);
  d_lha = CalculateLHA (d_last, d_RA);



  /* Convert hms UT to Local time then send everything to write to file */

  UT_to_Local (p_hour, p_month, p_day, p_year, d_correct);
  StructDate (p_file, d_JD, d_month, d_day, d_year, d_hour, d_min, d_sec,
	      d_last, d_lha);
}



void FindMins (FILE *p_file, double d_tmin, double d_period, double d_startJD,
	      double d_stopJD, double d_phase, double d_correct, double d_RA,
	      double d_lambda)
{
  int i_count, i_proceed;

  i_count = 0;
  d_tmin += (d_period * d_phase);   /* Set d_tmin to count at proper phase */
  if (d_startJD >= d_tmin)             /* If start date is later than tmin */
    {
      do                                                     /* Run a loop */
	d_tmin += d_period;             /* Incrementing tmin by the period */
      while (d_tmin < d_startJD);     /* Until it goes over the start date */
    }
  else                               /* Start date is earlier than tmin so */
    {
      do                                                     /* Run a loop */
	d_tmin -= d_period;             /* Decrementing tmin by the period */
      while (d_tmin > d_startJD);    /* Until it goes below the start date */
    }

  while ((d_tmin >= d_startJD) && (d_tmin <= d_stopJD))
  {                          /* While tmin is between start and stop dates */
    i_count++;
    JD_to_Date (p_file, d_tmin, d_phase, d_correct, d_RA, d_lambda);
					    /* Convert JD to calendar date */
    d_tmin += d_period;                    /* Increment tmin by the period */
    if (i_count == 20)                       /* If screen near full, pause */
    {
      printf ("\n\n<Press any key to continue calculations.>");
      i_proceed = getch();
      printf ("\n\n");
      i_count = 0;
    }
  }
}



double GetTmin (void)
{
  double d_tmin;                                        /* Set up variable */

  printf ("Please enter the time of minimum (in JD):  ");        /* Prompt */
  scanf ("%lf", &d_tmin);                                 /* Get the input */

  return (d_tmin);                                 /* Send it back to main */
}



double GetPeriod (void)
{
  double d_period;                                      /* Set up variable */

  printf ("Please enter the period (in days): ");                /* Prompt */
  scanf ("%lf", &d_period);                               /* Get the input */

  return (d_period);                               /* Send it back to main */
}



double StartDate (double d_correct)
{
  double d_month, d_day, d_year, d_hour, d_min, d_sec, d_UT, d_startJD;
  double *p_hour, *p_month, *p_day, *p_year;           /* Set up variables */

  p_hour = &d_hour;
  p_month = &d_month;                                   /* Set up pointers */
  p_day = &d_day;
  p_year = &d_year;

  printf ("\nPlease enter the month, day, and year of the start date.");
  printf ("\nUse the form MM-DD-YYYY: ");                        /* Prompt */
  scanf ("%lf-%lf-%lf", &d_month, &d_day, &d_year);            /* Get date */
  printf ("Now enter the time of day (in Local time).");
  printf ("\nUse the form HH-MM-SS: ");                          /* Prompt */
  scanf ("%lf-%lf-%lf", &d_hour, &d_min, &d_sec);                /* Get UT */
  Local_to_UT (p_hour, p_month, p_day, p_year, d_correct);

  d_UT = hmsUT_to_decUT (d_hour, d_min, d_sec);      /* Convert hms to dec */
  d_startJD = Date_to_JD (d_month, d_day, d_year, d_UT);  /* Convert to JD */

  return (d_startJD);                /* Send start date in JD back to main */
}



double StopDate (double d_correct)
{
  double d_month, d_day, d_year, d_hour, d_min, d_sec, d_UT, d_stopJD;
  double *p_hour, *p_month, *p_day, *p_year;           /* Set up variables */

  p_hour = &d_hour;
  p_month = &d_month;                                   /* Set up pointers */
  p_day = &d_day;
  p_year = &d_year;

  printf ("\nPlease enter the month, day, and year of the stop date.");
  printf ("\nUse the form MM-DD-YYYY: ");                        /* Prompt */
  scanf ("%lf-%lf-%lf", &d_month, &d_day, &d_year);            /* Get date */
  printf ("Now enter the time of day (in Local time).");
  printf ("\nUse the form HH-MM-SS: ");                          /* Prompt */
  scanf ("%lf-%lf-%lf", &d_hour, &d_min, &d_sec);                /* Get UT */
  Local_to_UT (p_hour, p_month, p_day, p_year, d_correct);

  d_UT = hmsUT_to_decUT (d_hour, d_min, d_sec);      /* Convert hms to dec */
  d_stopJD = Date_to_JD (d_month, d_day, d_year, d_UT);   /* Convert to JD */

  return (d_stopJD);                 /* Send start date in JD back to main */
}



double GetPhase (void)
{
  double d_phase;                                       /* Set up variable */

  printf ("\nEnter the phase you wish to find times for:  ");    /* Prompt */
  scanf ("%lf", &d_phase);                                    /* Get phase */

  return (d_phase);                                /* Send it back to main */
}



double GetTimeZone (void)
{
  double d_correct;                                     /* Set up variable */

  printf ("\nEnter the correction from UT to Local time.");
  printf ("\nTo convert to EST, enter -5.0; MST, -7.0.");       /* Prompts */
  printf ("\nThe correction, please:  ");
  scanf ("%lf", &d_correct);                                      /* Input */

  return (d_correct);                              /* Send it back to main */
}



double GetLongitude (void)
{
  double d_lambda;                                      /* Set up variable */

  printf ("Please enter your longitude (in degrees).");
  printf ("\nUse positive for East, negative for West.");       /* Prompts */
  printf ("\nLongitude, please:  ");
  scanf ("%lf", &d_lambda);                                   /* Get input */

  return (d_lambda);                                   /* Send back lambda */
}



double GetRA (void)
{
  double d_RA, d_hour, d_min, d_sec;                   /* Set up variables */

  printf ("Please enter the Right Ascension of the star being observed.");
  printf ("\nUse the form HH-MM-SS:  ");                        /* Prompts */
  scanf ("%lf-%lf-%lf", &d_hour, &d_min, &d_sec);             /* Get input */
  d_RA = d_hour + d_min / 60 + d_sec / 3600;
					 /* Convert hms into decimal hours */

  return (d_RA);                                           /* Send back RA */
}



void PrintEndAndInfo (double d_correct, double d_RA, double d_lambda,
	double d_tmin, double d_period, double d_start, double d_stop,
	double d_phase, FILE *p_file, char *sname)
{
  fprintf (p_file, "%lf\t%2.0lf-%2.0lf-%4.0lf\t%2.0lf:%2.0lf:%2.0lf\t%lf\t%lf\n",
	   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
  fprintf (p_file, "\nJulian Date\t   Date\t\tLocal time\tSidereal\tHour Angle");
  fprintf (p_file, "\n           \tdd-mm-yyyy\t hh:mm:ss\t hours\t          hours\n");
  fprintf (p_file, "\nConversion factor -- UT to Local = %lf hours", d_correct);
  fprintf (p_file, "\nRA of observed star = %lf hours", d_RA);
  fprintf (p_file, "\nLongitude of observer = %lf degrees", d_lambda);
  fprintf (p_file, "\nTime of minimum = %lf JD", d_tmin);
  fprintf (p_file, "\nPeriod of variation = %lf days", d_period);
  fprintf (p_file, "\nPhase to observe = %lf", d_phase);
  fprintf (p_file, "\nDates from %lf JD to %lf JD", d_start, d_stop);
  fprintf (p_file, "\nStar Name:  %s", sname);
  fprintf (p_file, "\n\n");
}



int main (void)
{
  double d_tmin, d_period, d_start, d_stop, d_phase, d_correct, d_RA, d_lambda;
  char filename[51], starname[20], *sname;
  int i_count1, i_count2;                                   /* Set up vars */
  struct HoldDate date_time;
  FILE *p_file;

  sname = &starname[0];                              /* Initialize pointer */
  Instructions ();                                   /* Print instructions */
  printf ("Path and filename, please:  ");          /* Prompt for filename */
  scanf ("%s51", filename);                              /* Input filename */
  p_file = fopen (filename, "a");                             /* Open file */
  if (p_file == NULL)
  {
    printf ("Cannot open %s.  Terminating program.", filename);
    return (1);                                   /* Notify if cannot open */
  }

  d_correct = GetTimeZone ();                /* Get UT to Local correction */
  d_RA = GetRA ();                                             /* Get R.A. */
  d_lambda = GetLongitude ();                  /* Get observer's longitude */
  d_tmin = GetTmin ();                                       /* Input tmin */
  d_period = GetPeriod ();                                 /* Input period */
  d_start = StartDate (d_correct);                     /* Input start date */
  d_stop = StopDate (d_correct);                        /* Input stop date */
  d_phase = GetPhase ();                                    /* Input phase */
  GetName (sname);
  PrintHeadings (d_phase, sname);                        /* Print headings */
  FindMins (p_file, d_tmin, d_period, d_start, d_stop, d_phase, d_correct,
	    d_RA, d_lambda);                                /* Find minima */



  PrintEndAndInfo (d_correct, d_RA, d_lambda, d_tmin, d_period, d_start,
		   d_stop, d_phase, p_file, starname);
  /* Print line of zeros for file end and add extra info at end of file */

  fclose (p_file);                                       /* Close the file */

  return (0);                                               /* End program */
}

