//Author: Alexander Sage
// Hw7
// Date: 11/30/12
//Purpose: to write a specific Ode integrator that allows a choice between 2nd order and 4th order Runge-Kutta methods. 
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

void probSize (int& neqns);

void probINTconds (double& t, double& tstop, double& dt, double y[]);

void probRHS (double t, double y[], double rhs[]);

int main()
{



//*


ofstream outStream;                          //Declare the output stream.
outStream.setf(ios::fixed);
outStream.setf(ios::showpoint);



//*



int i;                               //loop index
int neqns;                        //number of equations used
int choose;                // to give the option between second and fourth order 
double t, tstop, dt, thalf;              //initial conditions on which the program opperates
double *y;                            // the array y is used as space to represent x, y, vx, and vy
double *ytemp;                        //ytemp is used as slightly different values of y 
double *k1, *k2, *k3, *k4;      // to evaluate right hand values of ODE
double *rhs;

probSize(neqns);                  //retrieving the size of the arrays from the exterier program probSize
y = new double [neqns];                   //defining the size of all the arrays used
ytemp = new double [neqns];
k1 = new double [neqns];
k2 = new double [neqns];
k3 = new double [neqns];
k4 = new double [neqns];
rhs = new double [neqns];


cout << "problem size " << neqns << endl;              // checking to see if it worked correctly


//*


outStream.open ("results.dat");


//*


cout << "Would you like to use second order, or fourth order Runge-Kutta" << endl;  //letting user choose
cout << "type 2 for second order or type 4 for fourth order" << endl;              // between 2nd and 4th order
cin >> choose;

probINTconds(t, tstop, dt, y);      //take in initial values from expernal program

if (choose = 2)              //user chose second order Runge-Kutta
{
  while (t < tstop)
  {
    probRHS(t, y, rhs);      // evaluate values for array RHS
    for(i=0; i < neqns; i ++)  // using loops to calculate each individual element of the array
    k1[i] = dt*rhs[i];
    for(i=0; i < neqns; i ++)
    ytemp[i] = y[i] + 0.5*k1[i];
    thalf = t + 0.5*dt;
    probRHS(thalf, ytemp, rhs);

    for(i=0; i < neqns; i ++)
    y[i] = y[i] + dt*rhs[i];
    t = t + dt;


//*


    for(i=0; i < neqns; i ++)
    outStream << t << y[i] << endl;                //writing out initial and the y array 


//*

  }
}
if (choose = 4)           // user chose fourth oder Runge-Kutta
{
  while (t < tstop)
  {
    probRHS(t, y, rhs);      
    for(i=0; i < neqns; i ++)           //defining k2 with values of the RHS
    k1[i] = dt*rhs[i];
    for(i=0; i < neqns; i ++)
    ytemp[i] = y[i] + 0.5*k1[i];          //defining the array ytemp to be used in RHS
    thalf = t + 0.5*dt;
    probRHS(thalf, ytemp, rhs);

    for(i=0; i < neqns; i ++)
    k2[i] = dt*rhs[i];                    //defining k2 with values of the RHS
    for(i=0; i < neqns; i ++)
    ytemp[i] = y[i] + 0.5*k2[i];           //defining the array ytemp to be used in RHS
    thalf = t + 0.5*dt;
    probRHS(thalf, ytemp, rhs);

    for(i=0; i < neqns; i ++)               //defining k3 with values of the RHS
    k3[i] = dt*rhs[i];                    
    for(i=0; i < neqns; i ++)
    ytemp[i] = y[i] + 0.5*k3[i];            //defining the array ytemp to be used in RHS
    thalf = t + 0.5*dt;
    probRHS(thalf, ytemp, rhs);

    for(i=0; i < neqns; i ++)
    k4[i] = dt*rhs[i];                     //defining k4 with values of the RHS

    for(i=0; i < neqns; i ++)
    y[i] = y[i] + (k1[i] + 2.0*k2[i] + 2.0*k3[i] +k4[i])/6.0;
    t = t + dt;


*


    for(i=0; i < neqns; i ++)
    outStream << t << y[i] << endl;                //writing out initial t and the y array 


*


  }
}
else                                   // accounting for the user not inputting a 2 or a 4
{
  cout << choose << "Is an improper input" << endl;
}
return(0);
}