#include<math.h>
#include<iostream>
#include<fstream>

// This program simulates a particle in a central potential

#define PI 3.141592653589793238462643383279502884197

using namespace std;

class particle
{
  public:
    double _m;
    double _x;
    double _y;
    double _vx;
    double _vy;

    particle(): _m(0), _vx(0), _vy(0), _x(0), _y(0) {};
    particle(double m):_m(m){};

    void init(double x, double y, double vx, double vy) 
      { _x = x; _y = y; _vx= vx; _vy = vy; }
};

void get_force(particle &p, double &fx, double &fy)
{
  double GM=4*PI*PI; // We use astronomical units
  double r = sqrt(p._x*p._x+p._y*p._y);
  double r3 = r * r * r; 
  fx = -GM*p._x/r3;
  fy = -GM*p._y/r3;
}

void 
euler(particle& p, double dt)
{
  double fx, fy;
  get_force(p, fx, fy);
  p._vx += fx*dt;
  p._vy += fy*dt;
  p._x += p._vx*dt;
  p._y += p._vy*dt;
}

void 
verlet(particle& p, double dt)
{
  double fx, fy;
  get_force(p, fx, fy); // before I move to the new position
  p._x += p._vx*dt + 0.5*fx*dt*dt;
  p._y += p._vy*dt + 0.5*fy*dt*dt;
  p._vx += 0.5*fx*dt;
  p._vy += 0.5*fy*dt;
  get_force(p, fx, fy); // after I move to the new position
  p._vx += 0.5*fx*dt;
  p._vy += 0.5*fy*dt;
}

int
main()
{
  double dt;
  double GM=4*PI*PI; // We use astronomical units
  double x0, y0, vx0, vy0;
  double m; // we don't use it, we assume m=1
  particle p;
  int nsteps;

  ofstream fout("central.dat",std::ios::out);
  if(!fout) {
     cerr << "*** ERROR: Could not open file central.dat\n";
     exit(-1);
  }

  cout << "Input initial position (x0):" << endl;
  cin >> x0;
  cout << "Input initial position (y0):" << endl;
  cin >> y0;
  cout << "Input initial velocity (vx0):" << endl;
  cin >> vx0;
  cout << "Input initial velocity (vy0):" << endl;
  cin >> vy0;
  cout << "Input step (dt):" << endl;
  cin >> dt;
  cout << "Input number of steps:" << endl;
  cin >> nsteps;
// Alternatively:
// double tmax;
// cout << "Maximum time (tmax): " << endl;
// cin >> tmax;
// nsteps = tmax/dt;

  p.init(x0, y0, vx0, vy0);
  for(int i = 1; i <= nsteps; i++){
    double r = sqrt(p._x*p._x+p._y*p._y);
    double r3 = r * r * r; 
//    euler(p, dt);
    verlet(p, dt);
    double et = 0.5*(p._vx*p._vx+p._vy*p._vy)-GM/r;
    fout << i*dt << " " << p._x << " " << p._y << " " << p._vx << " " << p._vy << " " << et << endl;
  }

  fout.close();

  return 0;
}