#include <iostream.h>
#include <iomanip.h>
#include <math.h>

/*
Step 1. Document Code and Set Input Variables

This algorithm approximates the velocity and angular displacement of the bob of a pendulum using a time
stepping scheme to first update the velocity and then use the updated velocity to update the displacement
using Equations (1) and (2) of the previous section. 

INPUT 

     double g, gravitational acceleration (parameter) 
     double L, length of pendulum (input variable) 
     double dt, duration of time step (input variable) 
     integer IMAX, number of time steps (input variable) 
     double theta0, initial angular displacement (input variable) 
     double vel0, initial velocity (input variable) 

OUTPUT 

     integer i, time step counter for i = 0,1,...,IMAX 
     double t, time value at step i 
     double vel, approximate velocity at step i 
     double theta, approximate angular displacement at step i 
*/

const double g=9.8;  // set gravitational acceleration parameter



int main(){
	double L,dt,theta0,vel0,theta,vel,t;
	int IMAX,i;

	cout << "Enter length of pendulum (meters)" << endl;
    cin >> L;
	cout << "Enter initial angular displacement (radians)" << endl;
    cin >> theta0;
	cout << "Enter initial angular velocity (meters/second)" << endl;
    cin >> vel0;
	cout << "Enter duration of time step (seconds)" << endl;
    cin >> dt;
	cout << "Enter number of time steps)" << endl;
    cin >> IMAX;

// End of Step 1: Document and Set Input Variables and Parameters

// Step 2:  Initialize Numerical Solution 
    i = 0;
	t = 0;
	theta = theta0;
	vel = vel0;
	cout << "i= " << i << setprecision(3) << "  t=" << t << "\t"
		 << "theta=" << theta  << "\t" << "velocity=" << vel <<endl;

// End of Step 2

// Step 3:  Begin Time Stepping 
for(i=1;i<=IMAX;i++)
   {
// Step 4:  Advance Solution One Time Step 
	vel = vel - g*sin(theta)*dt;
    theta = theta + vel*dt/L;
    t = t + dt;

// Step 5. Output Numerical Solution  
	cout << "i= " << i << setprecision(3) << "  t=" << t << "\t"
		 << "theta=" << theta  << "\t" << "velocity=" << vel <<endl;	

	}
	return(1);
}