#ifndef _utl_RK4ODEIntegrator_h_
#define _utl_RK4ODEIntegrator_h_

namespace utl {

  /**
    \class RK4Iterator RK4ODEIntegrator.h "utl/RK4ODEIntegrator.h"

    This iterator class will advance the solution for the current point.
    See class RK4ODEIntegrator for details.

    \author Darko Veberic
    \date 20 Jun 2008
    \version $Id: RK4ODEIntegrator.h 14717 2009-09-17 20:24:36Z lukas $

    \ingroup math
  */
  template<class DerivativeFunctor, class Vector>
  class RK4Iterator {
  public:
    RK4Iterator(const double x, const Vector& y,
                DerivativeFunctor& d)
      : fX(x), fDerivativeFunctor(d)
    { Copy(d, y, fY); }

    double GetX() const { return fX; }

    Vector& GetY() { return fY; }

    const Vector& GetY() const { return fY; }

    void
    Advance(const double dx, const Vector& dYdX)
    {
      const unsigned int n = fDerivativeFunctor;
      // first step: x -> x + dx/2
      const double dx2 = 0.5 * dx;
      Vector y1;
      for (unsigned int i = 0; i < n; ++i)
        y1[i] = fY[i] + dx2 * dYdX[i];
      // second step
      const double x1 = fX + dx2;
      Vector dYdX1;
      fDerivativeFunctor(x1, y1, dYdX1);
      for (unsigned int i = 0; i < n; ++i)
        y1[i] = fY[i] + dx2 * dYdX1[i];
      // third step
      Vector dYdX2;
      fDerivativeFunctor(x1, y1, dYdX2);
      for (unsigned int i = 0; i < n; ++i) {
        y1[i] = fY[i] + dx * dYdX2[i];
        dYdX2[i] += dYdX1[i];
      }
      // fourth step
      fX += dx;
      fDerivativeFunctor(fX, y1, dYdX1);
      const double dx6 = (1./6) * dx;
      for (unsigned int i = 0; i < n; ++i)
        fY[i] += dx6 * (dYdX[i] + dYdX1[i] + 2.*dYdX2[i]);
    }

    void
    Advance(const double dx)
    {
      Vector dYdX;
      fDerivativeFunctor(fX, fY, dYdX);
      Advance(dx, dYdX);
    }

    RK4Iterator&
    operator+=(const double dx)
    { Advance(dx); return *this; }

    RK4Iterator&
    operator-=(const double dx)
    { Advance(-dx); return *this; }

  private:
    static
    void
    Copy(const unsigned int n, const Vector& a, Vector& b)
    {
      for (unsigned int i = 0; i < n; ++i)
        b[i] = a[i];
    }

    double fX;
    Vector fY;
    DerivativeFunctor& fDerivativeFunctor;
  };


  /**
    \class RK4ODEIntegrator RK4ODEIntegrator.h "utl/RK4ODEIntegrator.h"

    This class prepares the integration of the Ordinary Differential Equations
    (ODE) of the form
    \f[
      \frac{{\rm d}y_i}{{\rm d}x} = f_i(x, y_0, \ldots, y_{n-1})
    \f]
    for \f$i\in{0, \ldots, n-1}\f$.

    ODE is defined through an function class with operator() defined
    (usualy such classes are called functor in C++). This functor has to
    define the operator() and unsigned int cast operator:
    \code
      class MyODE {
      public:
        void
        operator()(const double x, const double y[10], double dydx[10])
        {
          dydx[0] = ...
          dydx[1] = ...
        }

        operator unsigned int() { return 10; }
      };
    \endcode
    Operator() calculates derivatives \f${\rm d}y_i/{\rm d}x\f$ of variables
    \f$y_i\f$ with respect to the independent variable x.

    The cast operator unsigned int() returns number of OD equations (in this
    example 10).

    The class RK4ODEIntegrator returns an iterator that can be advanced with
    operator += dx for some step dx in independent variable.
    \code
      typedef double Vector[10];
      RK4ODEIntegrator<MyODE, Vector> ode(MyODE());

      RK4Iterator<MyODE, Vector> it = ode.Begin(x0, y0);

      it += dx;
    \endcode

    Note that the suitable candidates for template parameter Vector
    are vector classes with fixed size (std::vector<> is not appropriate
    since the default constructor leaves the object with zero size) that
    support the operator[] for element access. You can use the C-style
    arrays double[n] or Offline class utl::SVector<n>.

    See also utl::RK5ODEIntegrator and testODEIntegrator.cc for a real
    examples.

    \author Darko Veberic
    \date 20 Jun 2008
    \version $Id: RK4ODEIntegrator.h 14717 2009-09-17 20:24:36Z lukas $

    \ingroup math
  */
  template<class DerivativeFunctor>
  class RK4ODEIntegrator {
  public:
    RK4ODEIntegrator(DerivativeFunctor& d)
      : fDerivativeFunctor(d) { }

    template<class Vector>
    RK4Iterator<DerivativeFunctor, Vector>
    Begin(const double x, const Vector& y)
    {
      return
        RK4Iterator<DerivativeFunctor, Vector>(x, y, fDerivativeFunctor);
    }

    DerivativeFunctor& GetDerivativeFunctor()
    { return fDerivativeFunctor; }

  private:
    DerivativeFunctor& fDerivativeFunctor;
  };

}

#endif