/*****************************************************************************
 * solar.c
 *
 * DESCRIPTION
 *    This file contains the information needed to compute sunrise/sunset
 * and isNightPeriod()
 *
 * sunlib.C: (comment from George Trojan) The formulas are from:
 *             Explanatory Supplement to the Astronomical Almanac,
 *             Kenneth Seidelmann (ed), University Science Books 1992,
 *             pp 484-487 (eqns 9.311 through 9.33.3)
 *             Assume const Sun semidiameter = 16' and refraction = 34'
 *             Astronomical Almanac, 2001
 *
 * HISTORY
 *   07/16/1998 George Trojan, GSC/TDL: Created sunlib.C
 *   03/04/2004 Joe Lang (MDL): "converted the minimal amount of C++ syntax
 *              into standard C syntax" creating srss.c
 *   5/2006 Arthur Taylor (MDL): Modified from srss.c
 *
 * NOTES
 *****************************************************************************
 */
#include <math.h>
#include <time.h>
#include "solar.h"
#include "clock.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/*
* Unix time for J2000.0 12h UT
*/
const time_t Epoch = 946728000;
const double LONG_OF_SUN = 280.460;
const double LONG_ABERRATION = 0.9856474;
const double ANOMALY = 357.528;
const double ANOMALY_ABERRATION = 0.9856003;
const double ECLIPTIC_ABERRATION = 0.0000004;
const float SEC_PER_DAY = 86400.0;
const float DAYS_PER_CENTURY = 36524.0;

enum { RISET, SUN_UP, SUN_DOWN };

double radian (double arg)
{
   return (arg * M_PI / 180.0);
}
double degree (double arg)
{
   return (arg * 180.0 / M_PI);
}

/*-------------------------------------------------------------------------*/
/*
* returns UNIX time in UTC
*/
time_t AbsTime (int year, int month, int day, int hour, int min)
{
   time_t now = time (0);
   struct tm *tms = gmtime (&now);

   tms->tm_sec = 0;
   tms->tm_min = min;
   tms->tm_hour = hour;
   tms->tm_mday = day;
   tms->tm_mon = month - 1;
   tms->tm_year = year - 1900;
   /* changed because timezone is not in ANSI C */
/*  return (mktime (tms) - timezone);  */
   return (mktime (tms) - Clock_GetTimeZone ());
}

/*-------------------------------------------------------------------------*/
/*
*  returns number of centuries from J2000.0 12h UT
*/
static double NumCenturies (time_t absTime)
{
   return (difftime (absTime, Epoch) / (SEC_PER_DAY * DAYS_PER_CENTURY));
}

/*-------------------------------------------------------------------------*/
/*
*  Mean longitude of Sun, corrected for aberration, [deg]
*  num_cent: number of centuries from J2000.0 12h UT
*/
static double MeanLonOfSun (double num_cent)
{
   double adjustment = LONG_ABERRATION * DAYS_PER_CENTURY;
   double ml = LONG_OF_SUN + adjustment * num_cent;

   while (ml >= 360)
      ml -= 360;
   while (ml < 0)
      ml += 360;
   return (ml);
}

/*-------------------------------------------------------------------------*/
/*
*  Mean anomaly [deg]
*  num_cent: number of centuries from J2000.0 12h UT
*/
static double MeanAnomaly (double num_cent)
{
   double adjustment = ANOMALY_ABERRATION * DAYS_PER_CENTURY;
   double g = ANOMALY + adjustment * num_cent;

   while (g >= 360)
      g -= 360;
   while (g < 0)
      g += 360;
   return (g);
}

/*-------------------------------------------------------------------------*/
/*
*  Ecliptic longitude [deg]
*  mean_anom: mean anomaly [deg]
*  mean_lon: mean longitude of Sun [deg]
*/
static double EclipticLon (double mean_anom, double mean_lon)
{
   return (mean_lon + 1.915 * sin (radian (mean_anom)) +
           0.020 * sin (radian (2.0 * mean_anom)));
}

/*-------------------------------------------------------------------------*/
/*
*  Obliquity of ecliptic [deg]
*  num_cent: number of centuries from J2000.0 12h UT
*/
static double OblOfEcliptic (double num_cent)
{
   double adjustment = ECLIPTIC_ABERRATION * DAYS_PER_CENTURY;

   return (23.439 - adjustment * num_cent);
}

/*-------------------------------------------------------------------------*/
/*
*  Right ascension
*  obl_ecl: obliquity of ecliptic [deg]
*  ecl_lon: ecliptic longitude [deg]
*/
static double Ascension (double ecl_lon, double obl_ecl)
{
   double t;

   ecl_lon = radian (ecl_lon);
   obl_ecl = radian (obl_ecl);
   t = tan (obl_ecl / 2.) * tan (obl_ecl / 2.);
   return (degree (ecl_lon - (t * sin (2. * ecl_lon)) +
                   (0.5 * t * t * sin (4. * ecl_lon))));
}

/*-------------------------------------------------------------------------*/
/*
*  Equation of time [deg]
*  mean_lon: mean_longitude [deg]
*  alpha: ascension [deg]
*/
static double EqnOfTime (double mean_lon, double alpha)
{
   return (mean_lon - alpha);
}

/*-------------------------------------------------------------------------*/
/*
*  declination [deg]
*  obl_ecl: obliquity of ecliptic [deg]
*  ecl_lon: ecliptic longitude [deg]
*/
static double Declination (double obl_ecl, double ecl_lon)
{
   return (degree (asin (sin (radian (obl_ecl)) * sin (radian (ecl_lon)))));
}

/*-------------------------------------------------------------------------*/
/*
*  Calculates distance of sun from earth in astronomical units
*  mean_anom: mean anomaly [deg]
*/
static double DistanceFromSun (double mean_anom)
{
   return (1.00014 - (0.01671 * cos (radian (mean_anom))) -
           (0.00014 * cos (radian (2. * mean_anom))));
}

/*-------------------------------------------------------------------------*/
/*
*  Hour angle: solves the equation wr to omega [deg]
*     cos(theta) = sin(delta)*sin(fi) + cos(delta)*cos(fi)*cos(omega)
*     where:
*        theta:  zenith angle
*        delta:  declination
*        fi:     latitude
*        omega:  hour angle
*  decl: declination [deg]
*  lat: latitude [deg]
*  h: correction for refraction + semidiameter [deg]
*/
static double HourAngle (double decl, double lat, double h)
{
   double tmp;

   h = radian (h);
   decl = radian (decl);
   lat = radian (lat);
   tmp = sin (h) / (cos (decl) * cos (lat)) - tan (decl) * tan (lat);
   if (tmp >= 1.0)
      return (0.0);
   else if (tmp <= -1.0)
      return (180);
   else
      return (degree (acos (tmp)));
}

/*-------------------------------------------------------------------------*/
/*
*  Calculates semidiameter of sun based on distance of sun from earth.
*  dist: distance of sun from earth [au]
*/
static double SemidiameterOfSun (double dist)
{
   return (0.2666 / dist);
}

/*-------------------------------------------------------------------------*/
/*  Calculate the sunrise or sunset time
*  T: Number of Centuries since Jan 1, 2001, 12Z
   lat: Station latitude [deg]
   lon: Station longitude [deg]
   sunrise: 0 for sunset, 1 for sunrise
   zone: adjustment to UTC in hours to obtain local time
*/
double CalcTime (double T, double lat, double lon, int sunrise, double UT)
{
   int n;
/*
 * precision: 6 s
 */
   static const double eps = 6.0 / 3600.0;
   double UT1;

/* 10 iterations max or until when the consecutive values differ < 6 s
*/
   for (n = 0; n < 10; n++) {
      double L = MeanLonOfSun (T);
      double G = MeanAnomaly (T);
      double lambda = EclipticLon (G, L);
      double epsilon = OblOfEcliptic (T);
      double alpha = Ascension (lambda, epsilon);
      double E = EqnOfTime (L, alpha);
      double delta = Declination (epsilon, lambda);
      double R = DistanceFromSun (G);
      double SD = SemidiameterOfSun (R);
/*
*  Assume 34 arcmin for refraction
*/
      double h = -(34. / 60. + SD);
      double t = HourAngle (delta, lat, h);
      if (fabs (t) == 0 || fabs (t) == 180)
         return (t);
      if (sunrise == 0) {
         t = -t;
      }
      UT1 = 12.0 - (E + lon + t) / 15.0;
      T += (UT1 - UT) / (DAYS_PER_CENTURY * 24.0);
      if (fabs (UT1 - UT) < eps)
         return (UT1);
      UT = UT1;
   }
   return (UT1);
}

/*-------------------------------------------------------------------------*/
/*
*  calculates sunrise and sunset [sidereal days]
*  year, month, day: date
*  lat: latitude [deg]
*  lon: longitude [deg]
*  riset: returned time [hour LT (local time)]
*  sunrise: 1 for sunrise, 0 for sunset
*  zone: number of hours from UTC
*/
int RiseSet (int year, int month, int day, double lat, double lon,
             double *riset, int sunrise, int zone)
{
   double UT = 12.0 - lon / 15.0; /* start value */
   int min = 0;
   double T = NumCenturies (AbsTime (year, month, day, (int) (UT), min));
   double LT;
/*   Calculate the sunrise or sunset */
   double UT1 = CalcTime (T, lat, lon, sunrise, UT);
/* if no sunrise or sunset, exit this routine
*/

   if (fabs (UT1) == 0)
      return (SUN_DOWN);
/*  We need to consider the possibility that the sun set before it rose today.
 *  This can happen at high latitudes, where the sunset may set close to midnight,
 *  which is shifted into the next day because the selection of time zone doesn't
 *  coincide well with the position of the sun.  In other words, we need to consider
 *  the possibility that the sunset associated with yesterday, actually occurred
 *  today.
*/
   else if (fabs (UT1) == 180) {
      T = T - (1.0 / DAYS_PER_CENTURY);
      UT1 = CalcTime (T, lat, lon, sunrise, UT);
/* if the sun is up all day or down all day as a result of the adjustment to LT,
   exit this routine
*/
      if (fabs (UT1) == 0)
         return (SUN_DOWN);
      else if (fabs (UT1) == 180)
         return (SUN_UP);
/* Otherwise, compute the sunrise/sunset time in local time coordinates
*/
      LT = UT1 + zone;
/*  The sunrise or sunset doesn't occur in this calendar day as a result of
    of the adjustment to LT
*/
      if (LT < 24)
         return (SUN_UP);
      else {
         *riset = LT - 24;
         return (RISET);
      }
   }

/*    Now check to determine if we've calculated a sunrise/sunset for the
      correct day in local time.  If not adjust the day and try again.
*/
   LT = UT1 + zone;

   if (LT > 24) {
      T = T - (1.0 / DAYS_PER_CENTURY);
      UT1 = CalcTime (T, lat, lon, sunrise, UT);
/* if the sun is up all day or down all day as a result of the adjustment to LT,
   exit this routine
*/
      if (fabs (UT1) == 0)
         return (SUN_DOWN);
      else if (fabs (UT1) == 180)
         return (SUN_UP);

/*  Otherwise, compute the sunrise/sunset time in local time coordinates
*/
      LT = UT1 + zone;

/*  The sunrise or sunset doesn't occur in this calendar day as a result of
    of the adjustment to LT
*/
      if (LT < 24) {
         if (sunrise == 0)
            return (SUN_DOWN);
         else
            return (SUN_UP);
      }
   } else if (LT < 0) {
      T = T + (1.0 / DAYS_PER_CENTURY);
      UT1 = CalcTime (T, lat, lon, sunrise, UT);
      LT = UT1 + zone;
/*  The sunrise or sunset doesn't occur in this calendar day as a result
    of the adjustment to LT
*/
      if (LT > 0) {
         if (sunrise == 0)
            return (SUN_DOWN);
         else
            return (SUN_UP);
      }
   }
   LT = UT1 + zone;
   while (LT >= 24.0)
      LT -= 24.0;
   while (LT < 0.0)
      LT += 24.0;
   *riset = LT;
   return (RISET);
}

/*-------------------------------------------------------------------------*/
/* Return's the time of the rising or setting sun (f_sunrise = 1 for sunrise)
 * in UTC on the date specified by year/mon/day in UTC at location lat/lon */
/* returns 1 if sun always down.
 * returns 2 if sun always up.
 * returns 0 if time is valid. */
int sunTime (sInt4 year, int mon, int day, double lat, double lon,
             int f_sunrise, int *hh, int *mm, int *ss)
{
   double riset;
   int rc;

   /* num_off_UTC = 0 is number of hours offset from UTC. */
   rc = RiseSet (year, mon, day, lat, lon, &riset, f_sunrise, 0);

   switch (rc) {
      case SUN_DOWN:
         *hh = 0;
         *mm = 0;
         *ss = 0;
         return 1;
      case SUN_UP:
         *hh = 0;
         *mm = 0;
         *ss = 0;
         return 2;
   }
   *ss = (int) (riset * 3600);
   *hh = *ss / 3600;
   *mm = (*ss - *hh * 3600) / 60;
   *ss = *ss % 60;
   /* *ss = *ss - *hh * 3600 - *mm * 60; */
   if (*hh == 24)
      *hh = 0;
   return 0;
}

/*-------------------------------------------------------------------------*/
/* Returns true if the time (in UTC) at lat/lon is during the night.
 * returns false if the time (in UTC) at lat/lon is during the day.
 * zone (EST is 4) = amount to shift from UTC to LST
 */
#define ISDARK 1
#define ISLIGHT 0
int isNightPeriod (double time, double lat, double lon, int zone)
{
   sInt4 year;
   int mon;
   int day;
   int hour;
   int min;
   double sec;
   double rise;
   double set;
   double now;
   int rc;

   time = time - zone * 3600;
   Clock_PrintDate (time, &year, &mon, &day, &hour, &min, &sec);
   now = hour + min / 60. + sec / 3600.;

   /* sunrise is 1 so we get sunrise */
   /* num_off_UTC = 0 is number of hours offset from UTC. */
   rc = RiseSet (year, mon, day, lat, lon, &rise, 1, -1 * zone);
   switch (rc) {
      case SUN_DOWN:
         return ISDARK;
      case SUN_UP:
         return ISLIGHT;
   }

   /* sunrise is 0 so we get sunset */
   /* num_off_UTC = 0 is number of hours offset from UTC. */
   rc = RiseSet (year, mon, day, lat, lon, &set, 0, -1 * zone);
   switch (rc) {
      case SUN_DOWN:
         return ISDARK;
      case SUN_UP:
         return ISLIGHT;
   }

   if (rise < set) {
      if (now < rise)
         return ISDARK;
      else if (now < set)
         return ISLIGHT;
      else
         return ISDARK;
   } else {
      if (now < set)
         return ISLIGHT;
      else if (now < rise)
         return ISDARK;
      else
         return ISLIGHT;
   }
}

#ifdef DEBUG_SOLAR
int main (int argc, char **argv)
{
   double time;
   double lat;
   double lon;
   sInt4 year;
   int month;
   int day;
   int hour;
   int min;
   double sec;
   int i;
   char f_isNight;

   lat = 38;
   lon = -78;
   time = Clock_Seconds();

   for (i = 0; i < 100; ++i) {
      Clock_PrintDate (time, &year, &month, &day, &hour, &min, &sec);
      f_isNight = isNightPeriod (time, lat, lon, 4);
      if (f_isNight) {
         printf ("Night: %d/%d/%d %d:%d\n", month, day, year, hour, min);
      } else {
         printf ("Day: %d/%d/%d %d:%d\n", month, day, year, hour, min);
      }
      /* Add 1 hour to time */
      time = time + 3600;
   }
}
#endif

#undef ISDARK
#undef ISLIGHT