/*------------------------------------------------------*/
/* main.cpp												*/
/* Author: Aaron Pulley									*/
/* (804) 714-9052 - awpulley@cavtel.net					*/
/* Created: 2 Apr 2010									*/
/* RK-4 algorithm.										*/
/* History of modifications:							*/
/* - 30 Mar 2011 Modified from lab9main.cpp for			*/
/*   ME EN 373-3.										*/
/* -													*/
/*------------------------------------------------------*/
#include <iostream>
#include <fstream>
#include <cmath>
#include <string>
using namespace std;

//function prototypes
double xdot(double v);
double vdot(double x, double v);

int main()
{
	//initialize
	const double time(10), t0(0), x0(0.3), v0(0); //calculate over 10s, initial values
	double h, n, t, x, k1x, k2x, k3x, k4x, v, k1v, k2v, k3v, k4v;
	t=t0; //set initial values
	x=x0;
	v=v0;
	string filename;
	
	//Tells computer to expect file output
	ofstream output;

	//intro
	cout<<"This program performs the RK-4 method for lab 10."<<endl;

	//create output file
	cout<<"Please input a name for the output file (including extension):"<<endl;
	cin>>filename;
	cout<<endl;
	output.open(filename.c_str());

	//get step size
	cout<<"Please input step size (in seconds): ";
	cin>>h;
	cout<<endl;
	n=time/h;

	//Output/write initial row of data separated by tabs
	cout << t << "\t" << x << "\t" << v << endl;
	output << t << "\t" << x << "\t" << v << endl;

	//RK-4 loop
	for(int i=1; i<=n; i++)
	{
		//k values
		k1x=xdot(v);
		k1v=vdot(x, v);

		k2x=xdot(v+k1v*h/2);
		k2v=vdot(x+k1x*h/2, v+k1v*h/2);

		k3x=xdot(v+k2v*h/2);
		k3v=vdot(x+k2x*h/2, v+k2v*h/2);

		k4x=xdot(v+k3v*h);
		k4v=vdot(x+k3x*h, v+k3v*h);

		//next t, x, v
		t=t0+i*h;
		x=x+h*(k1x+2*k2x+2*k3x+k4x)/6;
		v=v+h*(k1v+2*k2v+2*k3v+k4v)/6;

		//Output/write row of data separated by tabs
		cout << t << "\t" << x << "\t" << v << endl;
		output << t << "\t" << x << "\t" << v << endl;
	}

	//conclude
	cout<<"End of computation."<<endl;
	system("pause");

	return 0;
}

//functions
double xdot(double v)
{
	return v;
}

double vdot(double x, double v)
{
	return -0.8*v-x;
}