// Simple Harmonic Oscillator or Pendulum driven synnetrically (deadbeat escapement)
// numerical integration of its equation 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(ofstream& ouch,
        double const t,
        double const theta,
        double const thetadot,
        double const deltag,
        double const da
){
    ouch << setw(10) << t               << ","
       << setw(15) << theta             << ","
       << setw(15) << thetadot          << ","
       << setw(15) << deltag                 << ","
       << setw(15) << da                << ","
       << 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);
}

////////////////////////////////////////
// Main program
//
int main() {
  cout << "boiler plate for lin. (OR non) Damped Simple Pendulum" << endl;
  cout << "w/ deadbeat type drive and g change" << endl;
  cout << "this version dev. of Verlet" << endl;

  //double period = 1;
    double period;
  double theta = 88e-3;
  //if (0) {
    //cout << "initial displacement angle in radians 5 deg. ~=88mrad:  ";
    //cin >> theta;
  //}
  double thetadot = 0;
  double t = 0;
    cout << "period (unit second)?";
  cin >> period;
    cout << "steps per cycle?  ";
  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;
    double da = 0;
  cin >> d;
  //cout << "g increment; hint: 0 or 1e-5:  ";
  cout << "quake amplitude";
    double deltag_max;
  cin >> deltag_max;
  cout << "desired number of cycles * number of steps per cycle = total steps" << endl;
    cout << "above steps per period" << endl;
    cout << steps_per_cycle << endl;
  cout << "number of steps total?:  " << endl;
  double stop = 0;
  cin >> stop;
  cout << "Printouts per cycle; > nine usually fittable:  ";
  double po_per_cycle = 0;
  cin >> po_per_cycle;
  double mod_number = steps_per_cycle / po_per_cycle;
    //double period;
  cout << "Data output filenname hint! string.txt:  ";
  string out_filename;
  cin >> out_filename;

  ofstream ouch;
  ouch.open(out_filename.c_str());
  if (!ouch.good()) {
    cerr << "Failed to open output file '" << out_filename << "'" << endl;
    exit(1);
  }

// column headers:
  ouch << setw(10) << Q+"t/s"+Q               << ","
       << setw(15) << Q+"theta/rad"+Q         << ","
       << setw(15) << Q+"thetadot"+Q          << ","
       << setw(15) << Q+"deltag_now"+Q        << ","
       << setw(15) << Q+"drive amp"+Q         << ","
       << endl;
// print the initial state:
  p(ouch, t, theta, thetadot, 0/*deltag*/, da);

///////////////////////////////////////////////////
// main loop
  for (int N = 1; N <= stop; N++){
    double da = 0;


      if ((theta < 0.08 && theta > 0.0 && thetadot <0.0)  || (theta > -0.08 && theta < 0.0 && thetadot >0.0) ){
          da = -d;}  else {da = 0;}

      if (thetadot >0.0) { da=-da; }  //drive

      double deltag_now = 0;
      if ((t > 1) && (t < 1.2 )) deltag_now = deltag_max;        // quake!

// acceleration at beginning of interval:
//  A = d^2 (theta) / dt^2  or [F/m]
    double A =  - r*thetadot +da
                - ((4*pi*pi)/(period*period))*(1+deltag_now) * theta;
    thetadot = A * deltat/2 + thetadot;
    theta = theta + deltat * thetadot;

// recalculate acceleration and end of interval, using new theta:
    A =   - r*thetadot +da
          - ((4*pi*pi)/(period*period))*(1+deltag_now)  * theta;
    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(ouch, t, theta, thetadot, deltag_now, da);
    }
  } // end main loop
}