#include <math.h>
#include "project.h"

static int is_leap (unsigned year)
{
  if (year % 4)
    return 0;

  if (year % 100)
    return 1;

  if (year % 400)
    return 0;

  return 1;
}

EPOCH time_epoch_sub (EPOCH a, EPOCH b)
{
  EPOCH ret;

  ret.ns = a.ns - b.ns;
  ret.s = a.s - b.s;

  if (!ret.s) return ret;

  if (ret.ns < 0) {
    ret.s--;
    ret.ns += 1000000000;
  }

  return ret;
}


EPOCH time_epoch_add (EPOCH a, EPOCH b)
{
  EPOCH ret;

  ret.ns = a.ns + b.ns;
  ret.s = a.s + b.s;

  while (ret.ns > 1000000000) {
    ret.s++;
    ret.ns -= 1000000000;
  }

  return ret;
}


UTC time_epoch_to_utc (EPOCH epoch)
{
  UTC u = {0};
  uint64_t day;
  unsigned y400, y100, y4;

  day = epoch.s / 86400;
  epoch.s -= day * 86400;

  day += 134774;

  u.wday = day % 7;
  u.wday++;


  y400 = day / 146097;        /*146097 days in 400 years */
  day -= (y400 * 146097);

  y100 = day / 36524;         /*36524 days in 100 years */
  day -= (y100 * 36524);

  y4 = day / 1461;            /*1461 days in 4 years */
  day -= (y4 * 1461);

  /* This may look redundant but 31 Dec in year 4 is special case */
  if (day < 1095) {           /*1095 days in 3 years */
    u.year = day / 365;     /*365 days in a year */
    day -= (365 * u.year);
  } else {
    u.year = 3;
    day -= 1095;
  }


  /* Now put it all back together */
  u.year += 1601;
  u.year += 4 * y4;
  u.year += 100 * y100;
  u.year += 400 * y400;


  u.jday = day + 1;

  u.is_leap = is_leap (u.year);


  if (!u.is_leap) {
    /*Days and months for ordinary years */
    if (u.jday < 32) {
      u.month = 1;
      u.mday = u.jday;
    } else if (u.jday < 60) {
      u.month = 2;
      u.mday = u.jday - 31;
    } else if (u.jday < 91) {
      u.month = 3;
      u.mday = u.jday - 59;
    } else if (u.jday < 121) {
      u.month = 4;
      u.mday = u.jday - 90;
    } else if (u.jday < 152) {
      u.month = 5;
      u.mday = u.jday - 120;
    } else if (u.jday < 182) {
      u.month = 6;
      u.mday = u.jday - 151;
    } else if (u.jday < 213) {
      u.month = 7;
      u.mday = u.jday - 181;
    } else if (u.jday < 244) {
      u.month = 8;
      u.mday = u.jday - 212;
    } else if (u.jday < 274) {
      u.month = 9;
      u.mday = u.jday - 243;
    } else if (u.jday < 305) {
      u.month = 10;
      u.mday = u.jday - 273;
    } else if (u.jday < 335) {
      u.month = 11;
      u.mday = u.jday - 304;
    } else {
      u.month = 12;
      u.mday = u.jday - 334;
    }
  } else {
    /*And leap years */
    if (u.jday < 32) {
      u.month = 1;
      u.mday = u.jday;
    } else if (u.jday < 61) {
      u.month = 2;
      u.mday = u.jday - 31;
    } else if (u.jday < 92) {
      u.month = 3;
      u.mday = u.jday - 60;
    } else if (u.jday < 122) {
      u.month = 4;
      u.mday = u.jday - 91;
    } else if (u.jday < 153) {
      u.month = 5;
      u.mday = u.jday - 121;
    } else if (u.jday < 183) {
      u.month = 6;
      u.mday = u.jday - 152;
    } else if (u.jday < 214) {
      u.month = 7;
      u.mday = u.jday - 182;
    } else if (u.jday < 245) {
      u.month = 8;
      u.mday = u.jday - 213;
    } else if (u.jday < 275) {
      u.month = 9;
      u.mday = u.jday - 244;
    } else if (u.jday < 306) {
      u.month = 10;
      u.mday = u.jday - 274;
    } else if (u.jday < 336) {
      u.month = 11;
      u.mday = u.jday - 305;
    } else {
      u.month = 12;
      u.mday = u.jday - 335;
    }
  }

  u.hour = epoch.s / 3600;
  epoch.s -= u.hour * 3600;
  u.minute = epoch.s / 60;
  epoch.s -= u.minute * 60;
  u.second = epoch.s;

  u.nanosecond = epoch.ns;

  return u;
}

EPOCH  time_utc_to_epoch (UTC u)
{
  unsigned y400;
  unsigned y100;
  unsigned y4;

  EPOCH ret;

  static int const mdays[] =
  { 0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 };
  static int const lmdays[] =
  { 0, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 };

  u.is_leap = is_leap (u.year);

  if (u.year < 100) u.year += 2000;


  if (!u.jday) {
    if (u.is_leap)
      u.jday = u.mday + lmdays[u.month];
    else
      u.jday = u.mday + mdays[u.month];
  }

  u.year -= 1601;
  y400 = u.year / 400;
  u.year -= y400 * 400;
  y100 = u.year / 100;
  u.year -= y100 * 100;
  y4 = u.year / 4;
  u.year -= y4 * 4;



  ret.s = u.jday - 1;
  ret.s += u.year * 365;
  ret.s += y4 * 1461;
  ret.s += y100 * 36524;
  ret.s += y400 * 146097;

  ret.s -= 134774;

  ret.s *= 86400;

  ret.s += u.second;
  ret.s += u.minute * 60;
  ret.s += u.hour * 3600;

  ret.ns = u.nanosecond;

  return ret;
}

void utc_to_str (char *dst, UTC *u)
{
  const char *dname[] = {"Sun", "Mon", "Tue", "Wed", "Thr", "Fri", "Sat", "Sun"};
  const char *mname[] = {"", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
  sprintf (dst, "%s %04d-%s-%02d %02d:%02d:%02d.%09d", dname[u->wday], u->year, mname[u->month], u->mday, u->hour, u->minute, u->second, u->nanosecond);
}


void time_print_utc (const char *p, UTC *u, const char *t)
{
  const char *dname[] = {"Sun", "Mon", "Tue", "Wed", "Thr", "Fri", "Sat", "Sun"};
  const char *mname[] = {"", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
  printf ("%s%s %04d-%s-%02d %02d:%02d:%02d.%09d %s\r\n", p ? p : "", dname[u->wday], u->year, mname[u->month], u->mday, u->hour, u->minute, u->second, u->nanosecond, t ? t : "");
}

void time_print_epoch (const char *p, EPOCH e, const char *t)
{
  UTC u = time_epoch_to_utc (e);
  time_print_utc (p, &u, t);
}


double
time_utc_to_tjd (UTC u)
{
  unsigned y400;
  unsigned y100;
  unsigned y4;

  double ret;
  unsigned jd;

  static int const mdays[] =
  { 0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 };
  static int const lmdays[] =
  { 0, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 };

  u.is_leap = is_leap (u.year);

  if (u.year < 100)
    u.year += 2000;


  if (!u.jday) {
    if (u.is_leap)
      u.jday = u.mday + lmdays[u.month];
    else
      u.jday = u.mday + mdays[u.month];
  }

  u.year -= 1601;
  y400 = u.year / 400;
  u.year -= y400 * 400;
  y100 = u.year / 100;
  u.year -= y100 * 100;
  y4 = u.year / 4;
  u.year -= y4 * 4;

  jd = u.jday - 1;
  jd += u.year * 365;
  jd += y4 * 1461;
  jd += y100 * 36524;
  jd += y400 * 146097;


  jd += 2305813;
  jd -= 2451545;



  ret = (double) u.nanosecond;
  ret /= 1000000000.;
  ret += (double) u.second;
  ret /= 60.;
  ret += (double) u.minute;
  ret /= 60.;
  ret += (double) u.hour;
  ret /= 24.;

  ret += .5;

  ret += (double) jd;

  return ret;
}


EPOCH time_tjd_to_epoch (double tjd)
{
  EPOCH e;


  tjd -= 2440587.5;

  tjd *= 86400.;
  e.s = (uint64_t) floor (tjd);

  tjd -= (double) e.s;

  tjd *= 1000000000.;
  e.ns = (unsigned)  floor (tjd);

  return e;
}

UTC time_tjd_to_utc (double tjd)
{
  return time_epoch_to_utc (time_tjd_to_epoch (tjd));
}



RA ra_normalize (RA r)
{

  while (r.ra < 0.) {
    r.ra += 360.;
    r.days--;
  }

  if (r.ra >= 3600.0) {
    r.days = (unsigned) (floor (r.ra / 360.) + .5);
    r.ra = remainder (r.ra, 360.);
  } else while (r.ra >= 360.0) {
      r.ra -= 360.;
      r.days++;
    }

  return r;
}


RA time_utc_to_ra (UTC u)
{
  RA r;
  double tjd = time_utc_to_tjd (u);
  double T = tjd / 36525.;

  r.ra =
    280.46061837 + (360.98564736629 * tjd) + (0.000387933 * T * T) -
    (T * T * T / 38710000.0);

  r.days = 0;

  return ra_normalize (r);
}



RA time_st_to_ra (ST st)
{
  RA r;
  unsigned i;

  r.days = st.days;

  i = st.hour;
  i *= 60;
  i += st.minute;
  i *= 60;
  i += st.second;

  r.ra = (double) st.nanosecond;
  r.ra /= 1000000000.;
  r.ra += (double) i;

  r.ra /= 240.;

  return r;
}



ST
time_ra_to_st (RA r)
{
  ST ret;
  unsigned i;

  r = ra_normalize (r);

  ret.days = r.days;

  r.ra *= 240.0;

  i = (int) (floor (r.ra) + .5);
  r.ra -= (double) i;
  r.ra *= 1000000000.;

  ret.nanosecond = (unsigned) (r.ra + .5);
  ret.second = i % 60;
  i /= 60;
  ret.minute = i % 60;
  i /= 60;
  ret.hour = i;


  return ret;
}



ST time_utc_to_lst (UTC u, double lon)
{
  RA r;

  r = time_utc_to_ra (u);

  r.ra += lon;

  return time_ra_to_st (r);
}



#if 0
int time_ra_cmp (double a, double b)
{
  return (ra_normalize (a - b) < 180.) ? -1 : 1;;
}


/*Find the next occuring time for RA tra */

static int test_side (EPOCH m, double tra)
{

  double ra;
  UTC u;

  u = time_epoch_to_utc (m);
  ra = time_utc_to_ra (u);

  return time_ra_cmp (ra, tra);
}


EPOCH time_ra_to_next_epoch (EPOCH l, double tra)
{
  EPOCH r, m;
  unsigned n;

  int dog = 48;

  printf ("time_ra_to_next_epoch %.9f\n", tra);

  /* XXX: we should use newton raphson, but */
  /* 1) we're on team hooke */
  /* 2) we want to limit the domain of solutions */
  /* So IB works better */


  while ((test_side (l, tra) != 1)  && dog--)
    l.s += 3600;

  r = l;

  while ((test_side (r, tra) != -1) && dog--)
    r.s += 3600;

  for (n = 0; n < 64; ++n) {

    m = time_epoch_sub (r, l);

    if (m.s < 0) {
      m.s = 0;
      m.ns = 0;
      break;
    }

    if (m.s & 1)
      m.ns += 1000000000;

    m.s /= 2;
    m.ns /= 2;

    m = time_epoch_add (l, m);

    if (test_side (m, tra) < 0)
      r = m;
    else
      l = m;
  }

  return m;
}

#endif





/*
 * Returns the EPOCH for the start of the metric year which
 * starts in the custumary year y
 *
 * the metric year begins on the autumnal equinox (at the paris observatory natch.)
 */
static EPOCH start_of_year (unsigned y)
{
  double  tjd = autumnal_equinox (y);
  UTC u;

  tjd += (2.3372305555 / 360); /* Offset of paris meridian from greenwich */

  u = time_tjd_to_utc (tjd);

  u.hour =
    u.minute =
      u.second =
        u.nanosecond = 0;

  return time_utc_to_epoch (u);
}

MTIME time_to_metric (EPOCH e, UTC u)
{
  MTIME ret;

  static EPOCH year_start, year_end;
  static unsigned year;

  unsigned s;
  unsigned d;

  double fd;

  EPOCH equinox;

  if ((e.s > year_end.s) || (e.s < year_start.s)) {
    /*calculate the equinox in this customary year */

    equinox = start_of_year (u.year);

    if (e.s < equinox.s) {
      year_start = start_of_year (u.year - 1);
      year_end = equinox;
      year = u.year - 1792;
    } else {
      year_start = equinox;
      year_end = start_of_year (u.year + 1);
      year = u.year - 1791;
    }
  }

  ret.year = year;

  e = time_epoch_sub (e, year_start);
  s = e.s;
  d = s / 86400;
  s -= d * 86400;

  ret.jday = d;

  ret.month = (d / 30) + 1;
  ret.mday = (d % 30) + 1;
  ret.wday = d % 10;

  /* Now for the horror of time of day */

  fd = (double) e.ns;
  fd /= 1000000000.;
  fd += (double) s;
  fd /= 86400.;
  fd *= 100000;

  s = (int) fd;
  fd -= (double) s;
  fd *= 1000000000.;

  ret.nanosecond = (int) fd;

  ret.second = s % 100;
  s -= ret.second;
  s /= 100;

  ret.minute = s % 100;
  s -= ret.minute;
  s /= 100;

  ret.hour = s;

  return ret;
}


MTIME time_utc_to_metric (UTC u)
{
  return  time_to_metric (time_utc_to_epoch (u), u);
}


MTIME time_epoch_to_metric (EPOCH e)
{
  return  time_to_metric (e, time_epoch_to_utc (e));
}



void time_print_metric (const char *p, MTIME *m, const char *t)
{
  const char *dname[] = {"Pri", "Duo", "Tri", "Qua", "Qui", "Sex", "Sep", "Oct", "Non", "Dec" };
  const char *mname[] = {"", "Ven", "Bru", "Fri", "Niv", "Plu", "Ven", "Ger", "Flo", "Pra", "Mes", "The", "Fru", "San"};
  printf ("%s%s %03d-%s-%02d %01d:%02d:%02d.%09d %s\r\n", p ? p : "", dname[m->wday], m->year, mname[m->month], m->mday, m->hour, m->minute, m->second, m->nanosecond, t ? t : "");
}