// Calculate the 1-d spring pendulum energy and trajectory
// Isaac Lenton
// August 2017

#include <fstream>

const double mass = 1.0;             // Mass of pendulum (kg)
const double acceleration = 0.0;     // Gravitational acceleration (m/s/s)
const double springConst = 1.0;      // Spring constant (kg/s/s)
const double x0 = -1.0;              // Initial position (m)

const double dt = 1.0e-3;            // Simulation step size
const unsigned outputRate = 100;     // Number of steps between output cycles
const unsigned outputLength = 100;   // Length of simulation

double potential(double x) {
  return -mass*acceleration*x + 0.5*springConst*x*x; }
double eqmotion(double x) { return mass*acceleration - springConst*x; }
double kinetic(double v) { return 0.5*mass*v*v; }

int main(int argc, char** argv) {

  double x = x0, v = 0.0, t = 0.0;

  std::ofstream fp("output.dat");
  fp << "# t\tx\tv\tEk\tEp\tEtotal (should be const)\n";
  fp << t << '\t' << x << '\t' << v << '\t' << kinetic(v) << '\t'
      << potential(x) << '\t' << potential(x)+kinetic(v) << '\n';

  for (unsigned i = 0; i < outputLength; ++i) {
    for (unsigned k = 0; k < outputRate; ++k) {
      double oldv = v;
      v += eqmotion(x)/mass*dt;
      x += oldv*dt;
    }
    t += outputRate*dt;

    fp << t << '\t' << x << '\t' << v << '\t' << kinetic(v) << '\t'
        << potential(x) << '\t' << potential(x)+kinetic(v) << '\n';
  }

  return 0;
}