RungeKutta.hpp

//
// $Id: RungeKutta.hpp,v 1.1 2004/08/12 16:21:33 werner Exp $
//
// $Log: RungeKutta.hpp,v $
// Revision 1.1  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
// Revision 1.8  2004/05/12 14:36:22  werner
// Separation of chart definitions on a manifold and the tangential space.
// Introduced convenient coordinate transformations.
//
// Revision 1.7  2004/05/11 18:05:10  werner
//
// Revision 1.6  2004/05/11 16:53:47  werner
// Introduction of coordinate systems for improved type-safety.
// Will support easy coordinate transformations soon.
//
// Revision 1.5  2004/05/06 22:42:16  werner
// Towards a specification of a spacetime via the Acceleration structure.
//
// Revision 1.4  2004/05/05 15:56:52  werner
// Separation of DOP core routines with dynamic size into the ODE library.
//
// Revision 1.3  2004/05/03 13:33:33  werner
// integration improved
//
// Revision 1.2  2004/03/29 11:51:02  werner
// Common interface among simple integrators, DiffMe and Vecal, and preliminiary work on integrating dop853.
//
// Revision 1.1  2004/03/22 11:55:02  werner
// Schwarzschild geodesic integration.
//
// Revision 1.1  2004/02/13 16:36:21  werner
// Initial preliminiary version of the Vector Algebra Library.
//
#ifndef __RUNGE_KUTTA_Geodesic_HPP
#define __RUNGE_KUTTA_Geodesic_HPP "Created 20.07.2004 00:11:27 by werner"

#include <vecal/VVector.hpp>
#include <vecal/Matrix.hpp>

#include <ode/Integrator.hpp>
#include <spacetime/Geodesic.hpp>

namespace Traum
{

template <class Acceleration>
struct  RungeKuttaGeodesic : public IntegratorBase
{
        typedef typename Acceleration::Scalar_t     Scalar_t;
        typedef typename Acceleration::Point_t      Point_t;
        typedef typename Acceleration::Vector_t     Vector_t;
        typedef typename Acceleration::Christoffel_t Christoffel_t;

        typedef typename Point_t::Vector_t PointComponents_t;
        enum { Dims = PointComponents_t::SIZE };

        Point_t  x;
        Vector_t v;

        const Acceleration&Accel;
        Scalar_t ds;

        RungeKuttaGeodesic(const Acceleration&A)
        : Accel(A), ds(1)
        {}

        success_code advance(bool bk=false)
        {
        Scalar_t ds = bk?-this->ds:this->ds;

        Point_t  P = x + v*(ds/2.L);
        Vector_t V = v + v/2.L;

        Vector_t
                b1 = ds * Accel(x, v),
                b2 = ds * Accel(x + v*(ds/2.L) + b1*(ds/8.L), v + b1/2.L),
                b3 = ds * Accel(x + v*(ds/2.L) + b1*(ds/8.L), v + b2/2.L),
                b4 = ds * Accel(x + v*ds    + b3*(ds/2.L), v + b3);

                x += v*ds + (b1 +    b2 +    b3) * ds/6;
                v +=        (b1 +2.L*b2 +2.L*b3 + b4)/6;

                return StepOk;
        }
};


} /* namespace VecAl */ 

#endif /* __RUNGE_KUTTA_Geodesic_HPP */