// simple harmonic oscillator using Euler's method 
// P. Gorham, updated March 2013 for UH Physics 305

#include <iostream>
#include <iomanip>
#include <fstream>
#define USE_MATH_DEFINES
#include <cmath>
#include <cstdlib> 
using namespace std;

#define Tmax  50     // seconds for integration


main(int argc, char *argv[])
{
  double k,m,xt0,t0,vt0,w,vtii,xtii,t,dt,xti,vti, xtrue,vtrue;
   ofstream outfile;  
 
	outfile.open("osc1.dat");
	outfile.precision(5);         // set the precision for output

	k = 1.0;                   // spring constant
	m = 1.0;                   // mass in kg
	w = sqrt(k/m);
	xt0 = 1.0;                 // initial position
	t0 = 0.0;                  // initial time
	vt0 = 0.0;                 // initial velocity
	dt = 0.001;		   // time interval, 0.001=1ms is default

	xti = xt0;               // set & print out the initial conditions
	vti = vt0;
	t = t0;

/* Here is the loop that propagates the motion:
	vtii is the velocity at the new time, vti the previous step time;
	xtii is the new position, xti the previous;
	after each step is calculated, the old is set
	to the new, and the cycle is repeated  */

	for(t=t0; t<Tmax; t+= dt){
		vtii = vti + dt*(-k/m)*xti;
		xtii = xti + dt*vti;
		xti = xtii;  // set the old values to the new ones for next step
		vti = vtii;
		xtrue = xt0*cos(w*t);
		vtrue = vt0  + -w*sin(w*t);
		outfile << t << "  " << xtii << "  " << vtii << "  " << xtrue << "  " << vtrue << endl;
	   }


}