/*
 *  rk.cpp
 *  
 Programme to integrate the 2D system in BIO 341.
 
 dx/dt = f1(x, y)
 
 dy/dt = f2(x, y)
 
 using a Runge-Kutta scheme. The two functions f1, f2 are coded below. Here we use the linear system:
 
 dx/dt = a x + b y
 
dy/dt = c x  + d y
 
 Note that to integrate a simple ODE like dy/dt = f(y), you use x as the time variable and increment it manually in the loop,
 and y is the function to be integrated. The starting point has to be chosen accurately.
 
 *   
 */

#include <cmath>
#include <iostream>
#include <fstream>

// Useful type declarations

typedef std::istream                zInStream;
typedef std::ostream                zOutStream;
typedef std::ifstream               zInFileStream;
typedef std::ofstream               zOutFileStream;
typedef std::fstream                zFileStream;
typedef std::istringstream          zInStringStream;
typedef std::ostringstream          zOutStringStream;
typedef std::ios                    zIos;

#define zString                      std::basic_string<char>

using namespace std;


// Forward declaration

void Runge_Kutta(const double a[4], double dt, double xn, double yn, double* pxn1, double* pyn1);


int main(int argc, char* argv[])
{
// ****************************************
    // Input and output file specs
	
	zString m_OutFileName;     // Name of file attached to the output stream

	bool m_IOSuccess;		// true = no error, false = error
	bool m_bOpenFlag;		// Flag showing whether the file is open or not

//    zInFileStream  m_inStream;
    zOutFileStream  m_outStream;

	// Check if the input data has been passed in via the command line: note that argv[0] = executable name; if not, we prompt for input
	
    long NSTEPS = 1;
    double x0 = 0.0;    // initial values
    double y0 = 0.0;
    
    double params[4];  // matrix elements
            
	if(argc == 1)
	{
	    std::cout << "****************************************" << endl;
	    std::cout << "Enter file name, N, 4 matrix elements, and initial point x0, y0" << endl;
	    std::cin >> m_OutFileName >> NSTEPS >> params[0] >> params[1] >> params[2] >> params[3] >> x0 >> y0;
	}
	else
	{
        std::cout << "Error, enter data at prompt" << endl;
        return 1;
    }

	m_outStream.open(m_OutFileName.c_str());

	// If the file could not be opened flag an error

	if(m_outStream.fail())
	{
	    std::cout << "ERROR - unable to open output file " << m_OutFileName << endl;
		m_IOSuccess = false;
		m_bOpenFlag = false;
		return 1;
	}
	else
	{
		m_IOSuccess = true;
		m_bOpenFlag = true;
	}
// ****************************************
			
    double dt = 0.01;
    
    double xn, yn;
    double xn1, yn1;
    
    xn = x0;
    yn = y0;
    
    for(long i=0; i<NSTEPS; ++i)
    {
        m_outStream << xn << " " << yn << endl;
        
        Runge_Kutta(params, dt, xn, yn, &xn1, &yn1);
        
        xn = xn1;
        yn = yn1;
    }
	

// ****************************************
// Clean up
	
	m_outStream.close();

    return 0;
}

// **********************************************************************
// The pair of functions to return the time derivatives at a given time.
//

double f1(const double a[4], double x, double y)
{
    return a[0]*x + a[1]*y;
}

double f2(const double a[4], double x, double y)
{
    return a[2]*x + a[3]*y;
}

// **********************************************************************
// Function to integrate the two functions f1(x(t), y(y)), f2(x(t), y(t)) using 4th order Runge-Kutta by one time-step.
// It receives the current point (xn, yn) and returns the new point (xn1, yn1) as pointers.
// Note that the functions are assumed not to be explicitly time-dependent.
//
// It uses the two subsidiary functions f1(x,y), f2(x,y) to get the derivatives.
//
// k_ij = coefficient with i = 1, 2, 3, 4, and j = x, y.

void Runge_Kutta(const double a[4], const double dt, const double xn, const double yn, double *pxn1, double *pyn1)
{
    const double k11 = f1(a, xn,yn)*dt;
    const double k12 = f2(a, xn,yn)*dt;

    const double k21 = f1(a, xn + 0.5*k11,yn + 0.5*k12)*dt;
    const double k22 = f2(a, xn + 0.5*k11,yn + 0.5*k12)*dt;

    const double k31 = f1(a, xn + 0.5*k21,yn + 0.5*k22)*dt;
    const double k32 = f2(a, xn + 0.5*k21,yn + 0.5*k22)*dt;
    
    const double k41 = f1(a, xn + k31,yn + k32)*dt;
    const double k42 = f2(a, xn + k31,yn + k32)*dt;

    *pxn1 = xn + 0.1666666*(k11 + 2.0*k21 + 2.0*k31 + k41);
    *pyn1 = yn + 0.1666666*(k12 + 2.0*k22 + 2.0*k32 + k42);
}