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)
        {}
 
        success_code advance(bool bk=false) override
        {
                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;
        }
 
        void restart(const Point_t&x0, const Vector_t&v0) override
        {
                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;
                }
        }
 
 
        Point_t  position() const override
        {
                switch(myType)
                {
                case Euler:         return E.x;
                case RungeKutta:    return RK.x;
                case AdamsStoermer: return AS.x;
                case DOP853:        return DP853.position();
                }
        }
 
        Vector_t velocity() const override
        {
                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