AdamsStoermer.hpp

//
// $Id: AdamsStoermer.hpp,v 1.1 2004/08/12 16:21:33 werner Exp $
//
// $Log: AdamsStoermer.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 __AdamsStoermer_Geodesic_HPP
#define __AdamsStoermer_Geodesic_HPP "Created 27.02.2001 21:42:27 by werner"

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

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

namespace Traum
{

template <class Acceleration>
class   AdamsStoermerGeodesic : public RungeKuttaGeodesic<Acceleration>
{
public:
        typedef typename RungeKuttaGeodesic<Acceleration>::Scalar_t Scalar_t;
        typedef typename RungeKuttaGeodesic<Acceleration>::Point_t  Point_t;
        typedef typename RungeKuttaGeodesic<Acceleration>::Vector_t Vector_t;

        using RungeKuttaGeodesic<Acceleration>::x;
        using RungeKuttaGeodesic<Acceleration>::v;

private:
        int step;
        Vector_t b0, b1, b2, b3;
        
        void predict(Point_t&xp, Vector_t&vp, const Scalar_t&ds) const
        {
                xp = x + v*ds + ds*ds/360.*(323.*b0-264.*b1+159.*b2-38.*b3);
                vp = v        +    ds/24. *( 55.*b0- 59.*b1+ 37.*b2- 9.*b3);
        }

        void correct(Point_t&xp, Vector_t&vp, const Scalar_t&ds) const
        {
        Vector_t bp = -Accel(xp, vp);
                xp = x + v*ds
                       + ds*ds/1440.*(135.*bp+752.*b0-246.*b1+ 96.*b2-17.*b3);
                vp = v +    ds/720. *(251.*bp+646.*b0-264.*b1+106.*b2-19.*b3);
        }

        void shift()
        {
                b3 = b2;
                b2 = b1;
                b1 = b0;
                b0 = -Accel(x, v);
        }

public:
        AdamsStoermerGeodesic (const Acceleration&A)
        : RungeKuttaGeodesic<Acceleration>(A)
        , step(0)
        {}

        void    reset_step()
        {
                step = 0;
        }

        success_code advance(bool bk=false)
        {
        Scalar_t ds = bk?-this->ds:this->ds;
        
                if (step++ < 4) 
                {
                        RungeKuttaGeodesic<Acceleration>::advance(bk);
                        shift();
                        return StepOk;
                }
        Point_t xp;
        Vector_t vp;
                predict(xp, vp, ds);
                correct(xp, vp, ds);
                correct(xp, vp, ds);
                shift();
                x = xp;
                v = vp;

                return StepOk;
        }
};

} /* namespace VecAl */ 

#endif /* __AdamsStoermer_Geodesic_HPP */