#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include "sgdp4h.h"

#define LIM 384
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
#define XKMPAU 149597879.691	// AU in km

struct map
{
  int site_id;
  double mjd;
  float saltmin, alt;
  char nfd[LIM], observer[32], datadir[LIM];
  double lat, lng;
  float length;
} m;
void allnight (void);
void sunpos_xyz (double, xyz_t *, double *, double *);
double modulo (double, double);
void equatorial2horizontal (double, double, double, double *, double *);
void mjd2date (double mjd, char *date);
double gmst (double mjd);
void nfd_now (char *s);
double nfd2mjd (char *date);
void get_site (int site_id);
double date2mjd (int year, int month, double day);
void usage (void);

int
main (int argc, char *argv[])
{
  int arg = 0;
  char *env;

  // Default settings
  strcpy (m.observer, "Unknown");
  m.site_id = 0;
  m.mjd = -1.0;
  m.saltmin = -6.0;
  m.alt = 0.0;
  m.lat = 0.0;
  m.lng = 0.0;

  // Get default site
  env = getenv ("ST_DATADIR");
  if (env != NULL)
    {
      strcpy (m.datadir, env);
    }
  else
    {
      printf ("ST_DATADIR environment variable not found.\n");
    }
  env = getenv ("ST_COSPAR");
  if (env != NULL)
    {
      get_site (atoi (env));
    }
  else
    {
      printf ("ST_COSPAR environment variable not found.\n");
    }

  // Get current time
  nfd_now (m.nfd);
  m.mjd = nfd2mjd (m.nfd);

  // Decode options
  while ((arg = getopt (argc, argv, "t:s:S:h")) != -1)
    {
      switch (arg)
	{

	case 't':
	  strcpy (m.nfd, optarg);
	  m.mjd = nfd2mjd (m.nfd);
	  break;

	case 'S':
	  m.saltmin = atof (optarg);
	  break;

	case 's':
	  get_site (atoi (optarg));
	  break;

	case 'h':
	  usage ();
	  return 0;
	  break;

	default:
	  usage ();
	  return 0;
	}
    }

  // Compute set/rise times of sun
  allnight ();

  return 0;
}

// Usage
void
usage (void)
{
  printf ("allnight t:s:S:\n\n");
  printf ("t    date/time (yyyy-mm-ddThh:mm:ss.sss) [default: now]\n");
  printf ("S    Minimum sun altitude\n");
  printf ("s    site (COSPAR)\n");

  return;
}

// Convert equatorial into horizontal coordinates
void
equatorial2horizontal (double mjd, double ra, double de, double *azi,
		       double *alt)
{
  double h;

  h = gmst (mjd) + m.lng - ra;

  *azi =
    modulo (atan2
	    (sin (h * D2R),
	     cos (h * D2R) * sin (m.lat * D2R) -
	     tan (de * D2R) * cos (m.lat * D2R)) * R2D, 360.0);
  *alt =
    asin (sin (m.lat * D2R) * sin (de * D2R) +
	  cos (m.lat * D2R) * cos (de * D2R) * cos (h * D2R)) * R2D;

  return;
}

// Return x modulo y [0,y)
double
modulo (double x, double y)
{
  x = fmod (x, y);
  if (x < 0.0)
    x += y;

  return x;
}

// Solar position
void
sunpos_xyz (double mjd, xyz_t * pos, double *ra, double *de)
{
  double jd, t, l0, m, e, c, r;
  double n, s, ecl;

  jd = mjd + 2400000.5;
  t = (jd - 2451545.0) / 36525.0;
  l0 = modulo (280.46646 + t * (36000.76983 + t * 0.0003032), 360.0) * D2R;
  m = modulo (357.52911 + t * (35999.05029 - t * 0.0001537), 360.0) * D2R;
  e = 0.016708634 + t * (-0.000042037 - t * 0.0000001267);
  c = (1.914602 + t * (-0.004817 - t * 0.000014)) * sin (m) * D2R;
  c += (0.019993 - 0.000101 * t) * sin (2.0 * m) * D2R;
  c += 0.000289 * sin (3.0 * m) * D2R;

  r = 1.000001018 * (1.0 - e * e) / (1.0 + e * cos (m + c));
  n = modulo (125.04 - 1934.136 * t, 360.0) * D2R;
  s = l0 + c + (-0.00569 - 0.00478 * sin (n)) * D2R;
  ecl =
    (23.43929111 +
     (-46.8150 * t - 0.00059 * t * t + 0.001813 * t * t * t) / 3600.0 +
     0.00256 * cos (n)) * D2R;

  *ra = atan2 (cos (ecl) * sin (s), cos (s)) * R2D;
  *de = asin (sin (ecl) * sin (s)) * R2D;

  pos->x = r * cos (*de * D2R) * cos (*ra * D2R) * XKMPAU;
  pos->y = r * cos (*de * D2R) * sin (*ra * D2R) * XKMPAU;
  pos->z = r * sin (*de * D2R) * XKMPAU;

  return;
}

void
allnight (void)
{
  int flag;
  xyz_t sunpos;
  double ra, de, azi, alt, alt0;
  double mjd, mjdrise = -1.0, mjdset = -1.0;
  char nfd[32];

  // Find solar altitude at reference time
  sunpos_xyz (m.mjd, &sunpos, &ra, &de);
  equatorial2horizontal (m.mjd, ra, de, &azi, &alt);

  // Sun below limit, find rise, then set
  if (alt < m.saltmin)
    {
      for (flag = 0, mjd = m.mjd; mjd < m.mjd + 0.5; mjd += 1.0 / 86400)
	{
	  sunpos_xyz (mjd, &sunpos, &ra, &de);
	  equatorial2horizontal (mjd, ra, de, &azi, &alt);

	  if (flag != 0)
	    {
	      if (alt > m.saltmin && alt0 <= m.saltmin)
		mjdrise = mjd;
	    }

	  if (flag == 0)
	    flag = 1;

	  alt0 = alt;
	}
      for (flag = 0, mjd = m.mjd - 0.5; mjd < m.mjd; mjd += 1.0 / 86400)
	{
	  sunpos_xyz (mjd, &sunpos, &ra, &de);
	  equatorial2horizontal (mjd, ra, de, &azi, &alt);

	  if (flag != 0)
	    {
	      if (alt < m.saltmin && alt0 >= m.saltmin)
		mjdset = mjd;
	    }

	  if (flag == 0)
	    flag = 1;

	  alt0 = alt;
	}
      // Sun above limit, find set, and rise
    }
  else
    {
      for (flag = 0, mjd = m.mjd; mjd < m.mjd + 1.0; mjd += 1.0 / 86400)
	{
	  sunpos_xyz (mjd, &sunpos, &ra, &de);
	  equatorial2horizontal (mjd, ra, de, &azi, &alt);

	  if (flag != 0)
	    {
	      if (alt > m.saltmin && alt0 <= m.saltmin)
		mjdrise = mjd;
	      if (alt < m.saltmin && alt0 >= m.saltmin)
		mjdset = mjd;
	    }

	  if (flag == 0)
	    flag = 1;

	  alt0 = alt;
	}
    }

  m.mjd = mjdset;
  mjd2date (m.mjd, m.nfd);
  mjd2date (mjdrise, nfd);
  printf ("%s %s\n", m.nfd, nfd);

  return;
}

// Compute Date from Julian Day
void
mjd2date (double mjd, char *date)
{
  double f, jd, dday;
  int z, alpha, a, b, c, d, e;
  int year, month, day, hour, min;
  float sec, x;

  jd = mjd + 2400000.5;
  jd += 0.5;

  z = floor (jd);
  f = fmod (jd, 1.);

  if (z < 2299161)
    a = z;
  else
    {
      alpha = floor ((z - 1867216.25) / 36524.25);
      a = z + 1 + alpha - floor (alpha / 4.);
    }
  b = a + 1524;
  c = floor ((b - 122.1) / 365.25);
  d = floor (365.25 * c);
  e = floor ((b - d) / 30.6001);

  dday = b - d - floor (30.6001 * e) + f;
  if (e < 14)
    month = e - 1;
  else
    month = e - 13;

  if (month > 2)
    year = c - 4716;
  else
    year = c - 4715;

  day = (int) floor (dday);
  x = 24.0 * (dday - day);
  x = 3600. * fabs (x) + 0.0001;
  sec = fmod (x, 60.);
  x = (x - sec) / 60.;
  min = fmod (x, 60.);
  x = (x - min) / 60.;
  hour = x;
  sec = floor (1000.0 * sec) / 1000.0;

  sprintf (date, "%04d-%02d-%02dT%02d:%02d:%02.0f", year, month, day, hour,
	   min, sec);

  return;
}

// Greenwich Mean Sidereal Time
double
gmst (double mjd)
{
  double t, gmst;

  t = (mjd - 51544.5) / 36525.0;

  gmst =
    modulo (280.46061837 + 360.98564736629 * (mjd - 51544.5) +
	    t * t * (0.000387933 - t / 38710000), 360.0);

  return gmst;
}

// Present nfd
void
nfd_now (char *s)
{
  time_t rawtime;
  struct tm *ptm;

  // Get UTC time
  time (&rawtime);
  ptm = gmtime (&rawtime);

  sprintf (s, "%04d-%02d-%02dT%02d:%02d:%02d", ptm->tm_year + 1900,
	   ptm->tm_mon + 1, ptm->tm_mday, ptm->tm_hour, ptm->tm_min,
	   ptm->tm_sec);

  return;
}

// nfd2mjd
double
nfd2mjd (char *date)
{
  int year, month, day, hour, min, sec;
  double mjd, dday;

  sscanf (date, "%04d-%02d-%02dT%02d:%02d:%02d", &year, &month, &day, &hour,
	  &min, &sec);
  dday = day + hour / 24.0 + min / 1440.0 + sec / 86400.0;

  mjd = date2mjd (year, month, dday);

  return mjd;
}

// Compute Julian Day from Date
double
date2mjd (int year, int month, double day)
{
  int a, b;
  double jd;

  if (month < 3)
    {
      year--;
      month += 12;
    }

  a = floor (year / 100.);
  b = 2. - a + floor (a / 4.);

  if (year < 1582)
    b = 0;
  if (year == 1582 && month < 10)
    b = 0;
  if (year == 1582 && month == 10 && day <= 4)
    b = 0;

  jd =
    floor (365.25 * (year + 4716)) + floor (30.6001 * (month + 1)) + day + b -
    1524.5;

  return jd - 2400000.5;
}

// Get observing site
void
get_site (int site_id)
{
  int i = 0;
  char line[LIM];
  FILE *file;
  int id;
  double lat, lng;
  float alt;
  char abbrev[3], observer[64], filename[LIM];

  sprintf (filename, "%s/data/sites.txt", m.datadir);
  file = fopen (filename, "r");
  if (file == NULL)
    {
      printf ("File with site information not found!\n");
      return;
    }
  while (fgets (line, LIM, file) != NULL)
    {
      // Skip
      if (strstr (line, "#") != NULL)
	continue;

      // Strip newline
      line[strlen (line) - 1] = '\0';

      // Read data
      sscanf (line, "%4d %2s %lf %lf %f", &id, abbrev, &lat, &lng, &alt);
      strcpy (observer, line + 38);

      // Change to km
      alt /= 1000.0;

      if (id == site_id)
	{
	  m.lat = lat;
	  m.lng = lng;
	  m.alt = alt;
	  m.site_id = id;
	  strcpy (m.observer, observer);
	}

    }
  fclose (file);

  return;
}