/*
 *  COLA -- a satellite close-approach finder
 *  Copyright (C) 1996 Matthew Francey
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  Send bugs reports, comments, critique, etc, to mdf@angoss.com.
 */
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <ctype.h>

#include "sgp.h"

#define	J2000	2451545.0		/* you-know-when                */
#define	U1970	 -10957.5		/* unix epoch relative to J2000 */

/*
 * convert a numeric date:
 *
 *	YYMMDD.fraction
 * or
 *	 YYDDD.fraction
 *
 * into the number of solar days from J2000.
 */
long	sgp_century = 1950L;

double
sgp_date(double t)
{
	double	T = floor(t), JD;
	long	y, m, d;

	if(T < 99365.0) {
		y = (long) floor(fmod(T*(1.0/1000.0), 100.0));
		m = 1L;
		d = (long) floor(fmod(T,             1000.0));
	} else {
		y = (long) floor(fmod(T*(1.0/10000.0), 100.0));
		m = (long) floor(fmod(T*(1.0/100.0),   100.0));
		d = (long) floor(fmod(T,               100.0));
	}

	if((y += 100L*(sgp_century/100L)) < sgp_century)
		y += 100L;
	JD = (double) (367L*y - 7L*(y + (m + 9L)/12L)/4L
		- 3L*((y + (m - 9L)/7L)/100L + 1L)/4L
		+ 275L*m/9L + d + 1721029L);
	JD += fmod(t, 1.0) - 0.5;
	return JD - J2000;
}

/*
 * return the current date/time
 */
time_t sgp_dnow = 0L;

double
sgp_now(void)
{
	return U1970 + (1.0/86400.0)*(double) (time(NULL) + sgp_dnow);
}

/*
 * convert to standard unix time
 */
time_t
sgp_unix(double t)
{
	t = floor((t - U1970)*86400.0 + 0.5);
	if(t < 0.0)
		t = 0.0;
	if(t > 2147483648.0)
		t = 2147483648.0;
	return (time_t) t;
}

/*
 * parts of a day -- used by jd_str() below to convert times...
 */
static	double	parts[] = {
	(86400.0/86400.0),
	(36000.0/86400.0),
	( 3600.0/86400.0),
	(  600.0/86400.0),
	(   60.0/86400.0),
	(   10.0/86400.0),
	(    1.0/86400.0)
};
	
/*
 * convert a user-entered string into a date/time
 */
double
sgp_getdate(char *p)
{
	double	t;
	char	*q, c, c1;
	int	i;

	/*
	 * special cases to make life easier
	 */
	if(strcmp(p, "now") == 0)
		return sgp_now();
	if(*p == 'J')
		return atof(p + 1) - J2000;

	/*
	 * date/time separation characters
	 */
	if((q = strchr(p, ':')) == (char *)NULL)
		if((q = strchr(p, '-')) == (char *)NULL)
			if((q = strchr(p, '+')) == (char *)NULL)
				q = strchr(p, '.');

	/*
	 * if none of ":-+.", assume it is a simple "date" (e.g., 900420)
	 * and return the time (tz-corrected)
	 */
	if(q == (char *)NULL || (c1 = *q) == '.')
		return sgp_date(atof(p));

	/*
	 * otherwise, decode the 'date' portion -- UTC relative
	 */
	*q++ = '\0';
	if(*p == '\0')
		t = floor(sgp_now()) - 0.5;	/* today */
	else
		t = sgp_date(atof(p));		/* whenever */

	/*
	 * '+' advances a day, '-' reverses a day
	 */
	if(c1 == '+' || c1 == '-') {
		do {
			t += (c1 == '+') ? 1.0 : -1.0;
		} while((c1 = *q++) == '+' || c1 == '-');
		--q;
	}

	/*
	 * decode the 'time' portion
	 */
	i = 1;
	while((c = *q++) != '\0' && isdigit(c) && i < 7)
		t += (double) (c - '0')*parts[i++];
	if(c == '.')
		t += atof(--q)*parts[i-1];
	return t;
}

/*
 * And perform the inverse of the above 'jd' function.
 *
 *	sgp_ymd() returns the YYMMDD.xxx format
 *
 *	sgp_ydd() returns the YYDDD.xxx format
 */
static	int	months[] = {
	31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31
};
static	int	dt_y, dt_n, dt_m, dt_d;

static void
dt(int n)
{
	dt_y = 80 + 4*(n/1461);
	if(dt_y >= 100)
		dt_y -= 100;
	if((n %= 1461) > 365) {
		dt_y += 1 + (n -= 366)/365;
		if((n %= 365) > 58)
			++n;
	}
	dt_n = n;
	for(dt_m=1; dt_m<=12 && n >= months[dt_m-1]; n -= months[dt_m++ - 1])
		;
	dt_d = n + 1;
}

double
sgp_ymd(double t)
{
	dt((int) floor(t -= -7305.5));
	return 10000.0*(double)dt_y
	       + 100.0*(double)dt_m
	       +       (double)dt_d
	       + fmod(t, 1.0);
}

double
sgp_ydd(double t)
{
	dt((int) floor(t -= -7305.5));
	return 1000.0*(double)dt_y + (double)dt_n + fmod(t, 1.0);
}

/*
 * the Greenwich Mean Sideral Time -- returned as an angle
 */
double
sgp_gmst(double t)
{
	double	x;

	x = 0.7790572733 + 1.002737909350795*t;
	return (2.0*M_PI)*(x - floor(x));
}