// funcyion to calculate step through ODE by
// 4th and 5th order Runge Kutta Fehlberg.
// calculates both orders and gives difference back
// ak 2/2/2001
// all in double precision
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <iomanip.h>
#include <stdio.h>
void deri( double x , double y[] , double f[], double omega2 , double gam);
double omega2;
double gam;
int len=2;
	double y[2],y_old[2];
	double f0[2] , f1[2], f2[2] , f3[2] , f4[2], f[2];

main()
 {
 	// constants for the RKF method
	double const a0  =.25 , a1=3./8. , a2=12./13. ;
	double const b10  =.25 ;
	double const b20 =3./32. , b21 = 9./32. ;
	double const b30 =1932./2197. , b31 = -7200./2197. , b32 = 7296./2197. ;
	double const b40 =439./216. , b41 = -8. , b42=3680./513. ,b43=-845./4104.;
	double const r1 = 25./216. , r2 = 1408./2565. , r3 = 2197./4104. , r4=
	-1./5. ;
	// end of constants

	// user input
	double spring=0.; // spring constant
	double mass=0.;   // mass
	double damp=0.;	  // damping
	double v_0=0.;	  // initial condition for velcoity
	double x_0=0.;	  // initial condition for amplitude
	double step=0.;   // stepsize
	char filename[80]; // plotfile
	// end user input
	
	double y[2],y_old[2];
	double e_old, e_new;
	double f0[2] , f1[2], f2[2] , f3[2] , f4[2], f[2];
	double result=0.;
	double x=0 , x0;
	int len=2;
	

// *******************************************************************
//	here we deal with the user input
	
	cout << " give the filename for output : ";
	cin >> filename;
 
 // Now open the file
	
	ofstream out_file(filename) ; // Open file
		out_file.setf(ios::showpoint); // keep decimal point
	
	cout << " give the step size: ";
	cin >> step ;
	cout << "mass and damping ";
	cin >> mass >> damp;
	cout << " spring constant ";
	cin >> spring;
	cout << " initial position and velocity ";
	cin >> x_0 >> v_0 ;
	omega2=spring/mass;
	gam=damp/(mass*2.);
	
//********************************************************************
//	input is finished
// 	let's do the calculation
	
//	use the initial conditions:
	x=0.;
	 y[0]=x_0;
	 y[1]=v_0;
	 y_old[0]=y[0];
	 y_old[1]=y[1];

	while(x<20.)
	{
//	here starts the loop
	x0=x;	
	deri(x0,&y[0],&f[0],omega2,gam);
	  x=x0+a0*step;
	   y[0]=y_old[0]+b10*step*f[0];
	   y[1]=y_old[1]+b10*step*f[1];
	deri(x,&y[0],&f1[0],omega2,gam);
	  x=x0+a1*step;
	   y[0]=y_old[0]+b20*step*f[0]+b21*step*f1[0];
	   y[1]=y_old[1]+b20*step*f[1]+b21*step*f1[1];
	deri(x,&y[0],&f2[0],omega2,gam);
	  x=x0+a2*step;
	   y[0]=y_old[0]+b30*step*f[0]+b31*step*f1[0]+b32*step*f2[0];
	   y[1]=y_old[1]+b30*step*f[1]+b31*step*f1[1]+b32*step*f2[1];
	deri(x,&y[0],&f3[0],omega2,gam);
	  x=x0+step;
	   y[0]=y_old[0]+b40*step*f[0]+b41*step*f1[0]+b42*step*f2[0]+b43*step*f3[0];
	   y[1]=y_old[1]+b40*step*f[1]+b41*step*f1[1]+b42*step*f2[1]+b43*step*f3[1];
	deri(x,&y[0],&f4[0],omega2,gam);
	y[0]=y_old[0]+step*(r1*f[0]+r2*f2[0]+r3*f3[0]+r4*f4[0]);
	y[1]=y_old[1]+step*(r1*f[1]+r2*f2[1]+r3*f3[1]+r4*f4[1]);
	
//	cout << x << " " << y[0]<<" " << y[1] << "\n";
	out_file << x << "   "<<y[1] <<"\n";
	y_old[0]=y[0];
	y_old[1]=y[1];
	}
//	cout << f[0] << "   "<<f[1] <<"\n";
//	cout << f1[0] << "   "<<f1[1] <<"\n";
	out_file.close();
	
	}
	
void deri( double x , double y[ ] , double f[ ], double omega2,double gam)
{
	f[0]=y[1];
	f[1]=-omega2*y[0]-2*gam*y[1];
	return;
}
// compile g++ damp.cxx