/*
 *  rk.cpp
 *  
 Programme to integrate a pair of (non-linear) dynamical equations
 
 dx/dt = f1(x,y)
 
 dy/dt = f2(x,y)
 
 using a Runge-Kutta scheme. The two functions f1, f2 are coded below.
 
 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
	
    double amat[4];
    long NSTEPS = 100;
    double x0 = 0.0;
    double y0 = 0.0;
    
    // For 1d, we enter only the diagonal elements of the matrix as the equations are uncoupled. The x equation is just time
    // incrementing, and the y equation is the function to be integrated.
    
    double a, b;
    
    amat[0] = 0.0;
    amat[1] = 0.0;
    amat[2] = 0.0;
    amat[3] = 0.0;

//	if(argc == 1)
	{
	    std::cout << "****************************************" << endl;
	    std::cout << "Enter file name, N, a[4], and initial point x0, y0" << endl;
	    std::cin >> m_OutFileName >> NSTEPS >> amat[0] >> amat[1] >> amat[2] >> amat[3] >> x0 >> y0;
	}
/*	else
	{
        std::cout << "Error, only enter file name" << zEndl;
        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
	{
	    std::cout << "Write trajectory to data file: " << m_OutFileName  << endl;
		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(amat, 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 the given point.
//
// Following is the deformed larva portrait
// x equation   return -y + x*(1.0 - x*x - y*y) + (x-1.0)*(2.0 - x*x - y*y);
// y equation   return  x + y*(1.0 - x*x - y*y) + (y-1.0)*(2.0 - x*x - y*y);

// Following produces a ball and star
//return -y + x*(1.0 - x*x - y*y) + (x-4.0)*(1.0 - x*x - y*y);
//return x + y*(1.0 - x*x - y*y) + (y-4.0)*(1.0 - x*x - y*y);

// Following is the leaky limit cycle
//return -y + x*(1.0 - x*x - y*y) + (x-1.0)*(4.0 - x*x - y*y);
//return x + y*(1.0 - x*x - y*y) + (y-1.0)*(4.0 - x*x - y*y);


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

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

// **********************************************************************
// Function to integrate the two functions x(t),y(t) using 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);
}