#include <iostream>
#include <cmath>

using namespace std;

#include"const.h"

// Time integration parameter
#define TMAX 50
#define dt 0.01

int main (int argc, const char **argv)
{
  // Define the variables we need
  double x[3],v[3],a[3],t;
  double mass = msun / m_unit;   // Mass of central body in internal units (will be 1 in the chosen units)
  double m = mearth / m_unit;    // Mass of planet (only needed for energy calculation)
  
  // Lets output at the beginning some useful information ... 
  cout << "# L = " << l_unit << endl;
  cout << "# T = " << t_unit << endl;
  cout << "# M = " << m_unit << endl;
  cout << "# V = " << v_unit << endl;
  cout << "# G = " << G << endl;

// Setup the earth orbit (1 au, 1 year period) for the state which we start with
  x[0] = au / l_unit;
  v[1] = 2.0 * M_PI;
  // Eccentric orbit: v[1] = 0.8 * 2.0 * M_PI;
   
  // Here is another interesting way to simplify initialisation of several variables in C to the same value
  x[1] = x[2] = v[0] = v[2] = a[0] = a[1] = a[2] = t = 0.0;

  // This is loop which is done until integration of t reaches our defined TMAX
  while(t < TMAX)
    {
      t += dt;
      
      double r = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);

      // Calculate acceleration
      for (int i=0;i<3;i++)                     // Calculate the gravitational acceleration for the 3 spacial components
	a[i] = -G * mass / (r * r * r) * x[i];

      // Update position
      for (int i=0;i<3;i++)                     // Update the 3 position components
	x[i] += v[i] * dt;

      // Update velocities
      for (int i=0;i<3;i++)                     // Update the 3 velocity components
	v[i] += a[i] * dt;

      r = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);  // r now already exists, so we do not need to define it again
      double v2 = v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
      double Epot = m * mass * G / r;            // Potential energy
      double Ekin = 0.5 * m * v2;                // Kinetic energy

      // Lets output time , position and total energy 
      cout << t << " " << x[0] << " " << x[1] << " " << x[2] << " " << Ekin << " " << Epot << endl;

    }

  return(0);
}

/* How to get it running:
  1) g++ example_simple.cc
  2) ./a.out >& orbit.dat
  3) gnuplot
     gnuplot> plot 'orbit.dat' u ($2):($3)
*/