#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 _vx;
    double _fx;
    double _k;

    particle(): _m(1), _vx(0), _x(0), _fx(0), _k(0) {};
    particle(double k):_k(k), _m(1), _x(0), _vx(0), _fx(0) {};

    void init(double x, double vx) 
      { _x = x; _vx= vx; }
/*
    void move(double dt)  // Euler's algorithm
      {
        _x += _vx*dt;
        _vx += _fx*dt/_m;
      }
*/
    void move(double dt)  // Verlet algorithm
      {
        _x += _vx*dt+0.5*_fx*dt*dt/_m;
      }
    void accelerate(double dt)  // Verlet algorithm
      {
        _vx += 0.5*_fx*dt/_m;
      }

};

void
get_forces(int n, particle planet[])
{
///////////////////////////////////////////////////////////////
// For each planet, calculate forces acting on it 
///////////////////////////////////////////////////////////////
    for(int i = 1; i <= n; i++){
      particle &p = planet[i];
      particle &p_left = planet[i-1];
      particle &p_right = planet[i+1];
      p._fx = -p_right._k*(p._x-p_right._x)-p._k*(p._x-p_left._x);
    }
}


int
main()
{
  double dt;
  double x1, vx1;
  double x2, vx2;
  double m1, m2;
  particle planet[NPLANETS+2];
  int nsteps;

  
  for(int i = 1; i <= NPLANETS; i++)
  {
    cout << "Input initial position planet " << i << ":" << endl;
    cin >> planet[i]._x; 
    cout << "Input initial velocity planet " << i << ":" << endl;
    cin >> planet[i]._vx; 
  }
  for(int i = 1; i <= NPLANETS+1; i++)
  {
    cout << "Input spring constant " << i << ":" << endl;
    cin >> planet[i]._k; 
  }

  cout << "Input time step:" << endl;
  cin >> dt; 
// Alternatively:
  double tmax;
  cout << "Maximum time (tmax): " << endl;
  cin >> tmax;
  nsteps = int(double(tmax)/double(dt));

///////////////////////////////////////////////////////////////
// Do loop in time
///////////////////////////////////////////////////////////////
  ofstream fout;
  fout.open("springs.dat",std::ios::out);
  for(int i = 1; i <= nsteps; i++){

//////////////////////////////////////////////////
//  Move the planets
//////////////////////////////////////////////////
    for(int n = 1; n <= NPLANETS; n++){
      particle &p = planet[n];
      get_forces(NPLANETS, planet);
      p.move(dt);
      p.accelerate(dt);
      get_forces(NPLANETS, planet);
      p.accelerate(dt);      

    }
    fout << i*dt << " " << planet[1]._x << " " << planet[2]._x << " " << endl;
  }
  fout.close();
}