/**
 * pendulum.cpp
 *
 * The program assumes the form
 * 
 * theta_dot = omega
 * omega_dot = - (g/L) sin(theta)
 *
 * for the equation of motion of a pendulum.  We integrate
 * the equation using a forward Euler scheme.
 *
 * We estimate the period by doubling the time between
 * successive crossings of theta == 0.  To find a crossing,
 * we look for where theta and theta_old have opposite signs:
 * 
 * theta_k * theta_(k-1) < 0.
 * 
 * We then estimate the period using the associated times:
 * 
 * T_k = 2*(time_k - time_(k-1)).
 * 
 * The average period is
 * 
 * T_avg = (1/Ntot) (T_1 + T_2 + ... + T_Ntot),
 * 
 * with Ntot the total number of T evaluations.  The error bar
 * is given by error = s/sqrt(Ntot), with
 * 
 * s = sqrt( (1/(Ntot - 1) * ( (T_1 - Tavg)^2 + ... + (T_Ntot - T_avg)^2 ) ).
 *
 * We also check the energy per unit mass per unit length
 * squared at the initial and final times:
 * 
 * e/mL^2 = 0.5 theta_dot^2 - g/L * cos(theta) .
 * 
 */

#include <iostream>
#include <fstream>   // for writing to file
#include <cmath>     // for trig functions and sqrt

using namespace std;

int main(int argc, char** argv) {
  // ensure data has been fed through command line
  if(argc != 6) {
    cerr << "Usage: " << argv[0] << " n1 n2 n3 n4 n5" << endl;
    cerr << "where" << endl;
    cerr << "  n1: the initial angle in degrees" << endl;
    cerr << "  n2: the initial angular velocity" << endl;
    cerr << "  n3: the ratio of g and L" << endl;
    cerr << "  n4: the timestep" << endl;
    cerr << "  n5: the number of steps" << endl;
    return 1;
  }
  
  // declare parameters
  
  // initial angle
  // initial velocity
  // g/L
  // timestep
  // number of timesteps taken
  
  // time variable
  // time value for previous step
  // theta value for previous step
  // acceleration
  // number of reversals
  
  // read initial data
  th0    = atof(argv[1])*pi/180;       // convert to radians
  om0    = atof(argv[2]);
  gl     = atof(argv[3]);
  dt     = atof(argv[4]);
  nsteps = atoi(argv[5]);
  
  double *t_plot, *th_plot, *om_plot;  // plotting variables
  double *period;                      // for calculating avg. period
  t_plot  = new double [nsteps];
  th_plot = new double [nsteps];
  om_plot = new double [nsteps];
  period  = new double [nperiod];
  
  
  /*********************************
   *                               *
   * Integrate equations of motion *
   *                               *
   *********************************/
  // set initial state
  
  // principal integration loop
  for(int step = 0; step < nsteps; step++) {
    
    // increment time
    
    // calculate acceleration
    
    // Euler algorithm
    
    // test if pendulum has passed through theta == 0
    // if yes, estimate period
    if (theta*theta_old < 0) {
      // if first crossing
      // just record time
      
      // if not the first crossing
      // double the time b/w this crossing
      // and the last
      
      // increment counter for reversals
      
    }
  }
  
  /**********************************
   *                                *
   * Estimate period of oscillation *
   *                                *
   **********************************/
  // estimate average period, including error bar
  
  // take average of periods calculated above
  if(nperiod > 1) {
    
    // calculate std dev from above average
    
  }
  else {
  }
  
  /********************************
   *                              *
   * Estimate energy conservation *
   *                              *
   ********************************/
  
  // find initial and final energies, and relative error b/w them
  
  /**********************
   *                    *
   * Print data to file *
   *                    *
   **********************/
  ofstream data("pendulum.txt");
  for(int i = 0; i < nsteps; i++) {
    data << t_plot[i]  << "\t";
    data << th_plot[i] << "\t";
    data << om_plot[i] << endl;
  }
  // close file
  data.close();
  
  // release allocated memory
  delete [] t_plot;
  delete [] th_plot;
  delete [] om_plot;
  delete [] period;
}