// Simple Harmonic Oscillator or Pendulum
// numerical integration of the equations of motion

#include <iostream>
#include <fstream>              /* for output file */
#include <cmath>
#include <iomanip>              /* for setw() */
#include <stdlib.h>             /* for exit() */
using namespace std;

double const pi = 4*atan(1);

////////////
// printout routine
//
void p(ostream& ouch,
        double const t,
        double const theta,
        double const thetadot,
        double const g,
        double const steps_per_cycle
){
    ouch << setw(10) << t               << ","
       << setw(15) << theta             << ","
       << setw(15) << thetadot          << ","
       << setw(15) << g                 << ","
       << setw(15) << steps_per_cycle   << ","
       << endl;
}

// The expression "\"" is just too ugly.
// Define a nicer-looking synonym:
string const Q = "\"";

// C++ modulo operator "%" only works for integers
// Generalize to doubles
double mod(double const num, double const denom){
  return num - denom * floor(num / denom);
}

// forward references:
void doit(ostream &xout);

////////////////////////////////////////
// Tiny main program
//
int main() {
  cout << "Data output filenname [or "-"]:  ";
  string out_filename;
  cin >> out_filename;
  if (out_filename == "-") {
    doit(cout);
  }
  else {
    ofstream ouch;
    ouch.open(out_filename.c_str());
    if (!ouch.good()) {
      cerr << "Failed to open output file '" << out_filename << "'" << endl;
      exit(1);
    }
    doit(ouch);
  }
}

///////////////////////////////////////
// function to do most of the work
//
void doit(ostream &xout) {
  cout << "boiler plate for lin. (OR non) Damped Simple Pendulum" << endl;
  cout << "w/ Sync. type impulse and g change" << endl;
  cout << "this version dev. of Verlet" << endl;

  double period = 1;
  double theta = 88e-3;
  if (0) {
    cout << "initial displacement angle in radians 5 deg. ~=88mrad:  ";
    cin >> theta;
  }
  double thetadot = 0;
  double t = 0;
  cout << "steps per cycle (i.e. period / delta_t):  ";
  double steps_per_cycle;
  cin >> steps_per_cycle;
  double deltat = period / steps_per_cycle;
  // cout << deltat << endl;
  double g = 9.82;
  // adjust l to achieve the desired period:
  double l = period*period * g / 4 / pi / pi;
  cout << "period/s: " << period << endl;
  cout << "length/m: " << l << endl;
  cout << "hint: r=0.0126 means Q=500" << endl;
  cout << "damping coefficient (r):  ";
  double r = 0.0126;
  cin >> r;
  cout << "drive amplitude;  hint: 0 or 0.038416:  ";
  double d = 0;
  cin >> d;
  cout << "g increment; hint: 0 or 1e-5:  ";
  double deltag = 0;
  cin >> deltag;
  cout << "number of steps total:  ";
  double stop = 0;
  cin >> stop;
  cout << "Printouts per cycle;  hint: 12:  ";
  double po_per_cycle = 0;
  cin >> po_per_cycle;
  double mod_number = steps_per_cycle / po_per_cycle;

// column headers:
  xout << setw(10) << Q+"t/s"+Q               << ","
       << setw(15) << Q+"theta/rad"+Q         << ","
       << setw(15) << Q+"thetadot"+Q          << ","
       << setw(15) << Q+"g"+Q                 << ","
       << setw(15) << Q+"steps/cycle"+Q       << ","
       << endl;

// print the initial state:
  p(xout, t, theta, thetadot, g, period/deltat);

///////////////////////////////////////////////////
// main loop
  for (int N = 1; N <= stop; N++){
    double da = 0;
    if (theta < 0.05 && theta > 0.0 && thetadot > 0.0){
       da = d;                // sets when driven
    }

// acceleration at beginning of interval:
//  A = d^2 (theta) / dt^2  [F/m]
    double A = - (g/l) * theta  - r*thetadot +da;
    thetadot = A * deltat/2 + thetadot;
    theta = theta + deltat * thetadot;

// recalculate acceleration and end of interval, using new theta:
    A = - (g/l) * theta  - r*thetadot +da;
    thetadot = thetadot + A*deltat/2;
    t = t + deltat;            // keep track of elapsed time

// Print out something every mod_number steps,
// (or as close thereto as we can, if mod_number is not an integer).
// Also print out the very last result
// (which will not be evenly spaced if total duration
// is not an integer multiple of the spacing between printouts).
    if (N == stop || (mod(N + 0.5, mod_number)) < 1) {
      p(xout, t, theta, thetadot, g, period/deltat);
    }
  } // end main loop
}