//----------------------------------------------------------------------------
//	Copyright (C) 2002, 2003, 2004 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.h

 \author Michael Fiedler

 \date created at 2008/11/21

 \brief Declaration of GSLSolver

 \since 3.0
 */

#ifndef ODEMX_GSLSOLVER_INCLUDED
#define ODEMX_GSLSOLVER_INCLUDED

#include <odemx/base/continuous/ODESolver.h>
#include <gsl/gsl_odeiv.h>

/** give function in GSL form global function
 */
int gsl_func(double t, const double y[], double dydt[], void* params);

/** give jacobi in GSL form global function
 */
int gsl_jacobian(double t, const double y[], double * dfdy, double dfdt[], void* params);

namespace odemx {
	namespace base {
		namespace continuous {

			//forward declaration
			class Monitor;

			class Continuous;

			class VariableContainer;

			class ODEObject;

			class ODESolver;

			/** \class GSLSolver

				\ingroup base

				\brief An Implementation of ODESolver which uses the GNU scientific library

				\sa ODESolver

				\since 3.0

			 */
			class GSLSolver : public ODESolver {
			public:
				/** \brief Construction
					\param order
						This parameter selects the used numeric solving method.

					explicit:

						0..2 Runge-Kutta 2

						3 Runge-Kutta 4

						4 Runge-Kutta-Fehlberg 4 5

						5 Runge-Kutta Cash-Karp 4 5

						6..8 Runge-Kutta Prince-Dormand 8 9

					implicit:

						-1 Runge-Kutta 2

						-2 Runge-Kutta 4

						-3 Bulirsch-Stoer

						-4 Gear method M=1

						-5 Gear method M=2

					\note For positive values of order the number to select a method is nearly the consistency order.
					This means the error in one step is proportional to (step size)^(consistency order)
					\note Explicit methods are faster, but have problems with stiff equations.
						Implicit methods can deal with stiff equations.
				 */
				GSLSolver(int order);

				/// destruction
				~GSLSolver();

				/** \brief gives an VariablleContainer suitable for this Solver
					\param size
						initial size of VariableContainer
				 */
				VariableContainer* getVariableContainer(int size);

				/** \brief makes one integration step for monitor
					\param time
						internal time of assigned monitor plus \p time is the time from which to do the step
				 */
				void makeStep(double time);

				/** \brief initialize the  solver

					\note this is also called if monitor has changed
				 */
				void init();

			protected:
				/* the following functions are redefined with implementation linked to
				 base class for friend declaration

				 @{
				 */
				/** \brief gives derivatives of the system described by the assigned monitor

				 evaluates all equations of all continuous process as one function
				 as sum of partial functions at time internalTime + time
				 */
				void evalF(double time);

				/** \brief gives jacbian of the system described in assigned monitor

				 evaluates all jacobian of all continuous process as one matrix
				 as sum of partial matrices at time internalTime + time
				 */
				void evalJacobian(double time);
				//@}

				/** \brief gives the variables of the assigned Monitor
				 */
				VariableContainer* getVariableContainerOfMonitor();

			private:
				double nextStep; //for evolving function

				/** solving algorithm
				 */
				const gsl_odeiv_step_type * T;

				/** step object
				 */
				gsl_odeiv_step * step;

				/** error control object for stepsize adaption
				 */
				gsl_odeiv_control * control;

				/** evolve object for stepwise evolution of solution
				 */
				gsl_odeiv_evolve * evolve;

				//making this functions global is done, because of type mismatch otherwise
				/** give function in GSL form global function
				 */
				friend int ::gsl_func(double t, const double y[], double dydt[], void* params);

				/** give jacobi in GSL form global function
				 */
				friend int ::gsl_jacobian(double t, const double y[], double * dfdy, double dfdt[], void* params);
			};
		}
	}
}

#endif

