/* Pendulum program */

#include<math.h>
#define g 9.8 /* sets g for use throughout the program*/
#define deltat 0.001 /*sets a timestep value for use throughout the program*/
#include <iostream>
#include <fstream>
using namespace std;

float euler(float x,float xdot)
{
  /*This function just increments a variable by its derivative times */
  /*the time step.  This is Euler's method, which is a lot.*/
  /*like the rectangular rule for evaluating a simple integral. */  
  x=x+xdot*deltat;
  return(x);
}

void read_start_vals(float *theta, float *length)
{
  cout << "What is the initial value of theta in degrees?\n";
  cin >> *theta; /*theta is  a pointer so no & needed */
  /*The next line converts theta into radians */
  *theta=*theta*M_PI/180.0; /*M_PI is pi defined in math.h */
  /* An asterisk is needed for theta because we want the value*/
  /* at the memory address of the pointer, not the address itself*/
  cout << "What is the length of the pendulum in m?\n";
  cin >> *length;
  /*Like in the case where we read in theta, we already have a pointer */
  /*So no ampersand is needed.*/
  return;
}

float angular_acceleration(float theta,float length)
{
  float alpha,omega;
  alpha=-(g/length)*sin(theta);
  /* This is just first year physics, but in that context we make the small
     angle approximation. When solving numerically, we do not make that approximation.*/
  return alpha;
}


main()
{
  int i,imax;
  float alpha,theta,omega,length,t,kineticenergy,potentialenergy,totalenergy,velocity;
  ofstream outfile;
  /*Above lines define the variables*/
  /*Mass is ignored here, since it falls out of all equations*/
  /*To get energy in joules, etc. would have to add a few lines of code*/
  outfile.open("./pendulumoffset"); /* opens a file for writing the output*/
  t=0.0;
  omega=0.0;
  /*The two lines above make sure we start at zero.*/
  read_start_vals(&theta,&length); /*Reads initial offset and length of pendulum */
  for (i=0;i<50000; i=i+1) /*Starts a loop that will go through 50000 time steps */
  {
    alpha=angular_acceleration(theta,length);/*Computes angular accel. */
    omega=euler(omega,alpha);/*Changes the angular velocity by Euler meth. */
    theta=euler(theta,omega);/*Changes angle by Euler meth.*/
    potentialenergy=g*length*(1.0-cos(theta));/*Computes grav. pot en.*/
    velocity=omega*length;/* Computes linear velocity*/
    kineticenergy=0.5*velocity*velocity;/*Computes kin. en */
    totalenergy=potentialenergy+kineticenergy;/*Computes total energy*/
    t=t+deltat;/*Increments the time*/
    outfile << t <<" " << theta <<" " << omega <<" " << kineticenergy <<" " << potentialenergy <<" " << totalenergy << endl;/* Prints out the key variables*/
    /*Note: keeping track of the total energy is a good way to test a*/
    /*code like this where the energy should be conserved*/
  }  

}