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

#define NSTEPS 1000000
float r2(float x, float y);
float r(float x, float y);

int main(int argc, void *argv[])
{
  int i;
  float x[2], v[2];
  float dt = 1e-4;
  float a = 1, e = 0.8, M=1, G = 1, Eo, E, Lo, L;
  float k1x[2], k2x[2], k3x[2], k4x[2];
  float k1v[2], k2v[2], k3v[2], k4v[2];
  float xnew[2], f[2], vnew[2];
  FILE *out, *outrk4;
  
  outrk4 = fopen("outrk4.dat", "w");
  out = fopen("out.dat", "w");
  
  x[0] = a*(1+e);
  x[1] = 0;
  v[0] = 0;
  v[1] = sqrt((G*M/a)*(1-e)/(1+e));
  Eo = 0.5 * fabs(v[1]);
  Lo = x[0] * v[1];

  for(i=0; i<NSTEPS; i++)
    {
      k1x[0] = v[0];
      k1x[1] = v[1];
      
      f[0] = -G*M/r2(x[0],x[1])*x[0]/r(x[0],x[1]);
      f[1] = -G*M/r2(x[0],x[1])*x[1]/r(x[0],x[1]);
      
      k1v[0] = f[0];
      k1v[1] = f[1];
      
      xnew[0] = x[0] + 0.5 * dt * (v[0] + 0.5 * dt * k1v[0]);
      xnew[1] = x[1] + 0.5 * dt * (v[1] + 0.5 * dt * k1v[1]);
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);

      k2v[0] = f[0];
      k2v[1] = f[1];
      
      xnew[0] = x[0] + 0.5 * dt * k1x[0];
      xnew[1] = x[1] + 0.5 * dt * k1x[1];
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);
      
      k2x[0] = v[0] + 0.5 * dt * f[0];
      k2x[1] = v[1] + 0.5 * dt * f[1];
      
      xnew[0] = x[0] + 0.5 * dt * (v[0] + 0.5 * dt * k2v[0]);
      xnew[1] = x[1] + 0.5 * dt * (v[1] + 0.5 * dt * k2v[1]);
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);

      k3v[0] = f[0];
      k3v[1] = f[1];
      
      xnew[0] = x[0] + 0.5 * dt * k2x[0];
      xnew[1] = x[1] + 0.5 * dt * k2x[1];
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);
      
      k3x[0] = v[0] + 0.5 * dt * f[0];
      k3x[1] = v[1] + 0.5 * dt * f[1];
      
      xnew[0] = x[0] + dt * (v[0] + dt * k3v[0]);
      xnew[1] = x[1] + dt * (v[1] + dt * k3v[1]);
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);

      k4v[0] = f[0];
      k4v[1] = f[1];
      
      xnew[0] = x[0] + dt * k3x[0];
      xnew[1] = x[1] + dt * k3x[1];
      
      f[0] = -G*M/r2(xnew[0],xnew[1])*xnew[0]/r(xnew[0],xnew[1]);
      f[1] = -G*M/r2(xnew[0],xnew[1])*xnew[1]/r(xnew[0],xnew[1]);
      
      k4x[0] = v[0] + dt * f[0];
      k4x[1] = v[1] + dt * f[1];

      v[0] = v[0] + dt/6. * (k1v[0] + 2 * k2v[0] + 2 * k3v[0] + k4v[0]);
      v[1] = v[1] + dt/6. * (k1v[1] + 2 * k2v[1] + 2 * k3v[1] + k4v[1]);

      x[0] = x[0] + dt/6. * (k1x[0] + 2 * k2x[0] + 2 * k3x[0] + k4x[0]);
      x[1] = x[1] + dt/6. * (k1x[1] + 2 * k2x[1] + 2 * k3x[1] + k4x[1]);

      E = 0.5 * sqrt(v[0] * v[0] + v[1] * v[1]);
      L = x[0] * v[1] - x[1] * v[0];

      fprintf(outrk4, "%g %g \n", (E-Eo)/E, (L-Lo)/L);
      fprintf(out, "%g %g %g %g\n", x[0], x[1], v[0], v[1]);
      
    }

  fclose(out);  
  fclose(outrk4);

  return 0;
}

float r2(float x, float y)
{
  return x*x + y*y;
}

float r(float x, float y)
{
  return sqrt(x*x + y*y);
}