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

// This program simulates a mini-solar system with two planets

#define PI 3.141592653589793238462643383279502884197
#define NPLANETS 2

using namespace std;

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

    particle(): _m(0), _vx(0), _vy(0), _x(0), _y(0), _fx(0), _fy(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 move(double dt)  // Euler's algorithm
      {
        _x += _vx*dt;
        _y += _vy*dt;
        _vx += _fx*dt/_m;
        _vy += _fy*dt/_m;
      }

};

int
main()
{
  double y, dt;
  double x1, y1, vx1, vy1;
  double x2, y2, vx2, vy2;
  double m1, m2;
  double GM=4*PI*PI; // We use astronomical units
  particle planet[NPLANETS];
  int nsteps;

  ofstream fout;
  fout.open("planets.dat",std::ios::out);
  if(!fout) {
     cerr << "*** ERROR: Could not open file planets.dat\n";
     exit(-1);
  }
  
  // We use the Solar Mass as our unit of mass (Solar Units)
  cout << "Input mass of planet 1 (m1/M):" << endl;
  cin >> m1;
  cout << "Input initial position (x1):" << endl;
  cin >> x1;
  cout << "Input initial position (y1):" << endl;
  cin >> y1;
  cout << "Input initial velocity (vx1):" << endl;
  cin >> vx1;
  cout << "Input initial velocity (vy1):" << endl;
  cin >> vy1;
  cout << "Input mass of planet 2 (m2/M):" << endl;
  cin >> m2;
  cout << "Input initial position (x2):" << endl;
  cin >> x2;
  cout << "Input initial position (y2):" << endl;
  cin >> y2;
  cout << "Input initial velocity (vx2):" << endl;
  cin >> vx2;
  cout << "Input initial velocity (vy2):" << endl;
  cin >> vy2;
  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;


// Initializae all particles
  planet[0]._m = m1;
  planet[0].init(x1, y1, vx1, vy1);
  planet[1]._m = m2;
  planet[1].init(x2, y2, vx2, vy2);
///////////////////////////////////////////////////////////////
// Do loop in time
///////////////////////////////////////////////////////////////
  for(int i = 1; i <= nsteps; i++){
///////////////////////////////////////////////////////////////
// For each planet, calculate forces acting on it 
///////////////////////////////////////////////////////////////
    for(int n = 0; n < NPLANETS; n++){
      particle &p = planet[n];
      p._fx = 0.;
      p._fy = 0.;
      // Compute the interaction forces with the other planet(s)
      double f12x = 0, f12y = 0;
      for(int j = 0; j < NPLANETS; j++){
        if(j != n) {
          double x12 = p._x-planet[j]._x;
          double y12 = p._y-planet[j]._y;
          double r12 = sqrt(x12*x12 + y12*y12);
          r12 = r12 * r12 * r12;
          f12x += GM*p._m*planet[j]._m*x12/r12; 
          f12y += GM*p._m*planet[j]._m*y12/r12; 
        }
      }

      double r3 = sqrt(p._x*p._x+p._y*p._y);
      r3 = r3 * r3 * r3; 
      p._fx = -GM*p._m*p._x/r3 - f12x;
      p._fy = -GM*p._m*p._y/r3 - f12y;
    }
    cout << i*dt << " ";
//////////////////////////////////////////////////
//  Move the planets
//////////////////////////////////////////////////
    for(int n = 0; n < NPLANETS; n++){
      particle &p = planet[n];
      p.move(dt);
      fout << p._x << " " << p._y << " " << p._vx << " " << p._vy << " ";
    }
    fout << endl;
  }

  fout.close();

  return 0;
}