#include <iostream>
#include <cmath>

class equation {
public:
  double dxdt(double x, double y, double vx, double vy) { return vx; }
  double dydt(double x, double y, double vx, double vy) { return vy; }
  double dvxdt(double x, double y, double vx, double vy) {
    double r = sqrt(x * x + y * y);
    return -x / (r * r * r);
  }
  double dvydt(double x, double y, double u, double v) {
    double r = sqrt(x * x + y * y);
    return -y / (r * r * r);
  }
};

struct mass {
  double x, y, vx, vy;
  mass() {}
  mass(double x, double y, double vx, double vy) : x(x), y(y), vx(vx), vy(vy) {}
  friend double area(mass m1, mass m2) {
    return 0.5 * std::fabs(m1.x * m2.y - m2.x * m1.y);
  }
  friend std::ostream &operator<<(std::ostream &stream, mass m) {
    stream.width(9);
    stream.precision(6);
    stream << "(x, y, vx, vy) = (" << m.x << ", " << m.y << ", ";
    stream                         << m.vx << ", " << m.vy << ")";
    return stream;
  }
};

static int const N = 10000;
static double const T  = 10.0;

class RungeKutta2 {
  double k1x, k2x, k1y, k2y, k1vx, k2vx, k1vy, k2vy;
public:
  void operator()(mass &m, equation f) {
    double h = T / N;
    k1x = f.dxdt(m.x, m.y, m.vx, m.vy);
    k1y = f.dydt(m.x, m.y, m.vx, m.vy);
    k1vx = f.dvxdt(m.x, m.y, m.vx, m.vy);
    k1vy = f.dvydt(m.x, m.y, m.vx, m.vy);

    k2x = f.dxdt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);
    k2y = f.dydt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vy, m.vy + h * k1vy);
    k2vx = f.dvxdt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);
    k2vy = f.dvydt(m.x + h * k1x, m.y + h * k1y, m.vx + h * k1vx, m.vy + h * k1vy);

    m.x = m.x + h * (k1x + k2x) / 2.0;
    m.y = m.y + h * (k1y + k2y) / 2.0;
    m.vx = m.vx + h * (k1vx + k2vx) / 2.0;
    m.vy = m.vy + h * (k1vy + k2vy) / 2.0;
  }
};

int const K = 1000;
int main() {
  RungeKutta2 euler;

  equation f;
  mass m(0.0, 1.0, 1.2, 0.0), m0;
  int i = K + 1;
  double s = 0.0;
  for (int k = 0; k < T * N; k++) {
    m0 = m;
    euler(/*ref*/m, f);
    s += area(m0, m);
    if (--i == 0) {
      std::cout << "t = " << (double)k / N << ':' << m << ": s = " << s << std::endl;
      i = K;
      s = 0.0;
    }
  }
  return 0;
}
/* end */