/*  Utility programs for unit conversions and calendar
 */

/*double PI = 3.14159265358979323846;*/
#include "kep.h"
#if __STDC__
double caltoj (long, int, double);
#else
double caltoj();
#endif

char *intfmt = "%d";
char *lngfmt = "%ld";
char *dblfmt = "%lf";
char *strfmt = "%s";

/* Display Right Ascension and Declination
 * from input equatorial rectangular unit vector.
 * Output vector pol[] contains R.A., Dec., and radius.
 */
int showrd( msg, p, pol )
char *msg;
double p[], pol[];
{
double x, y, r;
int i;

r = 0.0;
for( i=0; i<3; i++ )
	{
	x = p[i];
	r += x * x;
	}
r = sqrt(r);

x = zatan2( p[0], p[1] );
pol[0] = x;

y = asin( p[2]/r );
pol[1] = y;

pol[2] = r;

if (prtflg != 0)
  {
    printf( "%s  R.A. ", msg );
    hms( x );
    printf( "Dec. " );
    dms( y );
    printf( "\n" );
  }
return(0);
}




/* Display magnitude of correction vector
 * in arc seconds
 */
int showcor( strng, p, dp )
char *strng;
double p[], dp[];
{
double p1[3], dr, dd;
int i;

if( prtflg == 0 )
	return(0);

for( i=0; i<3; i++ )
	p1[i] = p[i] + dp[i];

deltap( p, p1, &dr, &dd );
printf( "%s dRA %.3fs dDec %.2f\"\n", strng, RTS*dr/15.0, RTS*dd );
return(0);
}



/* Radians to degrees, minutes, seconds
 */
int dms( x )
double x;
{
double s;
int d, m;

s = x * RTD;
if( s < 0.0 )
	{
	printf( " -" );
	s = -s;
	}
else
	printf( "  " );
d = (int) s;
s -= d;
s *= 60;
m = (int) s;
s -= m;
s *= 60;
printf( "%3dd %02d\' %05.2f\"  ", d, m, s );
return(0);
}




/* Radians to hours, minutes, seconds
 */
#define RTOH (12.0/PI)

int hms( x )
double x;
{
int h, m;
long sint, sfrac;
double s;

s = x * RTOH;
if( s < 0.0 )
	s += 24.0;
h = (int) s;
s -= h;
s *= 60;
m = (int) s;
s -= m;
s *= 60;
/* Handle shillings and pence roundoff. */
sfrac = (long) (1000.0 * s + 0.5);
if( sfrac >= 60000L )
  {
    sfrac -= 60000L;
    m += 1;
    if( m >= 60 )
      {
	m -= 60;
	h += 1;
      }
  }
sint = sfrac / 1000;
sfrac -= sint * 1000;
printf( "%3dh %02dm %02ld.%03lds  ", h, m, sint, sfrac );
return(0);
}

/*		julian.c
 *
 * This program calculates Julian day number from calendar
 * date, and the date and day of the week from Julian day.
 * The Julian date is double precision floating point
 * with the origin used by astronomers.  The calendar output
 * converts fractions of a day into hours, minutes, and seconds.
 * There is no year 0.  Enter B.C. years as negative; i.e.,
 * 2 B.C. = -2.
 *
 * The approximate range of dates handled is 4713 B.C. to
 * 54,078 A.D.  This should be adequate for most applications.
 *
 * B.C. dates are calculated by extending the Gregorian sequence
 * of leap years and century years into the past.  This seems
 * the only sensible definition, but I don't know if it is
 * the official one.
 *
 * Note that the astronomical Julian day starts at noon on
 * the previous calendar day.  Thus at midnight in the morning
 * of the present calendar day the Julian date ends in .5;
 * it rolls over to tomorrow at noon today.
 *
 * The month finding algorithm is attributed to Meeus.
 *
 * - Steve Moshier
 */


char *months[12] = {
"January",
"February",
"March",
"April",
"May",
"June",
"July",
"August",
"September",
"October",
"November",
"December"
};

char *days[7] = {
"Sunday",
"Monday",
"Tuesday",
"Wednesday",
"Thursday",
"Friday",
"Saturday"
};

long cyear = 1986L;
extern long cyear;
static int month = 1;
static double day = 1.0;
short yerend = 0;
extern short yerend;

double zgetdate()
{
double J;

/* Get operator to type in a date.
 */
getnum( "Calendar date: Year", &cyear, lngfmt );

if( (cyear > 53994L) || (cyear < -4713L) )
	{
	printf( "Year out of range.\n" );
	goto err;
	}
if( cyear == 0 )
	{
	printf( "There is no year 0.\n" );
err:	J = 0.0;
	goto pdate;
	}

getnum( "Month (1-12)", &month, intfmt);
getnum( "Day.fraction", &day, dblfmt );

/* Find the Julian day. */
J = caltoj(cyear,month,day);
/*printf( "Julian day %.1f\n", J );*/

pdate:
/* Convert back to calendar date. */
/* jtocal( J ); */
return(J);
}



/*     Calculate Julian day from Gregorian calendar date
 */
double caltoj( year, month, day )
long year;
int month;
double day;
{
long y, a, b, c, e, m;
double J;


/* The origin should be chosen to be a century year
 * that is also a leap year.  We pick 4801 B.C.
 */
y = year + 4800;
if( year < 0 )
	{
	y += 1;
	}

/* The following magic arithmetic calculates a sequence
 * whose successive terms differ by the correct number of
 * days per calendar month.  It starts at 122 = March; January
 * and February come after December.
 */
m = month;
if( m <= 2 )
	{
	m += 12;
	y -= 1;
	}
e = (306 * (m+1))/10;

a = y/100;	/* number of centuries */
if( year <= 1582L )
	{
	if( year == 1582L )
		{
		if( month < 10 )
			goto julius;
		if( month > 10)
			goto gregor;
		if( day >= 15 )
			goto gregor;
		}
julius:
	printf( " Julian Calendar assumed.\n" );
	b = -38;
	}
else
	{ /* -number of century years that are not leap years */
gregor:
	b = (a/4) - a;
	}

c = (36525L * y)/100; /* Julian calendar years and leap years */

/* Add up these terms, plus offset from J 0 to 1 Jan 4801 B.C.
 * Also fudge for the 122 days from the month algorithm.
 */
J = b + c + e + day - 32167.5;
return( J );
}




/* Calculate month, day, and year from Julian date */

int jtocal( J )
double J;
{
int month, day;
long year, a, c, d, x, y, jd;
int BC;
double dd;

if( J < 1721425.5 ) /* January 1.0, 1 A.D. */
	BC = 1;
else
	BC = 0;

jd = (long) (J + 0.5); /* round Julian date up to integer */

/* Find the number of Gregorian centuries
 * since March 1, 4801 B.C.
 */
a = (100*jd + 3204500L)/3652425L;

/* Transform to Julian calendar by adding in Gregorian century years
 * that are not leap years.
 * Subtract 97 days to shift origin of JD to March 1.
 * Add 122 days for magic arithmetic algorithm.
 * Add four years to ensure the first leap year is detected.
 */
c = jd + 1486;
if( jd >= 2299160.5 )
	c += a - a/4;
else
	c += 38;
/* Offset 122 days, which is where the magic arithmetic
 * month formula sequence starts (March 1 = 4 * 30.6 = 122.4).
 */
d = (100*c - 12210L)/36525L;
/* Days in that many whole Julian years */
x = (36525L * d)/100L;

/* Find month and day. */
y = ((c-x)*100L)/3061L;
day = (int) (c - x - ((306L*y)/10L));
month = (int) (y - 1);
if( y > 13 )
	month -= 12;

/* Get the year right. */
year = d - 4715;
if( month > 2 )
	year -= 1;

/* Day of the week. */
a = (jd + 1) % 7;

/* Fractional part of day. */
dd = day + J - jd + 0.5;

/* post the year. */
cyear = year;

if( BC )
	{
	year = -year + 1;
	cyear = -year;
	if( prtflg )
		printf( "%ld B.C. ", year );
	}
else
	{
	if( prtflg )
		printf( "%ld ", year );
	}
day = (int) dd;
if( prtflg )
	printf( "%s %d %s", months[month-1], day, days[(int) a] );

/* Flag last or first day of year */
if( ((month == 1) && (day == 1))
	|| ((month == 12) && (day == 31)) )
	yerend = 1;
else
	yerend = 0;

/* Display fraction of calendar day
 * as clock time.
 */
a = (long) dd;
dd = dd - a;
if( prtflg )
	{
	hms( 2.0*PI*dd );
	if( J == TDT )
		printf( "TDT\n" ); /* Indicate Terrestrial Dynamical Time */
	else if( J == UT )
		printf( "UT\n" ); /* Universal Time */
	else
		printf( "\n" );
	}
return(0);
}


/* Reduce x modulo 360 degrees
 */
double mod360(x)
double x;
{
long k;
double y;

k = (long) (x/360.0);
y = x  -  k * 360.0;
while( y < 0.0 )
	y += 360.0;
while( y > 360.0 )
	y -= 360.0;
return(y);
}




/* Reduce x modulo 2 pi
 */
#define TPI (2.0*PI)
double modtp(x)
double x;
{
double y;

y = floor( x/TPI );
y = x - y * TPI;
while( y < 0.0 )
	y += TPI;
while( y >= TPI )
	y -= TPI;
return(y);
}


/* Get operator to type in hours, minutes, and seconds
 */
static int hours = 0;
static int minutes = 0;
static double seconds = 0.0;

double gethms()
{
double t;

getnum( "Time: Hours", &hours, intfmt );
getnum( "Minutes", &minutes, intfmt );
getnum( "Seconds", &seconds, dblfmt );
t = (3600.0*hours + 60.0*minutes + seconds)/86400.0;
return(t);
}

int getnum( msg, num, format )
char *msg;
void *num;
char *format;
{
char s[40];

printf( "%s (", msg );
if( format == strfmt )
	printf( format, (char *) num );
else if( format == dblfmt )
	printf( format, *(double *)num );
else if( format == intfmt )
	printf( format, *(int *)num );
else if( format == lngfmt )
	printf( format, *(long *)num );
else
	printf( "Illegal input format\n"  );
printf( ") ? ");
fgets(s, 40, stdin);
if( s[0] != '\0' )
	sscanf( s, format, num );
return(0);
}



/*
 * Convert change in rectangular coordinatates to change
 * in right ascension and declination.
 * For changes greater than about 0.1 degree, the
 * coordinates are converted directly to R.A. and Dec.
 * and the results subtracted.  For small changes,
 * the change is calculated to first order by differentiating
 *   tan(R.A.) = y/x
 * to obtain
 *    dR.A./cos**2(R.A.) = dy/x  -  y dx/x**2
 * where
 *    cos**2(R.A.)  =  1/(1 + (y/x)**2).
 *
 * The change in declination arcsin(z/R) is
 *   d asin(u) = du/sqrt(1-u**2)
 *   where u = z/R.
 *
 * p0 is the initial object - earth vector and
 * p1 is the vector after motion or aberration.
 *
 */

int deltap( p0, p1, dr, dd )
double p0[], p1[];
double *dr, *dd;
{
double dp[3], A, B, P, Q, x, y, z;
int i;

P = 0.0;
Q = 0.0;
z = 0.0;
for( i=0; i<3; i++ )
	{
	x = p0[i];
	y = p1[i];
	P += x * x;
	Q += y * y;
	y = y - x;
	dp[i] = y;
	z += y*y;
	}

A = sqrt(P);
B = sqrt(Q);

if( (A < 1.e-7) || (B < 1.e-7) || (z/(P+Q)) > 5.e-7 )
	{
	P = zatan2( p0[0], p0[1] );
	Q = zatan2( p1[0], p1[1] );
	Q = Q - P;
	while( Q < -PI )
		Q += 2.0*PI;
	while( Q > PI )
		Q -= 2.0*PI;
	*dr = Q;
	P = asin( p0[2]/A );
	Q = asin( p1[2]/B );
	*dd = Q - P;
	return(0);
	}


x = p0[0];
y = p0[1];
if( x == 0.0 )
	{
	*dr = 1.0e38;
	}
else
	{
	Q = y/x;
	Q = (dp[1]  -  dp[0]*y/x)/(x * (1.0 + Q*Q));
	*dr = Q;
	}

x = p0[2]/A;
P = sqrt( 1.0 - x*x );
*dd = (p1[2]/B - x)/P;
return(0);
}