//----------------------------------------------------------------------------
//	Copyright (C) 2002, 2004, 2007 Humboldt-Universitaet zu Berlin
//
//	This library is free software; you can redistribute it and/or
//	modify it under the terms of the GNU Lesser General Public
//	License as published by the Free Software Foundation; either
//	version 2.1 of the License, or (at your option) any later version.
//
//	This library is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//	Lesser General Public License for more details.
//
//	You should have received a copy of the GNU Lesser General Public
//	License along with this library; if not, write to the Free Software
//	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//----------------------------------------------------------------------------
/**	\file GSLSolver.cpp

	\author Michael Fiedler

	\date created at 2009/01/02

	\brief Implementation of class GSLSolver

	\sa GSLSolver.h

	\since 3.0
*/

#include <odemx/base/continuous/GSLSolver.h>
#include <odemx/base/continuous/GSLContainer.h>
#include <exception> //TODO: Implement exceptions.
extern "C" {
#include <gsl/gsl_errno.h>
}

#define PUT_IN_INTERVALL(a, b, c) c > b ? b : c < a ? a : c

using namespace odemx::base;
using namespace odemx::base::continuous;
using namespace std;

//------------------------------------------------------------------construction/destruction
GSLSolver::GSLSolver(int order)
:	nextStep(-1.0),
	step(0),
	control(0),
	evolve(0)
{
	switch (order) {
		case 0:
		case 1:
		case 2:
			T = gsl_odeiv_step_rk2; break;
		case 3:
			T = gsl_odeiv_step_rk4; break;
		case 4:
			T = gsl_odeiv_step_rkf45; break;
		case 5:
			T = gsl_odeiv_step_rkck; break;
		case 6:
		case 7:
		case 8:
			T = gsl_odeiv_step_rk8pd; break;
		case -1:
			T = gsl_odeiv_step_rk2imp; break;
		case -2:
			T = gsl_odeiv_step_rk4imp; break;
		case -3:
			T = gsl_odeiv_step_bsimp; break;
		case -4:
			T = gsl_odeiv_step_gear1; break;
		case -5:
			T = gsl_odeiv_step_gear2; break;
		default:
			T = gsl_odeiv_step_rkf45;
	}
	lastStep = -1.0;
	monitor = 0;
	//cout << "GSLSolverConstructor1" << endl;
}

GSLSolver::~GSLSolver() {
	if (evolve) {
		gsl_odeiv_evolve_free (evolve);
		evolve = 0;
	}
	if (control) {
		gsl_odeiv_control_free (control);
		control = 0;
	}
	if (step) {
		gsl_odeiv_step_free (step);
		step = 0;
	}
}

//-----------------------------------------------------------give suitable container
VariableContainer* GSLSolver::getVariableContainer(int size) {
	return new GSLContainer(size);
}

//--------------------------------------------------------------------redefinitions
void GSLSolver::evalF(double time) {
	//cout << "GSLSolver::evalF1" << endl;
	ODESolver::evalF(time);
}

void GSLSolver::evalJacobian(double time) {
	//cout << "GSLSolver::evalJacobian1" << endl;
	ODESolver::evalJacobian(time);
}

VariableContainer* GSLSolver::getVariableContainerOfMonitor() {
	return ODESolver::getVariableContainerOfMonitor();
}

//-----------------------------------give evaluation of monitor system in form for GSL
int gsl_func(double t, const double y[], double dydt[], void* params) {
	//cout << "GSLSolver::gsl_func1" << endl;
	GSLSolver* solv = (GSLSolver*)params;
	GSLContainer* container = (GSLContainer*)solv->getVariableContainerOfMonitor();
	container->setValues(y);
	//cout << "GSLSolver::gsl_func2" << endl;
	container->setDerivativesPointer(dydt);
	solv->evalF(t);
	//cout << "GSLSolver::gsl_func3" << endl;
	/*printf ("gsl_func: t = %.10f y = %+.10f %+.10f : dydt = %+.10f %+.10f \n",
			t, y[0], y[1], dydt[0], dydt[1]);
	 gsl_func ok in test
	 */
	return GSL_SUCCESS;
}

int gsl_jacobian(double t, const double y[], double * dfdy, double dfdt[], void* params) {
	cout << "# GSLSolver::gsl_jacobian1" << endl;
	GSLSolver* solv = (GSLSolver*)params;
	GSLContainer* container = (GSLContainer*)solv->getVariableContainerOfMonitor();
	container->setValues(y);
	container->setJacobiPointer(dfdy, dfdt);
	solv->evalJacobian(t);
	return GSL_SUCCESS;
}

//--------------------------------------------------------------------------calculate one step
void GSLSolver::makeStep(double time) {
	GSLContainer* container = (GSLContainer*)getVariableContainerOfMonitor();
	int size = container->size;

	//init system
	gsl_odeiv_system sys = {gsl_func, gsl_jacobian, size, this};

	double h, t = time;
	if (nextStep > 0)
		h = PUT_IN_INTERVALL(minStep, maxStep, nextStep);
	else
		h = (minStep + maxStep) / 2;

	int status;

		/*container->values first contains the state and after the function call contains
		the new state */

	status = gsl_odeiv_evolve_apply (evolve, control, step, &sys, &t, time + maxStep, &h, container->variables);
		/*status could be GSL_SUCCESS or an error code of the functions gsl_func or gsl_jacobian
		 errors in this functions must be handled via exception or here
		 error of calculation can be examined in e->yerr[] component by component
		*/
	if (status != GSL_SUCCESS) {
		//todo insert suitably exception
	}

	nextStep = h;

	lastStep = t - time;

	if (h < minStep) {
		/*throw exception*/
	}
}

//-------------------------------------------------------------------------init gsl structures
void GSLSolver::init() {
	GSLContainer* container = (GSLContainer*)getVariableContainerOfMonitor();
	int size = container->size;

	//init step
	if (evolve) {
		gsl_odeiv_evolve_free (evolve);
		evolve = 0;
	}
	if (control) {
		gsl_odeiv_control_free (control);
		control = 0;
	}
	if (step) {
		gsl_odeiv_step_free (step);
		step = 0;
	}
	step = gsl_odeiv_step_alloc (T, size);

	control = errorType == absolute ?
	gsl_odeiv_control_y_new(errorLimit, 0.0) : gsl_odeiv_control_y_new(0.0, errorLimit);

	evolve = gsl_odeiv_evolve_alloc (size);

	lastStep = -1.0;
	nextStep = -1.0;
}
