#include "moon.h"


/*********************************************
 *Declaration of functions/methods/macros used
 *********************************************/

#define SQR(x)	(x) * (x)
void Sid_Time(double JC, long secs, double *ST);
void moon_geoHP(double JD, long secs, double *ra, double *dec, double *dist);
void geo_topo(double LST, double *dist, double *ra, double *dec);
void eq_horiz(double dec, double h, double *az, double *alt);
void series(double T, double p[], double *u, double *v, double *w);


/************************************************************
 * This method call all other in this source file.  Takes inputs
 *      based on moon and calls functions to make corrections
 *		and calculate data.
 * Inputs: Modified Julian Date (Date since         )
 *		   Time in seconds
 * Outputs: Right Ascention and Declination, Alltitude and Azimuth, hour angle
 ************************************************************/
void moon (long MJD, long secs, double *dec, double *ra, double *alt, double *az, double *hAngle)

{
	double JD, JC, LST, dist;
	LST = dist = 0;

	//Modified Julian Date to Julian Date  
	JD = 2400000.5 + MJD;

	//time argument: Jul. Cent. from the epoch 2000.0
	JC = (JD - 2451545) / 36525;

	// get the local siderial time
	Sid_Time(JC, secs, &LST);
	LST *= 15*DEGRAD;

	// geocentric coordinates 
	moon_geoHP(JD, secs, ra, dec, &dist);

	// topocentric correction 
	geo_topo(LST, &dist, ra, dec);

	// horizontal coordinates (alt, az);
	*hAngle = LST - *ra;
	eq_horiz(*dec, *hAngle, az, alt);
	
	*ra  *= RADDEG;
	if (*ra<0) *ra += 360.0;
	*dec *= RADDEG;
	*az  *= RADDEG;
	if (*az<0) *az += 360.0;
	*alt *= RADDEG;
	*hAngle *= RADDEG/15.;
	if (*hAngle>24.0) *hAngle -= 24.0;
	if (*hAngle< 0.0) *hAngle += 24.0;
}



/*****************************************************
 * topocentric correction at your geographioc position
 * (geographic latitude)        input: Local siderial
 * return: distance of the moon in		ra and dec
 *****************************************************/
void geo_topo(double LST, double *dist, double *ra, double *dec)
{
	double     		x,y,z;
	x = *dist*cos(*dec)*cos(*ra) - cos(OGL)*cos(LST);
	y = *dist*cos(*dec)*sin(*ra) - cos(OGL)*sin(LST);
	z = *dist*sin(*dec) - sin(OGL);
	*dist = sqrt(x*x+y*y+z*z);
	*ra = atan2(y,x);
	*dec = asin(z/(*dist));
}


/*****************************************************
 * Sidereal day is one revolution w/ respect to stars.
 * Sidereal day is shorter than solar day.  To compensate,
 * add 3min 57sec to each hour
 * ( 0.002737962962963 = 3min 56.56sec / 24h )
 *****************************************************/
void Sid_Time(double JC, long secs, double *ST)
{
	double		h,hours;
	
	h = 0.2790572732621 + 100.0021390378 * JC;
	h += 0.0000010775925925 * JC * JC;
	h -= floor(h);
	*ST = 24.0 * h;
	hours = secs/3600.0;

	*ST = *ST + hours + hours * 0.002737962962963;
	*ST = *ST + OL/15;

	if ( *ST > 24.0 )
		*ST -= 24.0;
}


/**********************************************
 * Converts from equatorial to horizontal coordinates
 * Input: declinatin and angle
 * Outputs: azimuth and altitude
 **********************************************/
void eq_horiz(double dec, double h, double *az, double *alt)
{
	double         f,g;

	*alt = asin(cos(h)*cos(dec)*cos(OGL) + sin(dec)*sin(OGL));
	f=sin(h)*cos(dec)/cos(*alt);
	g=(cos(h)*cos(dec)*sin(OGL) - sin(dec)*cos(OGL))/cos(*alt);
	*az = atan2(f,g);

	return;
}


/************************************************************
 * Find the geocentric coordinates of the moon
 * Inputs: Julian Date and time in seconds
 * Outputs: Right Ascention, Declination and Distance of moon
 * Notes: Time argument T is Jul Cent from E 1900
 *        Time argument t is Jul date and time from
 *		  E 2000 Jan 1.5	== Jul Date 2451545.0
 ************************************************************/
void moon_geoHP(double JD, long secs, double *ra, double *dec, double *dist)
{
	double		T,t;
	double		p[20],h;
	double		u,v,w;
	long		i;
	u = v = w = 0;

	t =  	JD  - 2451545. + secs/(24 * 3600.);
	T =  	t / 36525 + 1;

	p[0]  = 0;
	p[1]	= 0.606434 + 0.03660110129*t;
	p[2]	= 0.374897 + 0.03629164709*t;
	p[3]	= 0.259091 + 0.03674819520*t;
	p[4]	= 0.827362 + 0.03386319198*t;
	p[5]	= 0.347343 - 0.00014709391*t;
	p[7]	= 0.779072 + 0.00273790931*t;
	p[8]	= 0.993126 + 0.00273777850*t;
	p[9]	= 0;
	p[10]	= 0;
	p[11]	= 0;
	p[12] = 0.505498 + 0.00445046867*t;
	p[13] = 0.140023 + 0.00445036173*t;
	p[14]	= 0;
	p[15]	= 0;
	p[16] = 0.053856 + 0.00145561327*t;
	p[17]	= 0;
	p[18]	= 0;
	p[19] = 0.056531 + 0.00023080893*t;

	for(i = 0; i < 20; i++) p[i] = (p[i] - ceil(p[i])) * TWOPI;

	series(T, p, &u, &v, &w);

	h 	= sqrt(u - v * v);
	*ra 	= p[1] + asin(w / h);
	h 	= sqrt(u);
	*dec 	= asin(v / h);
	*dist = 60.40974 * h;
}


/**********************************************************
 * This is the insane funcion that makes all of the corrections,
 *    giving the final position of the moon.
 * Inputs: Time, and an array that holds................
 * Outputs: (x, y, z) coordinates for the moon?????????????
 **********************************************************/
void series(double T, double p[], double *u, double *v, double *w)
{
	*u = 
		1.00000
		-0.10828 *   cos(p[2])
		-0.01880 *   cos(p[2]-2*p[4])
		-0.01479 *   cos(2*p[4])
		+0.00181 *   cos(2*p[2]-2*p[4])
		-0.00147 *   cos(2*p[2])
		-0.00105 *   cos(2*p[4]-p[8])
		-0.00075 *   cos(p[2]-2*p[4]+p[8])
		-0.00067 *   cos(p[2]-p[8])
		+0.00057 *   cos(p[4])
		+0.00055 *   cos(p[2]+p[8])
		-0.00046 *   cos(p[2]+2*p[4])
		+0.00041 *   cos(p[2]-2*p[3])
		+0.00024 *   cos(p[8])
		+0.00017 *   cos(2*p[4]+p[8])
		+0.00013 *   cos(p[2]-2*p[4]-p[8])
		-0.00010 *   cos(p[2]-4*p[4])
		-0.00009 *   cos(p[4]+p[8])
		+0.00007 *   cos(2*p[2]-2*p[4]+p[8])
		+0.00006 *   cos(3*p[2]-2*p[4])
		+0.00006 *   cos(2*p[3]-2*p[4])
		-0.00005 *   cos(2*p[4]-2*p[8])
		-0.00005 *   cos(2*p[2]-4*p[4])
		+0.00005 *   cos(p[2]+2*p[3]-2*p[4])
		-0.00005 *   cos(p[2]-p[4])
		-0.00004 *   cos(p[2]+2*p[4]-p[8])
		-0.00004 *   cos(3*p[2])
		-0.00003 *   cos(p[2]-4*p[4]+p[8])
		-0.00003 *   cos(2*p[2]-2*p[3])
		-0.00003 *   cos(2*p[3]);

	*v =
		 0.39558 *   sin(p[3]+p[5])
		+0.08200 *   sin(p[3])
		+0.03257 *   sin(p[2]-p[3]-p[5])
		+0.01092 *   sin(p[2]+p[3]+p[5])
		+0.00666 *   sin(p[2]-p[3])
		-0.00644 *   sin(p[2]+p[3]-2*p[4]+p[5])
		-0.00331 *   sin(p[3]-2*p[4]+p[5])
		-0.00304 *   sin(p[3]-2*p[4])
		-0.00240 *   sin(p[2]-p[3]-2*p[4]-p[5])
		+0.00226 *   sin(p[2]+p[3])
		-0.00108 *   sin(p[2]+p[3]-2*p[4])
		-0.00079 *   sin(p[3]-p[5])
		+0.00078 *   sin(p[3]+2*p[4]+p[5])
		+0.00066 *   sin(p[3]+p[5]-p[8])
		-0.00062 *   sin(p[3]+p[5]+p[8])
		-0.00050 *   sin(p[2]-p[3]-2*p[4])
		+0.00045 *   sin(2*p[2]+p[3]+p[5])
		-0.00031 *   sin(2*p[2]+p[3]-2*p[4]+p[5])
		-0.00027 *   sin(p[2]+p[3]-2*p[4]+p[5]+p[8])
		-0.00024 *   sin(p[3]-2*p[4]+p[5]+p[8])
		-0.00021 *T* sin(p[3]+p[5])
		+0.00018 *   sin(p[3]-p[4]+p[5])
		+0.00016 *   sin(p[3]+2*p[4])
		+0.00016 *   sin(p[2]-p[3]-p[5]-p[8])
		-0.00016 *   sin(2*p[2]-p[3]-p[5])
		-0.00015 *   sin(p[3]-2*p[4]+p[8])
		-0.00012 *   sin(p[2]-p[3]-2*p[4]-p[5]+p[8])
		-0.00011 *   sin(p[2]-p[3]-p[5]+p[8])
		+0.00009 *   sin(p[2]+p[3]+p[5]-p[8])
		+0.00009 *   sin(2*p[2]+p[3])
		+0.00008 *   sin(2*p[2]-p[3])
		+0.00008 *   sin(p[2]+p[3]+2*p[4]+p[5])
		-0.00008 *   sin(3*p[3]-2*p[4]+p[5])
		+0.00007 *   sin(p[2]-p[3]+2*p[4])
		-0.00007 *   sin(2*p[2]-p[3]-2*p[4]-p[5])
		-0.00007 *   sin(p[2]+p[3]+p[5]+p[8])
		-0.00006 *   sin(p[3]+p[4]+p[5])
		+0.00006 *   sin(p[3]-2*p[4]-p[8])
		+0.00006 *   sin(p[2]-p[3]+p[5])
		+0.00006 *   sin(p[3]+2*p[4]+p[5]-p[8])
		-0.00005 *   sin(p[2]+p[3]-2*p[4]+p[8])
		-0.00004 *   sin(2*p[2]+p[3]-2*p[4])
		+0.00004 *   sin(p[2]-3*p[3]-p[5])
		+0.00004 *   sin(p[2]-p[3]-p[8])
		-0.00003 *   sin(p[2]-p[3]+p[8])
		+0.00003 *   sin(p[3]-p[4])
   		+0.00003 *   sin(p[3]-2*p[4]+p[5]-p[8])
   		-0.00003 *   sin(p[3]-2*p[4]-p[5])
		+0.00003 *   sin(p[3]-p[8])
		-0.00003 *   sin(p[3]-p[4]+p[5]-p[8])
		-0.00002 *   sin(p[2]-p[3]-2*p[4]+p[8])
		-0.00002 *   sin(p[3]+p[8])
		+0.00002 *   sin(p[2]+p[3]-p[4]+p[5])
		-0.00002 *   sin(p[2]+p[3]-p[5])
		+0.00002 *   sin(3*p[2]+p[3]+p[5])
		-0.00002 *   sin(2*p[2]-p[3]-4*p[4]-p[5])
		+0.00002 *   sin(p[2]-p[3]-2*p[4]-p[5]-p[8])
		-0.00002 *T* sin(p[2]-p[3]-p[5])
		-0.00002 *   sin(p[2]-p[3]-4*p[4]-p[5])
		-0.00002 *   sin(p[2]+p[3]-4*p[4])
		-0.00002 *   sin(2*p[2]-p[3]-2*p[4])
		+0.00002 *   sin(p[2]+p[3]+2*p[4])
		+0.00002 *   sin(p[2]+p[3]-p[8])
		+0.00003 * sin(p[2]+p[3]-2*p[4]+p[5]-p[8]);
		
	*w =
   		 0.10478 *   sin(p[2])
		-0.04105 *   sin(2*p[3]+2*p[5])
		-0.02130 *   sin(p[2]-2*p[4])
		-0.01779 *   sin(2*p[3]+p[5])
		+0.01774 *   sin(p[5])
		+0.00987 *   sin(2*p[4])
		-0.00338 *   sin(p[2]-2*p[3]-2*p[5])
		-0.00309 *   sin(p[8])
		-0.00190 *   sin(2*p[3])
		-0.00144 *   sin(p[2]+p[5])
		-0.00144 *   sin(p[2]-2*p[3]-p[5])
		-0.00113 *   sin(p[2]+2*p[3]+2*p[5])
		-0.00094 *   sin(p[2]-2*p[4]+p[8])
		-0.00092 *   sin(2*p[2]-2*p[4])
		+0.00071 *   sin(2*p[4]-p[8])
		+0.00070 *   sin(2*p[2])
		+0.00067 *   sin(p[2]+2*p[3]-2*p[4]+2*p[5])
		+0.00066 *   sin(2*p[3]-2*p[4]+p[5])
		-0.00066 *   sin(2*p[4]+p[5])
		+0.00061 *   sin(p[2]-p[8])
		-0.00058 *   sin(p[4])
		-0.00049 *   sin(p[2]+2*p[3]+p[5])
		-0.00049 *   sin(p[2]-p[5])
		-0.00042 *   sin(p[2]+p[8])
		+0.00034 *   sin(2*p[3]-2*p[4]+2*p[5])
		-0.00026 *   sin(2*p[3]-2*p[4])
		+0.00025 *   sin(p[2]-2*p[3]-2*p[4]-2*p[5])
		+0.00024 *   sin(p[2]-2*p[3])
		+0.00023 *   sin(p[2]+2*p[3]-2*p[4]+p[5])
  		+0.00023 *   sin(p[2]-2*p[4]-p[5])
  		+0.00019 *   sin(p[2]+2*p[4])
		+0.00012 *   sin(p[2]-2*p[4]-p[8])
		+0.00011 *   sin(p[2]-2*p[4]+p[5])
		+0.00011 *   sin(p[2]-2*p[3]-2*p[4]-p[5])
		-0.00010 *   sin(2*p[4]+p[8])
		+0.00009 *   sin(p[2]-p[4])
		+0.00008 *   sin(p[4]+p[8])
		-0.00008 *   sin(2*p[3]+2*p[4]+2*p[5])
		-0.00008 *   sin(2*p[5])
		-0.00007 *   sin(2*p[3]+2*p[5]-p[8])
		+0.00006 *   sin(2*p[3]+2*p[5]+p[8])
		-0.00005 *   sin(p[2]+2*p[3])
		+0.00005 *   sin(3*p[2])
		-0.00005 *   sin(p[2]+16*p[7]-18*p[12])
		-0.00005 *   sin(2*p[2]+2*p[3]+2*p[5])
		+0.00004 *T* sin(2*p[3]+2*p[5])
		+0.00004 *   cos(p[2]+16*p[7]-18*p[12])
		-0.00004 *   sin(p[2]-2*p[3]+2*p[4])
		-0.00004 *   sin(p[2]-4*p[4])
		-0.00004 *   sin(3*p[2]-2*p[4])
		-0.00004 *   sin(2*p[3]+2*p[4]+p[5])
		-0.00004 *   sin(2*p[4]-p[5])
		-0.00003 *   sin(2*p[8])
		-0.00003 *   sin(p[2]-2*p[4]+2*p[8])
		+0.00003 *   sin(2*p[3]-2*p[4]+p[5]+p[8])
		-0.00003 *   sin(2*p[4]+p[5]-p[8])
		+0.00003 *   sin(2*p[2]+2*p[3]-2*p[4]+2*p[5])
		+0.00003 *   sin(2*p[4]-2*p[8])
		-0.00003 *   sin(2*p[2]-2*p[4]+p[8])
		+0.00003 *   sin(p[2]+2*p[3]-2*p[4]+2*p[5]+p[8])
		-0.00003 *   sin(2*p[2]-4*p[4])
		+0.00002 *   sin(2*p[3]-2*p[4]+2*p[5]+p[8])
		-0.00002 *   sin(2*p[2]+2*p[3]+p[5])
		-0.00002 *   sin(2*p[2]-p[5])
		+0.00002 *T* cos(p[2]+16*p[7]-18*p[12])
		+0.00002 *   sin(4*p[4])
		-0.00002 *   sin(2*p[3]-p[4]+2*p[5])
		-0.00002 *   sin(p[2]+2*p[3]-2*p[4])
		-0.00002 *   sin(2*p[2]+p[5])
		-0.00002 *   sin(2*p[2]-2*p[3]-p[5])
		+0.00002 *   sin(p[2]+2*p[4]-p[8])
		+0.00002 *   sin(2*p[2]-p[8])
		-0.00002 *   sin(p[2]-4*p[4]+p[8])
		+0.00002 *T* sin(p[2]+16*p[7]-18*p[12])
		-0.00002 *   sin(p[2]-2*p[3]-2*p[5]-p[8])
		+0.00002 *   sin(2*p[2]-2*p[3]-2*p[5])
		-0.00002 *   sin(p[2]+2*p[4]+p[5])
		-0.00002 *   sin(p[2]-2*p[3]+2*p[4]-p[5]);
}