#include <stdio.h>
#include <math.h>
/*
#include <getopt.h>
*/
#include <unistd.h>

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

double earth_mass = 5.95e24;	/* kg */
double earth_radius = 6370;	/* km */

double moon_mass = 7.35e22;	/* kg */
double moon_dist = 384400;	/* km */

double sun_mass = 1.99e30;	/* kg */
double sun_dist = 150e6;	/* km */

/* 24 hours of data, one record every 10 mins */
#define NPOINTS (24*6)

void usage(FILE*);

int main(int argc, char* argv[])
{
	int i;
	int c;
	bool	moon_only = false;
	bool	sun_only = false;

	while ((c = getopt(argc, argv, "ms?")) != EOF) {
		switch (c) {
		case 'm':
			moon_only = true;
			break;
		case 's':
			sun_only = true;
			break;
		case '?':
			usage(stdout);
			return 0;
		default:
			fprintf(stderr, "tidalvector: bad option\n");
			return 1;
		}
	}

	for (i = 0; i < NPOINTS; ++i) {

		char date[64];
		char time[64];
		double alt_s;
		double az_s;
		double alt_m;
		double az_m;

		if (fscanf(stdin, "%s %s %lf %lf %lf %lf", date, time, &alt_s, &az_s, &alt_m, &az_m) != 6) {
			fprintf(stderr, "tidalvector: bad input data\n");
			return 1;
		}

// fprintf(stderr, "%s %s %f %f %f %f\n", date, time, alt_s, az_s, alt_m, az_m);

		alt_s *= d2r;
		az_s *= d2r;
		alt_m *= d2r;
		az_m *= d2r;

		// from the pt at the center of the Earth, e, there is an 
		// altitude correction, but no azimuth correction

		// altitude correction for center of the earth;
		// alt_m_geoc should be = alt_m when alt_m = pi/2
		// and +mpar when alt_m = 0
		// mpar = moon parallax
		double mpar = earth_radius/moon_dist;
		double alt_m_geoc = alt_m + mpar*cos(alt_m);

		// em = earth-moon
// x is the east-west projection, with east < 0 and west > 0,
		double emx = cos(alt_m_geoc)*cos(3*M_PI/2 - az_m); 
// y is the north-south projection, with north > 0 and south < 0,
		double emy = cos(alt_m_geoc)*sin(az_m - 3*M_PI/2); 
		double emz = sin(alt_m_geoc);
		double emd = 1;

		// pm = point-moon
		// in this coord sys the point p is simply
		// +earth_distance in the z direction
		double pmx = emx;
		double pmy = emy;
		double pmz = emz - mpar;
		double pmd = sqrt(pmx*pmx + pmy*pmy + pmz*pmz);

		double emd3 = emd*emd*emd;
		double pmd3 = pmd*pmd*pmd;

		// am = acceleration due to moon
		// in units of G*moon_mass/moon_dist^2
		double amx = pmx/pmd3 - emx/emd3;
		double amy = pmy/pmd3 - emy/emd3;
		double amz = pmz/pmd3 - emz/emd3;

		// no altitude correction needed for center of the earth;
		// alt_s_geoc = alt_s
		double spar = earth_radius/sun_dist;
		double alt_s_geoc = alt_s + spar*cos(alt_s);

		// es = earth-sun
		double esx = cos(alt_s_geoc)*cos(3*M_PI/2 - az_s); 
		double esy = cos(alt_s_geoc)*sin(az_s - 3*M_PI/2); 
		double esz = sin(alt_s_geoc);
		double esd = 1;

		// ps = point-sun
		double psx = esx;
		double psy = esy;
		double psz = esz - spar;
		double psd = sqrt(psx*psx + psy*psy + psz*psz);

		double esd3 = esd*esd*esd;
		double psd3 = psd*psd*psd;

		// as = acceleration due to sun
		// in units of G*sun_mass/sun_dist^2
		double asx = psx/psd3 - esx/esd3;
		double asy = psy/psd3 - esy/esd3;
		double asz = psz/psd3 - esz/esd3;

		// g = G*mass_earth/(earth_radius)^2
		// t = G*mass_moon/(moon_distance)^2
		// t/g = mass_moon/mass_earth*(earth_radius/moon_dist)^2

		double f = earth_radius/moon_dist;
		f *= f;
		f *= moon_mass/earth_mass;
		f /= 1.12e-7; 

		double s = moon_dist/sun_dist;
		s *= s;
		s *= sun_mass/moon_mass;	// s = 177.807
		if (moon_only) {
			s = 0;
		}

		double m = 1;
		if (sun_only) {
			m = 0;
		}
		double ax = f*(m*amx + s*asx);
		double ay = f*(m*amy + s*asy);
		double az = f*(m*amz + s*asz);
		
		printf("%s %s %f %f %f\n", date, time, ax, ay, az);
	}

	return 0;
}

void usage(FILE* fp)
{
	fprintf(fp, "usage: tidalvector [-m] [-s] < altaz.output\n");
	fprintf(fp, "       -m moon influence only\n");
	fprintf(fp, "       -s sun influence only\n");
}