Geodesic4D.hpp

//
// $Id: Geodesic4D.hpp,v 1.2 2004/05/06 22:42:16 werner Exp $
//
// $Log: Geodesic4D.hpp,v $
// Revision 1.2  2004/05/06 22:42:16  werner
// Towards a specification of a spacetime via the Acceleration structure.
//
// Revision 1.1  2004/05/05 15:56:52  werner
// Separation of DOP core routines with dynamic size into the ODE library.
//
// Revision 1.2  2004/05/03 13:33:33  werner
// integration improved
//
// Revision 1.1  2004/04/28 14:06:37  werner
// Integrators may be safely derived from the dop853 template now.
// The template definition does not include the member functions, these
// need to be included explicitely when instantiating the dop853 integrator.
//
//
#ifndef __GEODESIC4D_HPP
#define __GEODESIC4D_HPP "Created 27.03.2004 13:48:45 by werner benger"

#include "Geodesic.hpp"
#include "Integrate4D.hpp"

namespace Traum
{
using namespace VecAl;

template <class Acceleration>
class   Geodesic4D : public IntegratePath4D
{
public: typedef IntegratePath4D::real real;

        typedef typename Acceleration::Scalar_t Scalar_t;
        typedef typename Acceleration::Point_t  Point_t;
        typedef typename Acceleration::Vector_t Vector_t;
        enum { dim = Acceleration::Dims };

        const Acceleration&myAcceleration;

        Geodesic4D(const Acceleration&A)
        : myAcceleration(A)
        {}

        void restart(const real&s, const Point_t&x0, const Vector_t&v0)
        {
                init(s, x0.const_ptr(), v0.const_ptr() );
        }

        void Accel(real, const real *x, const real *v, real *d2x_ds2)
        {
        Point_t  P;
        Vector_t V;
                for(int k=0; k<dim; k++)
                {
                        P[k] = x[k];
                        V[k] = v[k];
                }
        const Vector_t&a = myAcceleration(P, V);
                for(int i=0; i<nPaths(); i++)
                        d2x_ds2[i] = -a[i];
        }

        Point_t  position() const
        {
        Point_t  retval;
                for(int i=0; i<dim; i++)
                        retval[i] = q(i);
                return retval;
        }

        Vector_t velocity() const
        {
        Vector_t retval;
                for(int i=0; i<dim; i++)
                        retval[i] = dq(i);
                return retval;
        }

        Point_t  position(real s) const
        {
        Point_t  retval;
                for(int i=0; i<dim; i++)
                        retval[i] = q(i, s);
                return retval;
        }

        Vector_t velocity(real s) const
        {
        Vector_t retval;
                for(int i=0; i<dim; i++)
                        retval[i] = dq(i, s);
                return retval;
        }
};


} // namespace 

#endif // __GEODESIC4D_HPP