////////////////////////////
// Anharmonic Oscillator
// includes harmonic oscillator as a special case
//
// numerical integration the equations of motin

using namespace std;
#include <iostream>     // for cerr, cout
#include <iomanip>      // for setw()
#include <cmath>        // for M_PI and cos()

// acceleration = - (g/L) * x   // also knowns as torque
// assume units such that g/L equals 1

#define accel(x) (-(x))

const double pi(M_PI);

void dumpit(double const period, int const Nst, int const steps_per_cycle,
        double const t, double const x, double const v, double const resid){
  double denom(Nst);
  if (resid == 0. && denom == 0) denom = 1;  // 0/0 is harmless;  treat as 0/1
  cout << fixed      << setw(10) << t/period
       << scientific << setw(15) << x
       << scientific << setw(15) << v
       << scientific << setw(15) << (x*x + v*v -1.)
       << scientific << setw(15) << sqrt(resid / denom)
       << scientific << setw(15) << resid / denom
       << scientific << setw(15) << steps_per_cycle
       << endl;
}

int main(){
// principal configuration options:
  int verbosity(0);
  double Cmax(3);              // desired number of cycles
  int steps_per_cycle(1000);            // integrator steps per cycle

// various quantities derived from the configuration:
  double period(2*pi);
  double dt(period/steps_per_cycle);       // time per integrator step

// initial conditions:
  double t(0);
  double x(1);                  // the position, also known as theta
  double v(0);                  // the velocity
  double resid(0);
  double Nst(0);                // total number of integrator steps taken

// prime the pump:
  double w(v + accel(x) * dt / 2.);    // advanced velocity
  double oldw(w);

// column headers:
  cout  << setw(10) << "t/period"
        << setw(15) << "x    "
        << setw(15) << "v    "
        << setw(15) << "     Delta E  "
        << setw(15) << "   rms resid  "
        << setw(15) << "    ms resid  "
        << setw(15) << "   steps/cycle"
        << endl;

// outer loop over all cycles
  for (int cc = 0; cc < Cmax; cc++) {

// inner loop;  for explanation, see
// https://www.av8n.com/physics/symplectic-integrator.htm#sec-verlet
//
// At the beg of the step, we have x[ii],    t[ii],    and w[ii]
// At the end of the step, we have x[1+ii],  t[1+ii],  and w[1+ii]
    for (int ii = 0; ii < steps_per_cycle; ii++){
      if (verbosity) {
        dumpit(period, Nst, steps_per_cycle, t, x, v, resid);
      }
      oldw = w;                   // leave some breadcrumbs in the forest
      x += w * dt;
      w += accel(x) * dt;
      t += dt;
      Nst++;
      double delta = x - cos(t);
      resid += delta*delta;
// v is needed for printouts; otherwise wouldn't be needed:
      v = (oldw + w) / 2.;
    }

// At this point,
// oldw is w[Nsteps-1]
// x    is x[Nsteps]
// w    is w[Nsteps]

    v = (oldw + w) / 2.;
    if (cc == Cmax-1 || !verbosity)
      dumpit(period, Nst, steps_per_cycle, t, x, v, resid);
  } // end outer loop over cycles
}