#include <stdio.h>
#include <math.h>

#define pi 3.14159265358979323
#define dtor (pi/180.0)



struct orbel
{
   double O;
   double i;
   double w;
   double a;
   double e;
   double M;
   double Odot;  //time unit = 1 day
   double idot;
   double wdot;
   double adot;
   double edot;
   double Mdot;
};

double getd(int y, int m, int d, int h)
{
  return  (long)367 * y - 7 * (y + (m + 9) / 12) / 4
	   + 275 * m / 9
	   + d - 730531.5
	   + h / 24.0;
}
void getplanetels(struct orbel *p);
void getplanetphases(double *x, double *y, double *z, double *vx, double *vy, double *vz, 
                     double d, double dt);

void calcpos(double *x, double *y, double *z,
	     double O, double i, double w, double a, double e, double M);
void calcposel(double *x, double *y, double *z, struct orbel *elements, double d);
void calcphaseel(double *x, double *y, double *z, double *vx, double *vy, double *vz,
		struct orbel *elements, double d, double dt);


void getplanetels(struct orbel *p)
{
  p[1].O = 48.3313;   p[1].Odot = 3.24587E-5;
  p[1].i =  7.0047;   p[1].idot = 5.00E-8;
  p[1].w = 29.1241;   p[1].wdot = 1.01444E-5;
  p[1].a = 0.387098;  p[1].adot = 0;
  p[1].e = 0.205635;  p[1].edot = 5.59E-10;
  p[1].M = 168.6562;  p[1].Mdot = 4.0923344368;

  p[2].O =  76.6799;  p[2].Odot =  2.46590E-5;
  p[2].i = 3.3946;    p[2].idot =  2.75E-8;
  p[2].w =  54.8910;  p[2].wdot =  1.38374E-5;
  p[2].a = 0.723330;  p[2].adot =  0;
  p[2].e = 0.006773;  p[2].edot = -1.302E-9;
  p[2].M =  48.0052;  p[2].Mdot =  1.6021302244;

  p[3].O = 0;         p[3].Odot =  0;
  p[3].i = 0.0;       p[3].idot =  0;
  p[3].w = 102.9472;  p[3].wdot =  4.70935E-5;
  p[3].a = 1.000000;  p[3].adot =  0;
  p[3].e = 0.016709;  p[3].edot = -1.151E-9;
  p[3].M = 356.0470;  p[3].Mdot =  0.9856002585;
  //O=0, w= 282.9404, M = 356.0470

  p[4].O =  49.5574;  p[4].Odot =  2.11081E-5;
  p[4].i = 1.8497;    p[4].idot = -1.78E-8;
  p[4].w = 286.5016;  p[4].wdot =  2.92961E-5;
  p[4].a = 1.523688;  p[4].adot =  0;
  p[4].e = 0.093405;  p[4].edot =  2.516E-9;
  p[4].M =  18.6021;  p[4].Mdot =  0.5240207766;

  p[5].O = 100.4542;  p[5].Odot =  2.76854E-5;
  p[5].i = 1.3030;    p[5].idot = -1.557E-7;
  p[5].w = 273.8777;  p[5].wdot =  1.64505E-5;
  p[5].a = 5.20256;   p[5].adot =  0;
  p[5].e = 0.048498;  p[5].edot =  4.469E-9;
  p[5].M =  19.8950;  p[5].Mdot =  0.0830853001;

  p[6].O = 113.6634;  p[6].Odot =  2.38980E-5;
  p[6].i = 2.4886;    p[6].idot = -1.081E-7;
  p[6].w = 339.3939;  p[6].wdot =  2.97661E-5;
  p[6].a = 9.55475;   p[6].adot =  0;
  p[6].e = 0.055546;  p[6].edot = -9.499E-9;
  p[6].M = 316.9670;  p[6].Mdot =  0.0334442282;

  p[7].O =  74.0005;  p[7].Odot =  1.3978E-5;
  p[7].i = 0.7733;    p[7].idot =  1.9E-8;
  p[7].w =  96.6612;  p[7].wdot =  3.0565E-5;
  p[7].a = 19.18171;  p[7].adot = -1.55E-8;
  p[7].e = 0.047318;  p[7].edot =  7.45E-9;
  p[7].M = 142.5905;  p[7].Mdot =  0.011725806;

  p[8].O = 131.7806;  p[8].Odot =  3.0173E-5;
  p[8].i = 1.7700;    p[8].idot = -2.55E-7;
  p[8].w = 272.8461;  p[8].wdot = -6.027E-6;
  p[8].a = 30.05826;  p[8].adot =  3.313E-8;
  p[8].e = 0.008606;  p[8].edot =  2.15E-9;
  p[8].M = 260.2471;  p[8].Mdot =  0.005995147;
}

void getplanetphases(double *x, double *y, double *z, double *vx, double *vy, double *vz,
		     double d, double dt)
{
  int i;
  struct orbel p[9];

  getplanetels(p);

  for (i = 1; i <= 8; i++)
  {
    calcphaseel(&x[i], &y[i], &z[i], &vx[i], &vy[i], &vz[i], &p[i], d, dt); 
  }
}


/*void getmoonpos(double *x, double *y, double *z, double d)
{
    double O, i, w, a, e, M;
    double xe, ye, ze;
    O = 125.1228 - 0.0529538083 * d
    i = 5.1454
    w = 318.0634 + 0.1643573223 * d
    a = 60.2666 * 149598000.0 / 6378.137; 
    e = 0.054900
    M = 115.3654 + 13.0649929509 * d
    calcpos(x,y,z,O,i,w,a,e,M);
    getearthpos(xe,ye,ze,d);
    *x += *xe;   
    *y += *ye;   
    *z += *ze;   
}*/



void calcpos(double *x, double *y, double *z,
	     double O, double i, double w, double a, double e, double M)
{
  double E;  //eccentric anomaly
  double v;  //true anomaly
  double r;  //perihelion distance

  double Eold = 0.0;
  double xv, yv;

  O *= dtor;
  i *= dtor;
  w *= dtor;
  M *= dtor;

  E = M + e * sin(M) * ( 1.0 + e * cos(M) );
  if (e > 0.05)
  {
    do
    {
      Eold = E;
      E -= ( E - e * sin(E) - M ) / ( 1 - e * cos(E) );
    } while (E - Eold >= 0.00001);
 }

  xv = a * ( cos(E) - e );
  yv = a * ( sqrt(1.0 - e*e) * sin(E) );

  v = atan2( yv, xv );
  r = sqrt( xv*xv + yv*yv );

  *x = r * (cos(O) * cos(v+w) - sin(O) * sin(v+w) * cos(i));
  *y = r * (sin(O) * cos(v+w) + cos(O) * sin(v+w) * cos(i));
  *z = r * (sin(v+w) * sin(i));

  return;
}

void calcposel(double *x, double *y, double *z, struct orbel *elements, double d)
{
  calcpos(x,y,z,
	  elements->O + elements->Odot*d,
	  elements->i + elements->idot*d,
	  elements->w + elements->wdot*d,
	  elements->a + elements->adot*d,
	  elements->e + elements->edot*d,     //why minus?  who knows!
	  elements->M + elements->Mdot*d);
}


void calcphaseel(double *x, double *y, double *z, double *vx, double *vy, double *vz, struct orbel *elements, double d, double dt)
{
  //If dt = 1 timestep this will give exact answer with Euler's method for t=1...
  double x1,y1,z1;
  calcposel(x,y,z,elements,d);
  calcposel(&x1,&y1,&z1,elements,d+dt);
  *vx = (x1-*x)/dt;
  *vy = (y1-*y)/dt;
  *vz = (z1-*z)/dt;
}