/*  Planets
 *           A program to determine the position of
 *	     the Sun., Mer., Ven., Mars, Jup, Sat & Ura
 *
 *  reference Jean Meeus's book, " Astronomical Formulae for
 *  Calculators
 *			program by
 *
 *			F.T. Mendenhall
 *			ihnp4!inuxe!fred
 *
 *		Copyright 1985 by F.T. Mendenhall
 *		all rights reserved
 *
 * (This copyright prevents you from selling the program - the
 *  author grants anyone the right to copy and install the program
 *  on any machine it will run on)
 !
 ! Modified by Alan Paeth, awpaeth@watcgl, January 1987 for
 !   (1) non-interactive mode (uses current GMT)
 !   (2) output in simplified Yale format for use with star charter software
 !   (3) lines 365-380 rewritten as a simplier C expression (awpaeth, Apr '87)
 !
 ! Modified by Petri Launiainen, pl@intrin.FI, February 1987 for
 !   easier LOGFILE definition/modification in Makefile
 !
 ! Modified by Alan Paeth, September 1987, for command line flags.
 !    The program continues to provide reduced Yale output in the old style,
 !     but new versions of "starchart" accept either. The old format allows
 !     representation of magnitudes below -.99 (useful for planets), but no
 !     spectral or stellar identification, which is not useful for planets.
 !
 ! Modified by Jim Cobb, March 1988, to compute moon positions as well.
 !    Also, to  change slightly the computation of epli, the obliquity
 !     of the ecliptic.  Formerly epli had a term 
 !		+0.00256*cos(omeg*rad).
 !     This term is a compensation for computing the apparent position
 !     of the sun (Chapter 18 in Meeus).  But for converting
 !     ecliptical coordinates to right ascension and declination we want
 !     the mean value (that is, we don't want this term).
 ! 
 ! Modified by Bob Leivian, Sept 1989, to compute local azimuth and altitude
 !     given a lat/lon, also fixed time so you can replace a single item
 !     without having to supply all (i.e. for today at 9:00 put -t 9) 
 !     also ran the code through 'cb' to standardize the indention.
 !     I split the code in to more managible pieces.
 !     Also printed out time in local and UTC as well as JD
 !
 ! See "Time Notes" section below if you plan to rewrite the user interface.
 */

#include <stdio.h>
#include <math.h>

/* default latitude longitude (of Tempe, AZ) -- replace with your favorite place */
#define DEFLONG -111.94
#define DEFLAT 33.32
#define DEFZONE 7 /* MST */

#ifndef PLANETFILE
#define PLANETFILE "./planet.star"
#endif

double	pie, rad;
double	htod(), atof(), kepler(), truean();
double	longi(), lati(), poly(), to_int(), range();
int	cls(), helio_trans(), geo_trans(), speak();

double	planet_pos();
double	sidereal_time();
int	zone;
double	longitude = DEFLONG;
double	latitude = DEFLAT;
double	GHA_Aries;


FILE    *logfile;
char	*progname;

/*
 * Times Notes.
 * 
 * This software would ideally use a generic UNIX system call to which would
 * provide day and date info as integers, plus current time west of GMT.
 * These values would become the default settings for each of the -m -d -y -t
 * & -z. Unfortunately, time on across many UNIX systems does not (to my
 * knowledge) exist as a uniform subroutine call.
 * 
 * Instead, we presume that the very low level "gettimeofday" (GMT in seconds,
 * from 1 Jan 1970) exists on all UNIX installations, and use this subroutine
 * to use the present time of day when no command line parameters are
 * present, as defaults for the current time have not been filled in.
 * 
 * When any flag-style command line parameters are present, hardcoded defaults
 * are substitude for -m, -d, -y -t switches, which may then be overridden.
 * The -z flag is filled in by a companion return value from "gettimeofday".
 * This approach makes for poor (on account of hard coding) defaults.
 * 
 * A worthwhile change would be to simplify the command line control structure
 * and default mechanism by using a higher level system call, or (safer),
 * write the routine to convert Unix GMT into day, month, year, so that it
 * lives within this program. This would allow for more reasonable defaults,
 * and still merely one UNIX system call to get GMT time. Perhaps sources
 * from a UNIX system may be studied to do this correctly.
 * 
 * Expertise in "time" on a specific UNIX system is *NOT* what is required --
 * what *IS* required is knowledge and access to enough independent UNIX
 * systems to insure (to a high probability) that the code remains portable.
 * In other words, only low-level system calls known to work across a large
 * range of systems should be employed.
 * 
 */

#include <sys/time.h>		/* for getting current GMT */

double
current_time(m, d, y, t, z)
int	*m;
int	*d;
int	*y;
double	*t;
int	*z;
{
#define SECSPERDAY (60.0 * 60.0 * 24.0)

#ifdef UNIX
#define GMT1970 2440587.5
	struct timeval tv;
	struct timezone tz;
	struct tm *p;

	gettimeofday(&tv, &tz);
	p = localtime(&tv);

	*y = p->tm_year + 1900;			/* local year */
	*m = p->tm_mon + 1;				/* local month */
	*d = p->tm_mday;				/* local day of month */

	*t = p->tm_hour + (p->tm_min / 60.0);	/* local time (0.0-23.999) */
	*z = tz.tz_minuteswest;	    /* local time hrs west of Greenwich */

	return (GMT1970 + tv.tv_sec / SECSPERDAY);
#endif 

#ifdef AMIGA
	time_t amiga_now;
	struct tm *p;

#define GMT1978 2443509.5
	amiga_now = time(&amiga_now);

    p = localtime(&amiga_now);

	*t = p->tm_hour + (p->tm_min / 60.0);	/* local time (0.0-23.999) */
	*z = DEFZONE * 60;
    *y = p->tm_year + 1900;			/* local year */
    *m = p->tm_mon + 1;				/* local month */
    *d = p->tm_mday;				/* local day of month */
    return GMT1978 + (amiga_now / SECSPERDAY) + (DEFZONE / 24.0);
#endif

#ifdef NO_TIME_AVAIL
#define MST1989 2447528.291667
	/* no time is avail, pick a default i.e. Noon Jan 1 1989 */
	*y = 1989;
	*m = 1;
	*d = 1;

	*t = 12.0;
	*z = DEFZONE * 60.;
	return MST1989;
#endif
}

double
commandtime(argc, argv)
char	**argv;
{
	int	mm, dd, yy, aa1, bb1;
	double	jd, year, month, day, time, timez;
	int	j;
	double	now;
	static char	*usage = 
"\nusage:planet {current time used}\nor\tplanet juliandate\nor\tplanet [-t hh.mm -z zone -y year -m mon -d day ]";

	/* get current time for defaults */
	now = current_time(&mm, &dd, &yy, &time, &zone);
	timez = zone / 60.; /* time zone in hours */

	/* no args - use current time */
	if (argc == 1)
		return now;

	if ((argv[1][0] != '-')) {
		/* one numeric arg - use it as the julian date */
		if (argc == 2)
			return (atof(argv[1]));
		die("no switches - %s", usage);
	} else {
		/* flags present, set up defaults */

		/* process any user overrides */
		for (j = 1; j < argc; j++) {
			if (argv[j][0] != '-')
				die("unknown argument - %s", usage /*argv[j]*/);
			switch (argv[j][1]) {
			case 't': 
				time = htod(argv[++j]); 
				break;
			case 'z': 
				timez = htod(argv[++j]);
				if (timez < 0 || timez > 24)
					die("time zone must be 0 to 24 west of UTC\n");
                                zone = timez * 60.;
				break;
			case 'y': 
				yy = atoi(argv[++j]); 
				break;
			case 'm': 
				mm = atoi(argv[++j]);
				if (mm < 1 || mm > 12)
					die("month must be 1..12\n");
				break;
			case 'd': 
				dd = atoi(argv[++j]);
				if (dd < 1 || dd > 31)
					die("day must be 1..31\n");
				break;
			case 'l':
				switch (argv[j][2]) {
				case 'o':
					longitude = atof(argv[++j]);
					if (longitude < -180 || longitude > 180)
						die("longitude must be +-180\n");
					break;
				case 'a':
					latitude = atof(argv[++j]);
					if (latitude < -90 || latitude > 90)
						die("latitude must be +-90\n");
					break;
				default:
					die("specify lat or lon\n");
				}
				break;
			default:
				die("unknown switch - %s", usage /*argv[j]*/ );
			}
			if (j == argc)
				die("trailing command line flag - %s", argv[j - 1]);
		}
	}
	day = ((float) dd + (time + timez) / 24.0);
	if (mm > 2) {
		year = yy;
		month = mm;
	} else {
		year = yy - 1;
		month = mm + 12;
	}
	aa1 = (int) (year / 100);
	bb1 = 2 - aa1 + (int) (aa1 / 4);
	jd = to_int(365.25 * year) + to_int(30.6001 * (month + 1));
	jd = jd + day + 1720994.5;
	if ((year + month / 100) > 1582.10)
		jd = jd + bb1;
	return (jd);
}


caltim(jd_frac, hour, min)
double	jd_frac;
int	*hour;
int	*min;
{
	double	x;

	x = jd_frac * 24.0;
	*hour = (int) floor(x);
	x -= floor(x);
	x *= 60.0;
	*min = (int) floor(x);
}


main(argc, argv)
int	argc;
char	**argv;
{
	double	jd, T0, epli;
	double	temp;
        long    int_jd;
	int	y, m, d, h, min, lh, lmin;
	char	*str;

	progname = argv[0];
#ifdef AMIGA
	logfile = 0; /* don't write a log */
#else
	if ((logfile = fopen(PLANETFILE, "w")) == 0)
		/*if ((logfile = fopen("/tmp/planet.star", "w")) == 0) {
        	
			die("fail on opening logfile %s\n", "/tmp/planet.star");
		}*/
#endif

	cls();

	pie = 3.1415926535898;
	rad = pie / 180.0;

	/*  calculate time T0 from 1900 January 0.5 ET */
	jd = commandtime(argc, argv);

	/* jd to calendar day */
        int_jd = jd +.5001;
	caldat(int_jd, &m, &d, &y);

	/* julian day to UTC */
	temp = jd - .4999;
	caltim(temp - floor(temp), &h, &min);
	printf("At %d/%d %d %d:%02d UTC (Julian: %.3f)\n", m, d, y, h, min, jd);

	/* UTC to Local */
	temp -= (zone / 1440.);
	caltim(temp - floor(temp), &lh, &lmin);
        int_jd = jd + .5001 - (zone / 1440.);
	caldat(int_jd, &m, &d, &y);
	if (lh > 12) {
		lh -= 12;
		str = "pm";
	} else {
		str = "am";
		if (lh == 0)
			lh = 12;
	}
	printf("local %d/%d %d:%02d%s (", m, d, lh, lmin, str);

	/* compute exact local time for a given longitude */
	temp = jd - .5 + (longitude / 360.);
	caltim(temp - floor(temp), &lh, &lmin);
	GHA_Aries =  15. * sidereal_time(jd);

	printf("%2d:%02d true local at lat:%1.1f lon:%1.1f, GHA Aries:%1.1f)\n",
	    lh, lmin, latitude, longitude, GHA_Aries);

	T0 = (jd - 2415020.0) / 36525.0;
	printf("\nOBJECT      R.A.              DEC          DISTANCE   altitude   azimuth\n");

	epli = planet_pos(jd, T0);

	moon_pos(T0, epli);

	putchar('\n');
	if (logfile)
		close(logfile);
        logfile = 0;

    /* show local position of some stars also */
	speak(101.1714, -16.7013, 0., -1.46, " " , "Sirius");
	speak(213.7956, 19.236, 0., -0.04, " " , "Arcturus");
	speak(88.6508, 7.4057, 0., 0.5, " " , "Betelgeuse");
	exit(0);
}				/* end of program main */


int
cls()
{
#define CLEARCNT 1
	int	i;
	for (i = 0; i < CLEARCNT; i++)
		putchar('\n');
}


double
poly(a0, a1, a2, a3, t)
double	a0, a1, a2, a3, t;
{
	return (a0 + a1 * t + a2 * t * t + a3 * t * t * t);
}


double
to_int(z)
double	z;
{
	long	trunk;

	trunk = (long) z;
	z = (double) trunk;
	return (z);
}


double
kepler(e, M)
double	e, M;
{
	double	corr, e0, E0, E1;
	e0 = e / rad;
	corr = 1;
	E0 = M;
	while (corr > 0.000001) {
		corr = (M + e0 * sin(E0 * rad) - E0) / (1 - e * cos(E0 * rad));
		E1 = E0 + corr;
		if (corr < 0) {
			corr = -1.0 * corr;
		}
		E0 = E1;
	}

	return (E1);
}


double
truean(e, E)
double	e, E;
{
	double	nu;

	nu = 2.0 * atan(sqrt((1 + e) / (1 - e)) * tan(E * rad / 2));
	nu = nu / rad;
	if (nu < 0.0) {
		nu = nu + 360.0;
	}
	if (fabs(nu - E) > 90.0) {
		if (nu > 180.0) {
			nu = nu - 180.0;
		} else {
			nu = nu + 180.0;
		}
	}
	return (nu);
}


double
longi(w2, i, u)
double	w2, i, u;
{
	double	x, y, l;

	y = cos(i * rad) * sin(u * rad);
	x = cos(u * rad);
	l = atan2(y, x);
	l = l / rad;
	if (l < 0.0) {
		l = l + 360.0;
	}
	return (l + w2);
}


double
lati(u, i)
double	u, i;
{
	double	b;
	b = asin(sin(u * rad) * sin(i * rad)) / rad;
	return (b);
}


int
speak(ra, dec, dis, mag, sym, str)
double	ra, dec, dis, mag;
char	*sym, *str;
{
	int	h, m, s, deg;
	double	lha, sha, altitude, azimuth;
	;

	if (ra < 0) {
		ra = ra + 360.0;
	}
	sha = 360. - ra;
	ra = ra / 15.0;		/* convert from degs to hours */
	h = (int) to_int(ra);
	m = (int) ((ra - h) * 60);
	s = (int) (((ra - h) * 60 - m) * 60);
	printf("%-12s%2dh %2dm %2ds ", str, h, m, s);
	if (logfile) 
		fprintf(logfile, "%02d%02d%02d", h, m, s);
	deg = (int) to_int(dec);
	m = (int) ((dec - deg) * 60);
	s = (int) (((dec - deg) * 60 - m) * 60);
	if (m < 0) {
		m = m * -1;
	}
	if (s < 0) {
		s = s * -1;
	}
#ifdef AMIGA
	/* its printf doesn't know about signed */
	if (deg > 0)
		printf("   +%2ddeg %2dm %2ds ", deg, m, s);
	else
		printf("   %3ddeg %2dm %2ds ", deg, m, s);
#else
	printf("   %+3ddeg %2dm %2ds ", deg, m, s);
#endif

	if (dis > 0.0001)
		printf("  %6.3f", dis);
	else
		printf("        ");
	if (logfile)
		if (mag < 0.0)
			fprintf(logfile, "%+03d%02d-%2d%s       %s\n",
			    deg, m, -(int) (10 * mag), sym, str);
		else
			fprintf(logfile, "%+03d%02d%3d%s       %s\n",
			    deg, m, (int) (100 * mag), sym, str);

	/* if we have a valid angle */
	if (GHA_Aries > 0) {
		/* now tell where it is at a given place on earth */
		lha = range(GHA_Aries + sha + longitude);

		altitude = sin(latitude * rad) * sin(dec * rad) + 
		    cos(latitude * rad) * cos(dec * rad) * cos(lha * rad);
		altitude = asin(altitude) / rad;

		azimuth = sin(lha * rad) / 
		    (cos(lha * rad) * sin(latitude * rad) - 
		    tan(dec * rad) * cos(latitude * rad));

		/* correct for quadrant */
		if (lha <= 180.) {
			if (azimuth > 0)
				azimuth = 180. + atan(azimuth) / rad;
			else
				azimuth = 360. + atan(azimuth) / rad;
		} else {
			if (azimuth <= 0.)
				azimuth = 180. + atan(azimuth) / rad;
			else
				azimuth = atan(azimuth) / rad;
		}

		printf("      %5.1f     %5.1f\n", altitude, azimuth);
	} else
		printf("\n");
	return (0);
}


double
range(val)
double	val;
{
	while (val < 0.0) {
		val = val + 360.0;
	}
	if (val > 360.0) {
		val = val - (to_int(val / 360.0) * 360.0);
	}
	return (val);
}


/*
 * Helio_trans converts heliocentric ecliptical coordinates to geocentric
 * right ascension and declination.  It also outputs the converted value.
 */
int
helio_trans(r, b, ll, Stheta, Sr, epli, mag, sym, str)
double	r, b, ll, Stheta, Sr, epli, mag;
char	*sym, *str;
{
	double	N, D, lambda, delta, beta;

	N = r * cos(b * rad) * sin((ll - Stheta) * rad);
	D = r * cos(b * rad) * cos((ll - Stheta) * rad) + Sr;
	lambda = atan2(N, D) / rad;
	if (lambda < 0.0) {
		lambda = lambda + 360.0;
	}
	lambda = range(lambda + Stheta);
	delta = sqrt(N * N + D * D + (r * sin(b * rad)) * (r * sin(b * rad)));
	beta = asin(r * sin(b * rad) / delta) / rad;
	geo_trans(lambda, beta, epli, delta, mag, sym, str);
	return (0);
}


/*
 * Geo_trans converts geocentric ecliptical coordinates to geocentric right
 * ascension and declination.  It also outputs the converted value.  Note
 * that the output coordinates are referred to the mean equinox of date.  If
 * these coordinates are to be used in conjunction with a star database, use
 * the epoch program to convert the planet's coordinates to the epoch for the
 * star database.
 */
int
geo_trans(lambda, beta, epli, delta, mag, sym, str)
double	lambda, beta, epli, delta, mag;
char	*sym, *str;
{
	double	N, D, RA, DEC;

	N = sin(lambda * rad) * cos(epli * rad)
	-tan(beta * rad) * sin(epli * rad);
	D = cos(lambda * rad);
	RA = atan2(N, D) / rad;
	DEC = asin(sin(beta * rad) * cos(epli * rad)
	     + cos(beta * rad) * sin(epli * rad) * sin(lambda * rad)) / rad;
	speak(RA, DEC, delta, mag, sym, str);
	return (0);
}


double
htod(s)
char	*s;
{
	/*
     * htod(str) reads floating point strings of the form {+-}hh.{mm} thus
     * allowing for fractional values in base 60 (ie, degrees/minutes).
     */
	double	x, sign;
	int	full, frac;
	char	*t;
	t = s - 1;
	while (*++t) {
		if (((*t < '0') || (*t > '9')) && (*t != '.') && (*t != '+') && (*t != '-'))
			die("non-digit in dd.mm style numeric argument");
	}
	if (s == NULL)
		return 0.0;
	full = frac = 0;
	sign = 1.0;
	if (*s == '-') {
		sign = -1.0;
		s++;
	} else if (*s == '+')
		s++;
	while (*s && *s != '.')
		full = 10 * full + *s++ - '0';
	if (*s++ == '.') {
		if (*s)
			frac = 10 * (*s++ - '0');
		if (*s)
			frac += *s++ - '0';
		if (frac > 59)
			frac = 59;
	}
	x = (double) full + ((double) frac) / 60.0;
	return sign * x;
}


die(a, b)
char	*a, *b;
{
	fprintf(stderr, "%s: ", progname);
	fprintf(stderr, a, b);
	fprintf(stderr, "\n");
	fflush(stderr);
	exit(1);
}