#include "CalcEphem.h"

/* Moon.c */
double              Moon(double T, double *LAMBDA, double *BETA, double *R,
			 double *AGE);
double              NewMoon(double, double, double);

static double       kepler(double, double);
static double       jd(int, int, int, double);
static double       hour24(double);
static double       angle2pi(double);
static double       angle360(double);

static double
frac(double x)
{
   x -= (int)x;
   return ((x < 0) ? x + 1.0 : x);
}

static int
DayofYear(int year, int month, int day)
{
   return ((int)(jd(year, month, day, 0.0) - jd(year, 1, 0, 0.0)));
}

static int
DayofWeek(int year, int month, int day, char dowstr[])
{
   double              JD, A, Afrac;
   int                 n, iA;

   JD = jd(year, month, day, 0.0);
   A = (JD + 1.5) / 7.0;
   iA = (int)A;
   Afrac = A - (double)iA;
   n = (int)(Afrac * 7.0 + 0.5);
   switch (n)
     {
     case 0:
	strcpy(dowstr, "Sunday");
	break;
     case 1:
	strcpy(dowstr, "Monday");
	break;
     case 2:
	strcpy(dowstr, "Tuesday");
	break;
     case 3:
	strcpy(dowstr, "Wednesday");
	break;
     case 4:
	strcpy(dowstr, "Thursday");
	break;
     case 5:
	strcpy(dowstr, "Friday");
	break;
     case 6:
	strcpy(dowstr, "Saturday");
	break;
     }
   return (n);
}

void
CalcEphem(date, UT, c)
     long int            date;	/* integer containing the date (e.g. 960829) */
     double              UT;	/* Universal Time */
     CTrans             *c;	/* structure containing all the relevent coord trans info */
{
   int                 year, month, day;
   double              TU, TU2, TU3, T0;
   double              varep, varpi;
   double              eccen, epsilon;
   double              days, M, E, nu, lambnew;
   double              r0, earth_sun_distance;
   double              RA, DEC, RA_Moon, DEC_Moon;
   double              TDT;
   double              AGE;
   double              LambdaMoon, BetaMoon, R;
   double              Ta, Tb, Tc;
   double              SinGlat, CosGlat, Tau, lmst, x, y, z;
   double              SinTau, CosTau, SinDec, CosDec;

   c->UT = UT;
   year = (int)(date / 10000);
   month = (int)((date - year * 10000) / 100);
   day = (int)(date - year * 10000 - month * 100);
   c->year = year;
   c->month = month;
   c->day = day;

   c->doy = DayofYear(year, month, day);
   c->dow = DayofWeek(year, month, day, c->dowstr);

   /* 
    *  Compute Greenwich Mean Sidereal Time (gmst)
    *  The TU here is number of Julian centuries
    *  since 2000 January 1.5
    *  From the 1996 astronomical almanac
    */
   TU = (jd(year, month, day, 0.0) - 2451545.0) / 36525.0;
   TU2 = TU * TU;
   TU3 = TU2 * TU;
   T0 = (6.0 + 41.0 / 60.0 + 50.54841 / 3600.0) + 8640184.812866 / 3600.0 * TU
      + 0.093104 / 3600.0 * TU2 - 6.2e-6 / 3600.0 * TU3;
   T0 = hour24(T0);
   c->gmst = hour24(T0 + UT * 1.002737909);

   lmst = 24.0 * frac((c->gmst - c->Glon / 15.0) / 24.0);

   /*
    *
    *   Construct Transformation Matrix from GEI to GSE  systems
    *
    * 
    *   First compute:
    *          mean ecliptic longitude of sun at epoch TU (varep)
    *          elciptic longitude of perigee at epoch TU (varpi)
    *          eccentricity of orbit at epoch TU (eccen)
    *
    *   The TU here is the number of Julian centuries since
    *   1900 January 0.0 (= 2415020.0)
    */
   TDT = UT + 59.0 / 3600.0;
   TU = (jd(year, month, day, TDT) - 2415020.0) / 36525.0;
   varep = (279.6966778 + 36000.76892 * TU + 0.0003025 * TU * TU) * RadPerDeg;
   varpi = (281.2208444 + 1.719175 * TU + 0.000452778 * TU * TU) * RadPerDeg;
   eccen = 0.01675104 - 0.0000418 * TU - 0.000000126 * TU * TU;
   c->eccentricity = eccen;

   /*
    * Compute the Obliquity of the Ecliptic at epoch TU
    * The TU in this formula is the number of Julian
    * centuries since epoch 2000 January 1.5
    */
   TU = (jd(year, month, day, TDT) - jd(2000, 1, 1, 12.0)) / 36525.0;
   epsilon = (23.43929167 - 0.013004166 * TU - 1.6666667e-7 * TU * TU
	      - 5.0277777778e-7 * TU * TU * TU) * RadPerDeg;
   c->epsilon = epsilon;

   /*
    * Compute:
    *          Number of Days since epoch 1990.0 (days)
    *          The Mean Anomaly (M)
    *          The True Anomaly (nu)
    *      The Eccentric Anomaly via Keplers equation (E)
    *          
    *          
    */
   days = jd(year, month, day, TDT) - jd(year, month, day, TDT);
   M = angle2pi(2.0 * M_PI / 365.242191 * days + varep - varpi);
   E = kepler(M, eccen);
   nu = 2.0 * atan(sqrt((1.0 + eccen) / (1.0 - eccen)) * tan(E / 2.0));
   lambnew = angle2pi(nu + varpi);
   c->lambda_sun = lambnew;

   /*
    *  Compute distance from earth to the sun
    */
   r0 = 1.495985e8;		/* in km */
   earth_sun_distance =
      r0 * (1 - eccen * eccen) / (1.0 + eccen * cos(nu)) / 6371.2;
   c->earth_sun_dist = earth_sun_distance;

   /*
    * Compute Right Ascension and Declination of the Sun
    */
   RA =
      angle360(atan2(sin(lambnew) * cos(epsilon), cos(lambnew)) * 180.0 / M_PI);
   DEC = asin(sin(epsilon) * sin(lambnew)) * 180.0 / M_PI;
   c->RA_sun = RA;
   c->DEC_sun = DEC;

   /*
    * Compute Moon Phase and AGE Stuff. The AGE that comes out of Moon()
    * is actually the Phase converted to days. Since AGE is actually defined
    * to be time since last NewMoon, we need to figure out what the JD of the
    * last new moon was. Thats done below....
    */
   TU = (jd(year, month, day, TDT) - 2451545.0) / 36525.0;
   c->MoonPhase = Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
   LambdaMoon *= RadPerDeg;
   BetaMoon *= RadPerDeg;

   RA_Moon =
      angle360(atan2
	       (sin(LambdaMoon) * cos(epsilon) - tan(BetaMoon) * sin(epsilon),
		cos(LambdaMoon)) * DegPerRad);
   DEC_Moon =
      asin(sin(BetaMoon) * cos(epsilon) +
	   cos(BetaMoon) * sin(epsilon) * sin(LambdaMoon)) * DegPerRad;
   c->RA_moon = RA_Moon;
   c->DEC_moon = DEC_Moon;

   /*
    *  Compute Alt/Az coords
    */
   Tau = (15.0 * lmst - RA_Moon) * RadPerDeg;
   CosGlat = cos(c->Glat * RadPerDeg);
   SinGlat = sin(c->Glat * RadPerDeg);
   CosTau = cos(Tau);
   SinTau = sin(Tau);
   SinDec = sin(DEC_Moon * RadPerDeg);
   CosDec = cos(DEC_Moon * RadPerDeg);
   x = CosDec * CosTau * SinGlat - SinDec * CosGlat;
   y = CosDec * SinTau;
   z = CosDec * CosTau * CosGlat + SinDec * SinGlat;
   c->A_moon = DegPerRad * atan2(y, x);
   c->h_moon = DegPerRad * asin(z);
   c->Visible = (c->h_moon < 0.0) ? 0 : 1;

   /*
    * Compute accurate AGE of the Moon
    */
   Tb = TU - AGE / 36525.0;	/* should be very close to minimum */
   Ta = Tb - 0.4 / 36525.0;
   Tc = Tb + 0.4 / 36525.0;
   c->MoonAge = (TU - NewMoon(Ta, Tb, Tc)) * 36525.0;

   /*
    * Compute Earth-Moon distance
    */
   c->EarthMoonDistance = R;
}

static double
kepler(M, e)
     double              M, e;
{
   int                 n = 0;
   double              E, Eold, eps = 1.0e-8;

   E = M + e * sin(M);
   do
     {
	Eold = E;
	E = Eold + (M - Eold + e * sin(Eold)) / (1.0 - e * cos(Eold));
	++n;
     }
   while ((fabs(E - Eold) > eps) && (n < 100));
   return (E);
}

/*
 *  Compute the Julian Day number for the given date.
 *  Julian Date is the number of days since noon of Jan 1 4713 B.C.
 */
static double
jd(ny, nm, nd, UT)
     int                 ny, nm, nd;
     double              UT;
{
   double              A, B, C, D, JD, day;

   day = nd + UT / 24.0;

   if ((nm == 1) || (nm == 2))
     {
	ny = ny - 1;
	nm = nm + 12;
     }

   if (((double)ny + nm / 12.0 + day / 365.25) >=
       (1582.0 + 10.0 / 12.0 + 15.0 / 365.25))
     {
	A = ((int)(ny / 100.0));
	B = 2.0 - A + (int)(A / 4.0);
     }
   else
     {
	B = 0.0;
     }

   if (ny < 0.0)
     {
	C = (int)((365.25 * (double)ny) - 0.75);
     }
   else
     {
	C = (int)(365.25 * (double)ny);
     }

   D = (int)(30.6001 * (double)(nm + 1));

   JD = B + C + D + day + 1720994.5;
   return (JD);
}

static double
hour24(hour)
     double              hour;
{
   int                 n;

   if (hour < 0.0)
     {
	n = (int)(hour / 24.0) - 1;
	return (hour - n * 24.0);
     }
   else if (hour > 24.0)
     {
	n = (int)(hour / 24.0);
	return (hour - n * 24.0);
     }
   else
     {
	return (hour);
     }
}

static double
angle2pi(angle)
     double              angle;
{
   int                 n;
   double              a;

   a = 2.0 * M_PI;

   if (angle < 0.0)
     {
	n = (int)(angle / a) - 1;
	return (angle - n * a);
     }
   else if (angle > a)
     {
	n = (int)(angle / a);
	return (angle - n * a);
     }
   else
     {
	return (angle);
     }
}

static double
angle360(angle)
     double              angle;
{
   int                 n;

   if (angle < 0.0)
     {
	n = (int)(angle / 360.0) - 1;
	return (angle - n * 360.0);
     }
   else if (angle > 360.0)
     {
	n = (int)(angle / 360.0);
	return (angle - n * 360.0);
     }
   else
     {
	return (angle);
     }
}

#if 0
static void
Radec_to_Cart(ra, dec, r)
     double              ra, dec;	/* RA and DEC */
     Vector             *r;	/* returns corresponding cartesian unit vector */
{
   /*
    *  Convert ra/dec from degrees to radians
    */
   ra *= RadPerDeg;
   dec *= RadPerDeg;

   /*
    *  Compute cartesian coordinates (in GEI)
    */
   r->x = cos(dec) * cos(ra);
   r->y = cos(dec) * sin(ra);
   r->z = sin(dec);
}

int
LeapYear(year)
     int                 year;
{
   if ((year % 100 == 0) && (year % 400 != 0))
      return (0);
   else if (year % 4 == 0)
      return (1);
   else
      return (0);
}
#endif