IntegrateGeodesic.hpp

//
// $Id: IntegrateGeodesic.hpp,v 1.1 2004/08/12 16:21:33 werner Exp $
//
// $Log: IntegrateGeodesic.hpp,v $
// Revision 1.1  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
//
#ifndef __INTEGRATE_GEODESIC_HPP
#define __INTEGRATE_GEODESIC_HPP "Created 20.07.2004 00:19:45 by werner benger"

#include <spacetime/Geodesic853.hpp>
#include <spacetime/EulerGeodesic.hpp>
#include <spacetime/RungeKutta.hpp>
#include <spacetime/AdamsStoermer.hpp>

namespace Traum
{
using namespace VecAl;

class   GeodesicIntegratorBase
{
public:

        enum IntegratorType
        {
                Euler,
                RungeKutta,
                AdamsStoermer,
                DOP853
        };
        const IntegratorType myType;

        GeodesicIntegratorBase(IntegratorType t)
        : myType(t)
        {}
};

template <class TangentialSpaceType>
class   GeodesicIntegrator : public GeodesicIntegratorBase
{
public:
        typedef typename TangentialSpaceType::Point_t  Point_t;
        typedef typename TangentialSpaceType::Vector_t Vector_t;

        enum { Dims = TangentialSpaceType::Dims };

        GeodesicIntegrator(IntegratorType t)
        : GeodesicIntegratorBase(t)
        {}

        virtual success_code advance(bool bk=false) = 0;
        virtual void restart(const Point_t&x0, const Vector_t&v0) = 0;
        virtual Point_t  position() const = 0;
        virtual Vector_t velocity() const = 0;
};

template <class Acceleration>
class IntegrateGeodesic : public GeodesicIntegrator<typename Acceleration::TangentialSpace_t>
{
public:
        typedef typename Acceleration::TangentialSpace_t TangentialSpace_t;
        typedef GeodesicIntegrator<TangentialSpace_t>    Base_t;

        typedef typename Base_t::Point_t  Point_t;
        typedef typename Base_t::Vector_t Vector_t;

        typedef typename Acceleration::Point_t  GenericPoint_t;
        typedef typename Acceleration::Vector_t GenericVector_t;

        enum { Dims = Acceleration::Dims };

        EulerGeodesic<Acceleration>             E;
        RungeKuttaGeodesic<Acceleration>        RK;
        AdamsStoermerGeodesic<Acceleration>     AS;
        Geodesic853<Acceleration>               DP853;

        IntegrateGeodesic(const Acceleration&A, IntegratorType t)
        : Base_t(t)
        , E(A)
        , RK(A)
        , AS(A)
        , DP853(A)
        {}

        override success_code advance(bool bk=false)
        {
                switch(myType)
                {
                case Euler:         return E.advance(bk);
                case RungeKutta:    return RK.advance(bk);
                case AdamsStoermer: return AS.advance(bk);
                case DOP853:        return DP853.advance(bk);
                }
        }

        void setStepSize(double ds)
        {
                E.ds  = ds;
                RK.ds = ds;
                AS.ds = ds;
        }

        override void restart(const Point_t&x0, const Vector_t&v0)      
        {
                switch(myType)
                {
                case Euler:         E    .x = x0; 
                                    E    .v = v0; 
                                    break;

                case RungeKutta:    RK   .x = x0; 
                                    RK   .v = v0; 
                                    break;

                case AdamsStoermer: AS.reset_step();
                                    AS   .x = x0; 
                                    AS   .v = v0; 
                                    break;

                case DOP853:        DP853.restart1(0, x0, v0); 
                                    break;
                }
        }


        override Point_t  position() const
        {
                switch(myType)
                {
                case Euler:         return E.x;
                case RungeKutta:    return RK.x;
                case AdamsStoermer: return AS.x;
                case DOP853:        return DP853.position();
                }
        }

        override Vector_t velocity() const
        {
                switch(myType)
                {
                case Euler:         return E.v;
                case RungeKutta:    return RK.v;
                case AdamsStoermer: return AS.v;
                case DOP853:        return DP853.velocity();
                }
        }
};

} // namespace Vecal

#endif // __INTEGRATE_GEODESIC_HPP