//! \addtogroup Math
//! \{
//! \file pimathsolver.h
//! \brief
//! \~english Mathematical solver for differential equations
//! \~russian Математический решатель дифференциальных уравнений
//! \details
//! \~english Solver for ordinary differential equations using various numerical methods
//! \~russian Решатель обыкновенных дифференциальных уравнений с использованием различных численных методов
//! \}
/*
    PIP - Platform Independent Primitives
    PIMathSolver
    Ivan Pelipenko peri4ko@yandex.ru

    This program 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 3 of the License, or
    (at your option) any later version.

    This program 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 program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef PIMATHSOLVER_H
#define PIMATHSOLVER_H

#include "pimathmatrix.h"

/// Differential evaluations

//! Transfer function representation
//! \~english Structure representing transfer function with numerator and denominator coefficients
//! \~russian Структура, представляющая передаточную функцию с коэффициентами числителя и знаменателя
struct PIP_EXPORT TransferFunction {
	PIVector<double> vector_Bm, vector_An;
};

//! Mathematical solver for differential equations
//! \~\english Solver for ordinary differential equations using various numerical methods
//! \~russian Решатель обыкновенных дифференциальных уравнений
class PIP_EXPORT PIMathSolver {
public:
	//! Solving methods for differential equations
	enum Method {
		Global                    = -1, //!< Use global method
		Eyler_1                   = 01, //!< Euler method (first order)
		Eyler_2                   = 02, //!< Euler method (second order)
		EylerKoshi                = 03, //!< Euler-Cauchy method
		RungeKutta_4              = 14, //!< Runge-Kutta 4th order
		AdamsBashfortMoulton_2    = 22, //!< Adams-Bashforth-Moulton 2nd order
		AdamsBashfortMoulton_3    = 23, //!< Adams-Bashforth-Moulton 3rd order
		AdamsBashfortMoulton_4    = 24, //!< Adams-Bashforth-Moulton 4th order
		PolynomialApproximation_2 = 32, //!< Polynomial approximation 2nd order
		PolynomialApproximation_3 = 33, //!< Polynomial approximation 3rd order
		PolynomialApproximation_4 = 34, //!< Polynomial approximation 4th order
		PolynomialApproximation_5 = 35  //!< Polynomial approximation 5th order
	};

	//! Constructs an empty solver
	PIMathSolver();

	//! Solve differential equation with step h
	//! \~english Solve differential equation at point u with step h
	//! \~russian Решить дифференциальное уравнение в точке u с шагом h
	void solve(double u, double h);

	//! Initialize from transfer function
	//! \~english Set up solver from transfer function coefficients
	//! \~russian Инициализировать решатель из коэффициентов передаточной функции
	void fromTF(const TransferFunction & TF);
	//! Set solving method
	//! \~english Set numerical method for solving
	//! \~russian Установить численный метод решения
	void setMethod(Method m) { method = m; }
	//! Set simulation time
	//! \~english Set simulation time
	//! \~russian Установить время моделирования
	void setTime(double time);

	//! Solve using Euler method (1st order)
	void solveEyler1(double u, double h);
	//! Solve using Euler method (2nd order)
	void solveEyler2(double u, double h);
	//! Solve using Runge-Kutta 4th order
	void solveRK4(double u, double h);
	//! Solve using Adams-Bashforth-Moulton 2nd order
	void solveABM2(double u, double h);
	//! Solve using Adams-Bashforth-Moulton 3rd order
	void solveABM3(double u, double h);
	//! Solve using Adams-Bashforth-Moulton 4th order
	void solveABM4(double u, double h);
	//! Solve using polynomial approximation
	void solvePA(double u, double h, uint deg);
	//! Solve using polynomial approximation 2nd order
	void solvePA2(double u, double h);
	//! Solve using polynomial approximation 3rd order
	void solvePA3(double u, double h);
	//! Solve using polynomial approximation 4th order
	void solvePA4(double u, double h);
	//! Solve using polynomial approximation 5th order
	void solvePA5(double u, double h);

	//! Solution vector
	PIMathVectord X;
	//! Global default method
	static Method method_global;
	//! Description of available methods
	static const char methods_desc[];

private:
	void moveF();

	PIMathMatrixd A, M;
	PIMathVectord d, a1, b1;
	PIMathVectord k1, k2, k3, k4, xx;
	PIMathVectord XX, Y, pY;
	PIVector<PIMathVectord> F;
	PIVector<double> times;
	uint size = 0, step = 0;
	Method method = Global;
	double sum = 0., td = 0., ct = 0., lp = 0., dh = 0., t = 0., x1 = 0., x0 = 0.;
	bool ok = false;
};

#endif // PIMATHSOLVER_H