#include <stdio.h>
#include <stdlib.h>	/* atoi, atof */
#include <math.h>
#include <time.h>
#include <string.h>

#define pi M_PI
#define r2d (180/M_PI)
#define d2r (M_PI/180)

#define DEFLAT 9999
#define DEFLON 9999

/* function declarations */

int timescale(int y, int m, int d);
double eccentric_anomaly(double M, double e);
void sunposition(double d, double* Ms, double* ws, double* RAs, double* Decs);
void moonposition(double d,
					double Ms,
					double ws,
					double* rm,
					double* RAm,
					double* Decm);
double topocentric_correction(double r, double alt_geoc, double lat);
void altitude(double UT,
				double lon,
				double lat,
				double RA,
				double Decl,
				double Ms,
				double ws,
				double* AZ,
				double* Alt);
double modpi(double a);
static void usage();

/*-----------------------------------------------------------------------
Main orbital elements
    N = longitude of the ascending node
    i = inclination to the ecliptic (plane of the Earth's orbit)
    w = argument of perihelion
    a = semi-major axis, or mean distance from Sun
    e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
    M = mean anomaly (0 at perihelion; increases uniformly with time)

Related orbital elements
    w1 = N + w = longitude of perihelion
    L  = M + w1 = mean longitude
    q  = a*(1-e) = perihelion distance
    Q  = a*(1+e) = aphelion distance
    P  = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
    T  = Epoch_of_M - (M(deg)/360) / P  = time of perihelion
    v  = true anomaly (angle between position and perihelion)
    E  = eccentric anomaly
--------------------------------------------------------------------------*/

int main(int argc, char* argv[])
{
	int	i;

	double Ms;			/* Mean Anomaly of the Sun */
	double ws;			/* Argument of perihelion for the Sun */
	double RAs;
	double Decs;
	double RAm;
	double Decm;
	double rm;			/* Moon distance in Earth radii */

	double lon;
	double lat;
	int year;
	int month;
	int day;

	char*	datearg;

	if (argc != 4) {
		usage();
		return 1;
	}

	lon = DEFLON;
	lat = DEFLAT;
	datearg = 0;

	for (i = 1; i < 4; ++i) {
		if (strncmp(argv[i], "lat=", 4) == 0) {
			lat = atof(argv[i] + 4);
		} else if (strncmp(argv[i], "lon=", 4) == 0) {
			lon = atof(argv[i] + 4);
		} else if (strncmp(argv[i], "date=", 5) == 0) {
			datearg = argv[i] + 5;
		}
	}

	if (lon == DEFLON || lat == DEFLAT || datearg == 0) {
		usage();
		return 1;
	}

	year = 2000;
	month = 1;
	day = 1;

	switch (sscanf(datearg, "%d-%d-%d", &year, &month, &day)) {
	case 3:
		break;
	case 2:
		day = 1;
		break;
	case 1:
		month = 1;
		day = 1;
		break;
	case 0:
	default:
		usage();
		return 1;
	}

	if (month < 1 || month > 12) {
		fprintf(stderr, "bad month value\n");
		usage();
		return 1;
	}

	/*---
	note: allow arbitrary day, including forays into next or previous months

	if (day < 1 || day > 31) {
		fprintf(stderr, "bad day value\n");
		usage();
		return 1;
	}
	----*/

	printf("longitude %f\n", lon);
	printf("latitude %f\n", lat);
	printf("yyyy-mm-dd %04d-%02d-%02d\n", year, month, day);
	printf("\n");

	lon *= d2r;
	lat *= d2r;
	double d = timescale(year, month, day);


	printf("            Sun            Moon\n");
	printf("UTC      Alt      AZ    Alt      AZ\n");

	for (int i = 0; i < 24*6; ++i) {
		double UT = 0 + ((double)i)/(24*6);
		sunposition(d + UT, &Ms, &ws, &RAs, &Decs);
		moonposition(d + UT, Ms, ws, &rm, &RAm, &Decm);
		double AZs, Alts;
		double AZm, Altm;
		altitude(UT*24, lon, lat, RAs, Decs, Ms, ws, &AZs, &Alts);
		altitude(UT*24, lon, lat, RAm, Decm, Ms, ws, &AZm, &Altm);
		Altm = topocentric_correction(rm, Altm, lat);
		AZs = modpi(AZs);
		AZm = modpi(AZm);

		/* express the day into the correct month.
			if the day is already in the right month, this code does nothing */
		struct tm stm;
		memset(&stm, 0, sizeof(stm));
		stm.tm_year = year - 1900;
		stm.tm_mon = month - 1;
		stm.tm_mday = day;
		time_t mkt = mktime(&stm);
		struct tm* gmt = gmtime(&mkt);

		printf("%04d-%02d-%02d ", gmt->tm_year + 1900, gmt->tm_mon + 1, gmt->tm_mday);
		printf("%02d:%02d  %5.1f  %6.1f  %5.1f  %6.1f\n",
				i/6,
				10*(i%6),
				Alts*r2d,
				AZs*r2d,
				Altm*r2d,
				AZm*r2d);

	}
}

static void usage()
{
	fprintf(stderr, "usage: altaz lat={lat} lon={lon} date={yyyy-mm-dd}\n");
}

int timescale(int y, int m, int d)
{
	/*-----------
	The time scale used here is a "day number" from 2000 Jan 0.0 TDT, which
	is the same as 1999 Dec 31.0 TDT, i.e. precisely at midnight TDT at the
	start of the last day of this century.  With the modest accuracy we
	strive for here, one can usually disregard the difference between
	TDT (formerly canned ET) and UT.
    d  =  JD - 2451543.5  =  MJD - 51543.0
	We can also compute d directly from the calendar date like this:
    d = 367*Y - (7*(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530
	------------*/

	/* days since jan 01 2000 */
	return 367*y - 7*(y + (m+9)/12)/4 + 275*m/9 + d - 730530;
}

double eccentric_anomaly(double M, double e)
{ 
	/* E and M in radians */
    double E0 = M + e * sin(M) * ( 1.0 + e * cos(M) );
    double E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) );
	int iteration = 0;
	if (e > 0.05) {
		while (fabs(E1 - E0) > 0.001*d2r) {
    		E0 = E1;
    		E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) );
			if (++iteration >= 10) {
				fprintf(stderr, "eccentric anomaly does not converge\n");
				exit(1);
			}
		}
	}
	return E1;
}

void sunposition(double d, double* Ms, double* ws, double* RAs, double* Decs)
{
	/* obliquity of the ecliptic */

    double ecl = 23.4393 - 3.563E-7 * d;
	ecl *= d2r;

	/* Sun orbital elements (angles in degrees)
		N = longitude of the ascending node
    	i = inclination to the ecliptic (plane of the Earth's orbit)
    	w = argument of perihelion
    	a = semi-major axis, or mean distance from Earth
    	e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
    	M = mean anomaly (0 at perigee; increases uniformly with time)
	*/

    double N = 0.0;
    double i = 0.0;
    double w = 282.9404 + 4.70935E-5 * d;
    double a = 1.000000;		  /* AU */
    double e = 0.016709 - 1.151E-9 * d;
    double M = 356.0470 + 0.9856002585 * d;

	N *= d2r;
	i *= d2r;
	w *= d2r;
	M *= d2r;

	/* Excentric anomaly from the mean anomaly */

    double E = M + e*sin(M)*(1.0 + e*cos(M));

	/*
		Compute Sun's distance r and its true anomaly v using
    		xv = r * cos(v) = a * (cos(E) - e)
    		yv = r * sin(v) = a * sqrt(1.0 - e*e) * sin(E)
	*/

    double xv = cos(E) - e;
    double yv = sqrt(1.0 - e*e) * sin(E);
    double v = atan2(yv, xv);
    double r = sqrt(xv*xv + yv*yv);

	/* the Sun's true longitude: */

    double lonsun = v + w;

	/* Convert lonsun,r to ecliptic rectangular geocentric coordinates
	   xs,ys: */

    double xs = r*cos(lonsun);
    double ys = r*sin(lonsun);
	/* zs is zero */

	/* convert this to equatorial, rectangular, geocentric coordinates: */

    double xe = xs;
    double ye = ys*cos(ecl);
    double ze = ys*sin(ecl);

	/* compute the Sun's Right Ascension (RA) and Declination (Dec): */

    double RA  = atan2(ye, xe);
    double Dec = atan2(ze, sqrt(xe*xe+ye*ye));

/*--
	printf("Sun's RA %g\n", RA*r2d);
	printf("Sun's Dec %g\n", Dec*r2d);
--*/

	*Ms = M;
	*ws = w;
	*RAs = RA;
	*Decs = Dec;
}

void moonposition(double d,
					double Ms,
					double ws,
					double* rm,
					double* RAm,
					double* Decm)
{
	/* I. The first order calculation is similar for all planets */

	/* obliquity of the ecliptic */

    double ecl = 23.4393 - 3.563E-7 * d;
	ecl *= d2r;

	/* Moon orbital elements (angles in degrees)
		N = longitude of the ascending node
    	i = inclination to the ecliptic (plane of the Earth's orbit)
    	w = argument of perihelion
    	a = semi-major axis, or mean distance from Earth
    	e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
    	M = mean anomaly (0 at perigee; increases uniformly with time)
	*/

    double N = 125.1228 - 0.0529538083 * d;
    double i = 5.1454;
    double w = 318.0634 + 0.1643573223 * d;
    double a = 60.2666;				  /* Earth radii */
    double e = 0.054900;
    double M = 115.3654 + 13.0649929509 * d;

	N *= d2r;
	i *= d2r;
	w *= d2r;
	M *= d2r;

	/*
		First, compute the eccentric anomaly, E, from M, the mean anomaly,
		and e, the eccentricity.
	*/

	double E = eccentric_anomaly(M, e);

	/*
		Now compute the Moon's distance and true anomaly, using
    	xv = r * cos(v) = a * ( cos(E) - e )
    	yv = r * sin(v) = a * ( sqrt(1.0 - e*e) * sin(E) )
	*/
    double xv = a * ( cos(E) - e );
    double yv = a * ( sqrt(1.0 - e*e) * sin(E) );
    double v = atan2( yv, xv );
    double r = sqrt( xv*xv + yv*yv );

	/* Compute the Moon's position in 3-dimensional space: */

    double xg = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) );
    double yg = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) );
    double zg = r * ( sin(v+w) * sin(i) );

	/*
		As a check one can compute sqrt(xh*xh+yh*yh+zh*zh), which of course
		should equal r (except for small round-off errors).
	*/
/*--
		printf("\tcheck 1: r %f = sqrt(xg*xg+yg*yg+zg*zg) %f\n",
					r,
					sqrt(xg*xg+yg*yg+zg*zg));
--*/

	/*
		This is the geocentric (Earth-centered) position in the ecliptic
		coordinate system.
		If one wishes, one can compute the ecliptic longitude and
		latitude (this must be done if one wishes to correct for
		perturbations, or if one wants to precess the position to a standard
		epoch.
	*/

	/* the ecliptic longitude and latitude: */

    double lonecl = atan2( yg, xg );
    double latecl = atan2( zg, sqrt(xg*xg + yg*yg) );

	/* II. Perturbations of the Moon */

	/*
		If the position of the Moon is computed, and one wishes a better
		accuracy than about 2 degrees, the most important perturbations has
		to be taken into account.  If one wishes 2 arc minute accuracy, all
		the following terms should be accounted for.  If less accuracy is
		needed, some of the smaller terms can be omitted.
	*/

	/* First compute:

    	Ms, Mm             Mean Anomaly of the Sun and the Moon
    	Nm                 Longitude of the Moon's node
    	ws, wm             Argument of perihelion for the Sun and the Moon
	*/

	double Mm = M;
	double Nm = N;
	double wm = w;

    double Ls = Ms + ws;       /* Mean Longitude of the Sun  (Ns=0) */
    double Lm = Mm + wm + Nm;  /* Mean longitude of the Moon */
    double D = Lm - Ls;        /* Mean elongation of the Moon */
    double F = Lm - Nm;        /* Argument of latitude for the Moon */

	/* Add these terms to the Moon's longitude (degrees): */

	lonecl += (
    	-1.274 * sin(Mm - 2*D)          /* the Evection */
    	+0.658 * sin(2*D)               /* the Variation */
    	-0.186 * sin(Ms)                /* the Yearly Equation */
    	-0.059 * sin(2*Mm - 2*D)
    	-0.057 * sin(Mm - 2*D + Ms)
    	+0.053 * sin(Mm + 2*D)
    	+0.046 * sin(2*D - Ms)
    	+0.041 * sin(Mm - Ms)
    	-0.035 * sin(D)                 /* the Parallactic Equation */
    	-0.031 * sin(Mm + Ms)
    	-0.015 * sin(2*F - 2*D)
    	+0.011 * sin(Mm - 4*D))*d2r;

	/* Add these terms to the Moon's latitude (degrees): */

	latecl += (
    	-0.173 * sin(F - 2*D)
    	-0.055 * sin(Mm - F - 2*D)
    	-0.046 * sin(Mm + F - 2*D)
    	+0.033 * sin(F + 2*D)
    	+0.017 * sin(2*Mm + F))*d2r;

	/* Add these terms to the Moon's distance (Earth radii): */

	r +=
    	-0.58 * cos(Mm - 2*D)
    	-0.46 * cos(2*D);

	/*
		All perturbation terms that are smaller than 0.01 degrees in
		longitude or latitude and smaller than 0.1 Earth radii in distance
		have been omitted here.  A few of the largest perturbation terms even
		have their own names!  The Evection (the largest perturbation) was
		discovered already by Ptolemy a few thousand years ago (the Evection
		was one of Ptolemy's epicycles). The Variation and the Yearly
		Equation were both discovered by Tycho Brahe in the 16'th
		century.
	*/

	/*
		Now we have computed the geocentric (Earth-centered) coordinate of
		the Moon, and we have included the most important perturbations.
		We should convert the perturbed lonecl, latecl, r to (perturbed)
		xg, yg, zg:
	*/

    xg = r * cos(lonecl) * cos(latecl);
    yg = r * sin(lonecl) * cos(latecl);
    zg = r               * sin(latecl);

	/*
		We now have the Moons's geocentric (Earth centered) position in
		rectangular, ecliptic coordinates.
	*/

	/* Equatorial coordinates */

	/*
		Let's convert our rectangular, ecliptic coordinates to rectangular,
		equatorial coordinates: simply rotate the y-z-plane by ecl, the angle
		of the obliquity of the ecliptic:
	*/

    double xe = xg;
    double ye = yg * cos(ecl) - zg * sin(ecl);
    double ze = yg * sin(ecl) + zg * cos(ecl);

	/*
		Finally, compute the planet's Right Ascension (RA)
		and Declination (Dec):
	*/

    double RA  = atan2( ye, xe );
    double Dec = atan2( ze, sqrt(xe*xe+ye*ye) );

/*--
	printf("Moon's RA %g\n", RA*r2d);
	printf("Moon's Dec %g\n", Dec*r2d);
--*/

	/* Compute the geocentric distance:
		rg = sqrt(xg*xg+yg*yg+zg*zg) = sqrt(xe*xe+ye*ye+ze*ze)
	*/

/*--
	printf("\tcheck 2: r %f = sqrt(xg*xg+yg*yg+zg*zg) %f\n",
				r,
				sqrt(xg*xg+yg*yg+zg*zg));
	printf("\tcheck 3: r %f = sqrt(xe*xe+ye*ye+ze*ze) %f\n",
				r,
				sqrt(xe*xe+ye*ye+ze*ze));
--*/

	/* Topocentric coordinates */
	/* I ignore this small effect */

	*rm = r;
	*RAm = RA;
	*Decm = Dec;
}

double topocentric_correction(double r, double alt_geoc, double lat)
{
	/* The Moon's topocentric position */

	/*
		The Moon's position, as computed earlier, is geocentric, i.e. as seen
		by an imaginary observer at the center of the Earth.  Real observers
		dwell on the surface of the Earth, though, and they will see a
		different position - the topocentric position.  This position can
		differ by more than one degree from the geocentric position.  To
		compute the topocentric positions, we must add a correction to the
		geocentric position.
	*/

	/*
		Let's start by computing the Moon's parallax, i.e. the apparent
		size of the (equatorial) radius of the Earth, as seen from the Moon:
	*/

    double mpar = asin( 1/r );

	/*
		where r is the Moon's distance in Earth radii.  It's simplest to apply
		the correction in horizontal coordinates (azimuth and altitude):
		within our accuracy aim of 1-2 arc minutes, no correction need to be
		applied to the azimuth.  One need only apply a correction to the
		altitude above the horizon:
	*/

    double alt_topoc = alt_geoc - mpar * cos(alt_geoc);
	return alt_topoc;
}

void altitude(double UT,
			double lon,
			double lat,
			double RA,
			double Decl,
			double Ms,
			double ws,
			double* AZ,
			double* Alt)
{
	/*
		First we compute the Sun's RA and Decl for this moment.  Now we need
		to know our Local Sidereal Time.  We start by computing the sidereal
		time at Greenwich at 00:00 Universal Time, let's call this quantity
		GMST0:
	*/

    double L = Ms + ws;		/* the Sun's mean longitude */
    double GMST0 = L*r2d + 180;

	/*
		Note that we express GMST in degrees here to simplify the
		computations.  360 degrees of course corresponds to 24 hours, i.e.
		each hour corresponds to 15 degrees.
	*/

	/*
		Now we can compute our Local Sidereal Time (LST):
	*/

    double LST = GMST0 + UT*15.0 + lon*r2d;

	/*
		UT is the Universal Time, expressed in hours+decimals, the remaining
		quantities are expressed in degrees.  To convert UT to degrees we
		must multiply it by 15 above.  long is our local longitude in
		degrees, where east longitude counts as positive, and west longitude
		as negative.  (this is according to the geographic standard, and
		the recent astronomical standard; if you prefer to use the older
		astronomical standard where west longitude counts as positive, then
		you must change the '+' in front of 'lon' to a '-' above).
	*/

	/*
		Next let's compute the Sun's Local Hour Angle (LHA), i.e. the angle
		the Earth has turned since the Sun last was in the south:
	*/

    double LHA = LST - RA*r2d;

	/*
		A negative hour angle means the Sun hasn't been in the south yet, this
		day.  The angle -10 degrees is of course the same as 350 degrees,
		i.e. adding or subtracting even multiples of 360 degrees does not
		change the angle.
	*/

	/*
		We also need to know our latitude (lat), where north latitude
		counts as positive and south latitude as negative.  Now we
		can compute the Sun's altitude above the horizon:
	*/

	/*
		We compute sin(h), and then take the arcsine of this to get h, the
		Sun's altitude above the horizon.
	*/

    double sin_h = sin(lat) * sin(Decl) + cos(lat) * cos(Decl) * cos(LHA*d2r);
	double h = asin(sin_h);
	*Alt = h;

	/* formula for the azimuth, reckoned eastward from north: */
	double azim = atan2(-sin(LHA*d2r)*cos(Decl),
						cos(lat)*sin(Decl) - sin(lat)*cos(Decl)*cos(LHA*d2r)); 
	*AZ = azim;
}

double modpi(double a)
{
	while (a < 0) {
		a += 2*pi;
	}
	while (a > 2*pi) {
		a -= 2*pi;
	}
	return a;
}