// planet.cc     by George Saklatvala 
//
// Shows how to implement dynamics
// in an elegant advanced C++ way.
//
// The structure xyvect, declared at
// the top of the file, uses arithmetic operator overloading
//
// Usage:
//   planet > tmp
//   set size ratio -1 ; plot 'tmp' u 2:3 w linesp 3 4, 'tmp' u 2:3:4:5 w vector lt 5

#include<iostream>
#include<fstream>

#include<cfloat>
#include<cstdlib>
#include<cmath>

using namespace std;

struct xyvect{

  double x,y;
  xyvect(): x(), y() {}
  xyvect(double xx, double yy): x(xx), y(yy) {}
  double normsq() const {return x*x + y*y;}
  double norm() const {return sqrt(normsq());}
  double ang() const {return atan2(y,x);}

  xyvect& operator+=(const xyvect& r)
  {
    x += r.x;
    y += r.y;
    return *this;
  }

  xyvect& operator-=(const xyvect& r)
  {
    x -= r.x;
    y -= r.y;
    return *this;
  }

  xyvect& operator*=(double r)
  {
    x *= r;
    y *= r;
    return *this;
  }

  xyvect& operator/=(double r)
  {
    x /= r;
    y /= r;
    return *this;
  }

};

xyvect operator* (double d, const xyvect& v)
{
  xyvect a(v);
  a*=d;
  return a;
}

xyvect operator* (const xyvect& v, double d)
{
  xyvect a(v);
  a*=d;
  return a;
}

xyvect operator/ (const xyvect& v, double d)
{
  xyvect a(v);
  a/=d;
  return a;
}

struct impact{}; //exception class

// This is where the planet bit starts....

const double G = 1; // gravitational constant
const double M = 1; // mass of sun

class particle{

private:

  const double m;
  xyvect x;
  xyvect v;
  double t;
  xyvect a;

  double r() const {return x.norm();}
  double theta() const {return x.ang();}
  double J() const {return m*(x.x*v.y-x.y*v.x);}
  double T() const {return m*v.normsq();}
  double V() const {return -G*M*m/r();}
  double E() const {return T() + V();}

public:

  particle(double mass, xyvect x0, xyvect v0, double t0=0): 
    m(mass), x(x0), v(v0), t(t0), a(){}

  void x_inc(double dt) // step x in the direction of v.
  {
    x += v * dt;
    if(r()<2*dt*v.norm() || r() < DBL_EPSILON)throw impact();
  }

  void v_inc(double dt) // step v in the direction of the acceleration.
  {
    v += a * dt;
  }

  void t_inc(double dt)
  {
    t += dt;
  }

  void a_calc() // find the acceleration
  {
    a = -G*M*x/pow(r(),3);
  }

  ostream& write_state(ostream& of) const
  {
    of<<t<<" "
      <<x.x<<" "<<x.y<<" "
      <<v.x<<" "<<v.y<<" "
      <<r()<<" "<<theta()<<" "
      <<J()<<" "<<T()<<" "
      <<V()<<" "<<E()<<endl;
    return of;  
  }
};

void euler_step(particle& p, double dt)
{
  p.a_calc();
  p.x_inc(dt);
  p.v_inc(dt);
  p.t_inc(dt);
}

void leapfrog_step(particle& p, double dt)
{
  p.x_inc(0.5*dt);
  p.a_calc();
  p.v_inc(dt);
  p.x_inc(0.5*dt);
  p.t_inc(dt);
}

void simulate(particle& p, void(*step)(particle&, double),
	      double t, int n, int frames, ostream& of)
{ // number of steps made is n (total duration is t)
  // number of frames written out is frames
  const double dt = t / n;
  for(int i=0;i<frames-1;++i){
    p.write_state(of);
    for(int j=0;j<n/frames;++j)
      step(p,dt);
  }
}

int main()
{
  double time = 11.0; // time to simulate
  int n = 1000;      // total number of steps to make
  int frames = 100; // number of frames to record
  double m = 1;
  const xyvect x0(1.5,2.0);
  const xyvect v0(0.4,0);
  particle p(m,x0,v0);
  try{
    simulate(p,leapfrog_step,time,n,frames,cout);
    // edit the above command to change the step style
  }catch(impact&){
    cout << "Collision occurred\n";
  }
  return 0; 
}