// 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() */
#include <time.h>               /* for gmtime_r */
#include <sys/types.h>          /* for stat() */
#include <sys/stat.h>           /* for stat() */
#include <unistd.h>             /* for stat() */

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
){
    ouch << setw(10) << t               << ","
       << setw(15) << theta             << ","
       << setw(15) << thetadot          << ","
       << setw(15) << g                 << ","
       << endl;
}

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

// Same thing, plus a comma:
string const Qc = "\",";

// 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, int const argc, char * const * argv);
string datime(tm const * time);
string datime(time_t const time);

////////////////////////////////////////
// Tiny main program
//
int main(int const argc, char * const * argv) {
  cout << "Data output filenname [or '-']:  ";
  string out_filename;
  cin >> out_filename;
  if (out_filename == "-") {
    doit(cout, argc, argv);
  }
  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, argc, argv);
  }
}

///////////////////////////////////////
// function to do most of the work
//
void doit(ostream &xout, int const argc, char * const * argv) {
  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 length to achieve the desired period:
  double length = period*period * g / 4 / pi / pi;
  cout << "period/s: " << period << endl;
  cout << "length/m: " << length << 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 steps_total = 0;
  cin >> steps_total;
  cout << "Printouts per cycle;  hint: 12:  ";
  double po_per_cycle = 0;
  cin >> po_per_cycle;
  double mod_number = steps_per_cycle / po_per_cycle;

// leave some breadcrumbs in the forest:

// Time that the program was compiled:
  struct stat statbuf;
  string filetime;
  int rslt = stat(argv[0], &statbuf);
  if (rslt) filetime = "?stat failed?";
  else {
    filetime = datime(statbuf.st_mtime);
  }

// current time:
  time_t now;
  time(&now);

// The command line:
  xout << ",,,,," << setw(20) << Q+"cmdline:"+Qc;
     for (int ii = 0; ii < argc; ii++){
       xout << Q+argv[ii]+Qc;
     }
     xout << endl;
  xout << ",,,,," << setw(20) << Q+"compiled/UTC:"+Qc
                  << setw(15) << Q+filetime+Qc << endl;
  xout << ",,,,," << setw(20) << Q+"run/UTC:"+Qc
                  << setw(15) << Q+datime(now)+Qc << endl;
  xout << ",,,,," << setw(20) << Q+"g:"+Qc
                  << setw(15) << g << "," << endl;
  xout << ",,,,," << setw(20) << Q+"length:"+Qc
                  << setw(15) << length << "," << endl;
  xout << ",,,,," << setw(20) << Q+"period:"+Qc
                  << setw(15) << period << "," << endl;
  xout << ",,,,," << setw(20) << Q+"steps_per_cycle:"+Qc
                  << setw(15) << steps_per_cycle << "," << endl;
  xout << ",,,,," << setw(20) << Q+"step_size:"+Qc
                  << setw(15) << deltat << "," << endl;
  xout << ",,,,," << setw(20) << Q+"steps_total:"+Qc
                  << setw(15) << steps_total << "," << endl;
  xout << ",,,,," << setw(20) << Q+"po_per_cycle:"+Qc
                  << setw(15) << po_per_cycle << "," << endl;

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

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

///////////////////////////////////////////////////
// main loop
  for (int N = 1; N <= steps_total; 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/length) * theta  - r*thetadot +da;
    thetadot = A * deltat/2 + thetadot;
    theta = theta + deltat * thetadot;

// recalculate acceleration and end of interval, using new theta:
    A = - (g/length) * 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 == steps_total || (mod(N + 0.5, mod_number)) < 1) {
      p(xout, t, theta, thetadot, g);
    }
  } // end main loop
}

// Convert a "tm" structure into a string
// of a kind that a spreadsheet program will
// recognize as a date.
// It's a UTC date.
string datime(tm const & time){
  int const size(100);
  char text[size];
  int rslt = strftime (text, size, "%Y/%m/%d %T", &time);
  if (rslt) return text;
  else return "?bogus time?";
}

// Same as above, but accets a "time_t" variable as input.
string datime(time_t const time){
  tm fancy;
  gmtime_r(&time, &fancy);  // convert from one format to the other
  return datime(fancy);
}