#include <iostream>
#include <cmath>

class equation {
public:
  double dxdt(double x, double y, double z, double vx, double vy, double vz) { return vx; }
  double dydt(double x, double y, double z, double vx, double vy, double vz) { return vy; }
  double dzdt(double x, double y, double z, double vx, double vy, double vz) { return vz; }
  double dvxdt(double x, double y, double z, double vx, double vy, double vz) {
    double r = std::sqrt(x * x + y * y + z * z);
    return -x / (r * r * r);
  }
  double dvydt(double x, double y, double z, double vx, double vy, double vz) {
    double r = std::sqrt(x * x + y * y + z * z);
    return -y / (r * r * r);
  }
  double dvzdt(double x, double y, double z, double vx, double vy, double vz) {
    double r = std::sqrt(x * x + y * y + z * z);
    return -z / (r * r * r);
  }
};

struct mass {
  double x, y, z, vx, vy, vz;
  mass() {}
  mass(double x, double y, double z, double vx, double vy, double vz) : x(x), y(y), z(z), vx(vx), vy(vy), vz(vz) {}
  friend double area(mass m1, mass m2) {
    double t1, t2, t3;
    t1 = m1.y * m2.z - m1.z * m2.y;
    t2 = m1.z * m2.x - m1.x * m2.z;
    t3 = m1.x * m2.y - m1.y * m2.x;
    return 0.5 * std::sqrt(t1 * t1 + t2 * t2 + t3 * t3);
  }
  friend std::ostream &operator<<(std::ostream &stream, mass m) {
    stream << "(x, y, z, vx, vy, vz) = (" << m.x << ", " << m.y << ", " << m.z << ", ";
    stream                                << m.vx << ", " << m.vy << ", " << m.vz << ")";
    return stream;
  }
};

static int const N = 30000;
static double const T  = 3.0;

class RungeKutta4 {
public:
  void operator()(mass &m, equation f) {
    double k1x, k1y, k1z, k1vx, k1vy, k1vz;
    double k2x, k2y, k2z, k2vx, k2vy, k2vz;
    double k3x, k3y, k3z, k3vx, k3vy, k3vz;
    double k4x, k4y, k4z, k4vx, k4vy, k4vz;
    double h = T / N;

    k1x = h * f.dxdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
    k1y = h * f.dydt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
    k1z = h * f.dzdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
    k1vx = h * f.dvxdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
    k1vy = h * f.dvydt(m.x, m.y, m.z, m.vx, m.vy, m.vz);
    k1vz = h * f.dvzdt(m.x, m.y, m.z, m.vx, m.vy, m.vz);

    k2x = h * f.dxdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                     m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
    k2y = h * f.dydt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                     m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
    k2z = h * f.dzdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                     m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
    k2vx = h * f.dvxdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                       m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
    k2vy = h * f.dvydt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                       m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);
    k2vz = h * f.dvzdt(m.x + k1x / 2.0, m.y + k1y / 2.0, m.z + k1z / 2.0,
                       m.vx + k1vx / 2.0, m.vy + k1vy / 2.0, m.vz + k1vz / 2.0);

    k3x = h * f.dxdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                     m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
    k3y = h * f.dydt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                     m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
    k3z = h * f.dzdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                     m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
    k3vx = h * f.dvxdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                       m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
    k3vy = h * f.dvydt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                       m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);
    k3vz = h * f.dvzdt(m.x + k2x / 2.0, m.y + k2y / 2.0, m.z + k2z,
                       m.vx + k2vx / 2.0, m.vy + k2vy / 2.0, m.vz + k2vz / 2.0);

    k4x = h * f.dxdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
    k4y = h * f.dydt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
    k4z = h * f.dzdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
    k4vx = h * f.dvxdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
    k4vy = h * f.dvydt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);
    k4vz = h * f.dvzdt(m.x + k3x, m.y + k3y, m.z + k3z, m.vx + k3vx, m.vy + k3vy, m.vz + k3vz);

    m.x = m.x + (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
    m.y = m.y + (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;
    m.z = m.z + (k1z + 2.0 * k2z + 2.0 * k3z + k4z) / 6.0;
    m.vx = m.vx + (k1vx + 2.0 * k2vx + 2.0 * k3vx + k4vx) / 6.0;
    m.vy = m.vy + (k1vy + 2.0 * k2vy + 2.0 * k3vy + k4vy) / 6.0;
    m.vz = m.vz + (k1vz + 2.0 * k2vz + 2.0 * k3vz + k4vz) / 6.0;
  }
};

int const K = 2000;
int main() {
  RungeKutta4 e;

  equation f;
  mass m(0.0, 0.70710678, 0.70710678, 1.0, 0.0, 0.0), m0;
  std::cout << "t = 0.0 : " << m << std::endl;
  int i = K + 1;
  double s = 0.0;
  for (int k = 0; k < T * N; k++) {
    m0 = m;
    e(/*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 */