/*

	Astronomical calculation routines

*/

#include "screensv.h"

/*  JTIME  --  Convert a Unix date and time (tm) structure to astronomical
			   Julian time (i.e. Julian date plus day fraction,
			   expressed as a double).  */

double jtime(struct tm FAR *t)
{
	double j;
	
    j = ucttoj(t->tm_year + 1900, t->tm_mon, t->tm_mday, t->tm_hour, t->tm_min, t->tm_sec);
	return j;
}

/*  UCTTOJ  --	Convert GMT date and time to astronomical
				Julian time (i.e. Julian date plus day fraction,
				expressed as a double).  */

double ucttoj(long year, int mon, int mday,
			  int hour, int min, int sec)
{

	// Algorithm as given in Meeus, Astronomical Algorithms, Chapter 7, page 61

	int a, b, m;
	long y;
	double j;

    assert(mon  >= 0 && mon  < 12);
    assert(mday >  0 && mday < 32);
    assert(hour >= 0 && hour < 24);
    assert(min  >= 0 && min  < 60);
    assert(sec  >= 0 && sec  < 60);

    m = mon + 1;
    y = year;

	if (m <= 2) {
		y--;
		m += 12;
	}

	/* Determine whether date is in Julian or Gregorian calendar based on
	   canonical date of calendar reform. */

	if ((year < 1582) || ((year == 1582) && ((mon < 9) || (mon == 9 && mday < 5)))) {
		b = 0;
	} else {
		a = ((int) (y / 100));
		b = 2 - a + (a / 4);
	}

	j = (((long) (365.25 * (y + 4716))) + ((int) (30.6001 * (m + 1))) +
				mday + b - 1524.5) +
			((sec + 60L * (min + 60L * hour)) / 86400.0);
	return j;
}

/*  JYEAR  --  Convert	Julian	date  to  year,  month, day, which are
			   returned via integer pointers to integers (note that year
			   is a long).  */

void jyear(double td, long *yy, int *mm, int *dd)
{
	double z, f, a, alpha, b, c, d, e;

	td += 0.5;
	z = floor(td);
	f = td - z;

	if (z < 2299161.0) {
		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);

	*dd = (int) (b - d - floor(30.6001 * e) + f);
	*mm = (int) ((e < 14) ? (e - 1) : (e - 13));
	*yy = (long) ((*mm > 2) ? (c - 4716) : (c - 4715));
}

/*  JHMS  --  Convert Julian time to hour, minutes, and seconds.  */

void jhms(double j, int *h, int *m, int *s)
{
    long ij;

    j += 0.5;			      /* Astronomical to civil */
    ij = (long) ((j - floor(j)) * 86400.0);
    *h = (int) (ij / 3600L);
    *m = (int) ((ij / 60L) % 60L);
    *s = (int) (ij % 60L);
}
