#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>

typedef struct AberrationCoefficient
{
  int L2;
  int L3;
  int L4;
  int L5;
  int L6;
  int Ldash;

  int xsin;
  int xcos;
  int ysin;
  int ycos;
} AberrationCoefficient;

const AberrationCoefficient g_AberrationCoefficients[] =
{ 
  //L2   L3   L4  L5  L6  Ldash xsin xcos ysin ycos
  {  0,  0,   0,  1,  0,  0,  715,   0, -16, -715 }, // ysin actually zsin
  {  0,  0,   0,  0,  1,  0,  159,   0,  -6, -164 }, // ysin actually zsin
  {  0,  0,   0,  2,  0,  0,   33, -10, -11,  -32 },
  {  0,  2,   0,  -1, 0,  0,   31,   1,   1,  -30 },
  {  0,  3,  -8,  3,  0,  0,    8, -28,  27,    9 },
  {  0,  5,  -8,  3,  0,  0,    8, -28, -27,   -9 },
  {  0,  1,   0,  -2, 0,  0,   16,   0,   1,   16 },
  {  0,  1,   0,  1,  0,  0,   11,  -1,  -2,  -11 },
  {  0,  1,   0,  -1, 0,  0,  -11,  -2,  -2,   10 },
  {  0,  3,   0,  -2, 0,  0,  -10,   0,   0,   10 },
  {  0,  0,   0,  0,  2,  0,    0,  -9,  -9,    0 },
  //{  1, -66  }, // F, zcos
  //{  1,  17,  -17 }, // L7, xsin, ycos
  //{  1,  16,  -17 }, // L8, xsin, ycos
};

//  L7  L8  xsin ycos   zsin   zcos
//{  1,  0,   17, -16,     0,    -7    },
//{  0,  1,   16, -15,    -3,    -6    },

//////////////////////////////// Implementation ///////////////////////////////

double EarthVelocity(double JD)
{
  //double X, Y, Z;
  struct Coord3d { double X, Y; } velocity;
  double Z, tmp;

  double T = (JD - 2451545) / 36525;
  double L3 = 1.7534703 + 628.3075849 * T;
  double L4 = 6.2034809 + 334.0612431 * T;
  double L5 = 0.5995465 + 52.9690965 * T;
  double L6 = 0.8740168 + 21.3299095 * T;
  double Ldash = 3.8103444 + 8399.6847337 * T;
  //double l7 = 5.4812939 + 7.4781599 * t;
  //double l8 = 5.3118863 + 3.8133036 * t;
  double F = 1.6279052 + 8433.4661601 * T;

  //   L3       xsin   xsint    xcos  xcost    ysin   ysint ycos     ycost zsin   zsint zcos    zcost
  //{  1,   -1719914,     -2,    -25,     0,     25, -13,  1578089,  156,    10,  32, 684185, -358    },
  //{  2,       6434,    141,  28007,  -107,  25697, -95,    -5904, -130, 11141, -48,  -2559,  -55    },
  //{  3,        486,     -5,   -236,    -4,   -216,  -4,     -446,    5,   -94,   0,   -193,    0    },
    velocity.X  = -(1719914 +   2 * T) * sin(L3);
    velocity.X +=       -25            * cos(L3);
    velocity.Y  =        27            * sin(L3);
    velocity.Y +=  1720021 * cos(L3);
    Z  =   -66       * cos(F);

    velocity.X +=  ( 6434 + 141 * T) * sin(2*L3);
    velocity.X +=  (28007 - 107 * T) * cos(2*L3);
    velocity.Y +=  (28007 - 106 * T) * sin(2*L3);
    velocity.Y += -( 6435 + 141 * T) * cos(2*L3);

    velocity.X +=  (486 - 5 * T) * sin(3*L3);
    velocity.X += -(236 + 4 * T) * cos(3*L3);
    velocity.Y += -(236 + 4 * T) * sin(3*L3);
    velocity.Y += -(499 - 5 * T) * cos(3*L3);

    velocity.X += 7 * sin(4*L3);
    velocity.X -= 8 * cos(4*L3);
    velocity.Y -= 9 * sin(4*L3);
    velocity.Y -= 7 * cos(4*L3);

    velocity.X += (715 * sin(Ldash));
    velocity.Y -= (715 * cos(Ldash));

  int i, nAberrationCoefficients = sizeof(g_AberrationCoefficients) / sizeof(AberrationCoefficient);
  for (i=0; i<nAberrationCoefficients; i++)
  {
    double Argument = g_AberrationCoefficients[i].L3*L3 + 
                      g_AberrationCoefficients[i].L4*L4 +
                      g_AberrationCoefficients[i].L5*L5 + 
                      g_AberrationCoefficients[i].L6*L6; 
    velocity.X += (g_AberrationCoefficients[i].xsin * sin(Argument));
    velocity.X += (g_AberrationCoefficients[i].xcos * cos(Argument));
    velocity.Y += (g_AberrationCoefficients[i].ysin * sin(Argument));
    velocity.Y += (g_AberrationCoefficients[i].ycos * cos(Argument));
  }

  velocity.X *= 1.49597870691;
  velocity.Y *= 1.49597870691;
  Z *= 1.49597870691;

    velocity.X *= 0x1.7b425ed097b42p-7;
    velocity.Y *= 0x1.7b425ed097b42p-7;
    Z *= 0x1.7b425ed097b42p-7;
    printf("Velocity: X: %.5f, Y: %.5f, Z: %.5f\n", velocity.X, velocity.Y, Z);
    return sqrt((velocity.X*velocity.X + velocity.Y*velocity.Y));
}

int main(int argc, char **argv) {
    double JD = strtod(argv[1], NULL)+0.5;
    double RR;

    setvbuf(stdout, NULL, _IONBF, 0);

    RR = EarthVelocity(JD);
    printf("Range Rate: %.7f\n", RR);
    return 0;
}