/********************************************************/
/*                        DaTime                        */
/*  Sundry utility functions concerning date and time.  */
/********************************************************/
#include <math.h>

#include "os.h"
#include "coords.h"
#include "menu.h"
#include "sv_header.h"
#include "datime.h"

/*    ---------     Global variables.     ----------    */

/* Data array for Julian Date algorithm, based on that  */
/* given in NAO Technical Note No. 46 (January 1978):   */
  static int zfun[2][13] = {
    {0, -1, 30, 58, 89, 119, 150, 180, 211, 242, 272, 303, 333},
    {0, -2, 29, 58, 89, 119, 150, 180, 211, 242, 272, 303, 333}
  };

/*     ---------        Functions.       ----------     */

/*------------------------------------------------------*/
/* Function to convert time of day in hours to integer  */
/*                  hours and minutes.                  */
/*------------------------------------------------------*/
 void datime_hourmin(REAL time, int *hourptr, int *minptr)
{
  REAL realmin;

/* Ensure that result is a valid time.                  */
  if (time <  (REAL)0.0)
  {
    *hourptr = 0;
    *minptr  = 0;
    return;
  }

  if (time >= (REAL)24.0)
  {
    *hourptr = 23;
    *minptr  = 59;
    return;
  }

  *hourptr   = (int)time;
  realmin = (REAL)60.0 * (time - (REAL)*hourptr) + (REAL)0.5;
  *minptr    = (int)realmin;

  return;
}

/*------------------------------------------------------*/
/*   Function to calculate the days and fraction of a   */
/*  day since noon GMT on 00/01/1900.  The date is that */
/* specified in the observerstr pointed to by ptr, and  */
/*  the time of day is given by the parameters.  Time   */
/*        offset (time zone) is taken from *ptr.        */
/*------------------------------------------------------*/
  double datime_julian(observerstr *ptr, int hour, int min)
{
/* Based on algorithm in Nautical Almanac Office Tech-  */
/* nical Note No. 46 (January 1978).                    */
/* J = int(365.25*(year-1900)) + zfun[leap][month] +    */
/*     0.5 + day + ut/24                                */
/* where ut is universal time (=GMT) in hours,          */
/*       leap is 1 for a leap year and 0 otherwise,     */
/*   and zfun was defined in the section on global vars.*/

  int y = (1461 * (ptr->year - 1900)) / 4;
  double time = (double)hour + (double)min/60.0;

  return (double)y +
         (double)zfun[leap(ptr->year)][ptr->month] +
         0.5 +
         (double)ptr->day +
         (time - (double)ptr->offset)/24.0;
}

/*------------------------------------------------------*/
/*  Function to calculate local siderial time, for the  */
/* date, time and location specified. Time offset (time */
/*  zone) and longitude are taken from the observerstr  */
/* pointed to by ptr.  Other info is given by params.   */
/*   Result is in radians, in range 0 <= time < 2*PI.   */
/*------------------------------------------------------*/
  REAL datime_siderial(observerstr *ptr, double jul, int hour, int min)
{
/* Based on algorithm in Nautical Almanac Office Tech-  */
/* nical Note 46 (January 1978).                        */

/* observerstr fields offset and longit assumed valid.  */

  double t   = jul/36525.0;
  double tim = (double)hour + (double)min/60.0
              -(double)ptr->offset;
  double ghaa;
  REAL   lst;

/* Get Greenwich Hour Angle of Aries, 0 <= ghaa < 2*PI. */
  ghaa = fmod(99.6910 + 36000.7689*t + 15.0*tim + 360.0, 360.0);

/* Get local siderial time, 0 <= lst < 2*PI.            */
  lst = (REAL)ghaa*CONV - ptr->longit;
  if (lst <  ZERO  ) lst += TWO_PI;
  if (lst >= TWO_PI) lst -= TWO_PI;
  if (lst <  ZERO  ) lst  = ZERO;

  return lst;
}

/*------------------------------------------------------*/
/*  Function to calculate (one) civil time during the   */
/*     day at which the local siderial time equals      */
/*        'sid_time' (which must be in radians).        */
/*------------------------------------------------------*/
  void datime_time_today(REAL sid_time, observerstr *ptr,
                         int *hourptr, int *minptr)
{
/* Hour, min and sid fields of observerstr are assumed  */
/* to be consistent with each other.                    */

  REAL siderial_diff;
  REAL solar_diff;
  REAL tim;

/* Get sid_time into range 0 <= sid_time < 2*PI.        */
  while (sid_time <  ZERO  ) sid_time += TWO_PI;
  while (sid_time >= TWO_PI) sid_time -= TWO_PI;
  if    (sid_time <  ZERO  ) sid_time  = ZERO;

/* Find time difference (in siderial hours) between     */
/* siderial time now and required siderial time.        */
/* Diff is in range -24 < diff < 24                     */
  siderial_diff = (sid_time - ptr->sid)/CONV/(REAL)15.0;

/* Convert difference to solar hours.                   */
  solar_diff = SID_TO_SOL * siderial_diff;

/* Calculate required civil time in decimal hours.      */
  tim = (REAL)ptr->hour + (REAL)ptr->min/(REAL)60.0 + \
       solar_diff;

/* Bring time within range by adding or subtracting one */
/* siderial day, if necessary.                          */
  while (tim <  (REAL) 0.0) tim += (REAL)24.0 * SID_TO_SOL;
  while (tim >= (REAL)24.0) tim -= (REAL)24.0 * SID_TO_SOL;

/* Convert tim to integer hours and minutes.            */
  datime_hourmin(tim, hourptr, minptr);

  return;
}
