Geodesic853.hpp

//
// $Id: Geodesic853.hpp,v 1.7 2004/08/12 16:21:33 werner Exp $
//
// $Log: Geodesic853.hpp,v $
// Revision 1.7  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
// Revision 1.6  2004/05/13 12:06:37  werner
// Removed no-op code.
//
// Revision 1.5  2004/05/11 18:05:10  werner
// *** empty log message ***
//
// Revision 1.4  2004/05/06 22:42:16  werner
// Towards a specification of a spacetime via the Acceleration structure.
//
// Revision 1.3  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 __GEODESIC853_HPP
#define __GEODESIC853_HPP "Created 27.03.2004 13:48:45 by werner benger"

#include "Geodesic.hpp"
#include <ode/Integrate853.hpp>

namespace Traum
{
using namespace VecAl;

template <class Acceleration>
//class Geodesic853 : public IntegratePath853<long double>
class   Geodesic853 : public IntegratePath853<float>
{
        int components;
public:
        const Acceleration&myAcceleration;

        typedef typename Acceleration::Point_t  Point_t;
        typedef typename Acceleration::Vector_t Vector_t;

        enum { Dims = Acceleration::Dims };

        Geodesic853(const Acceleration&A)
        : myAcceleration(A)
        , components(0)
        {}

        void restart(int N, long double s, const Point_t*x0, const Vector_t*v0)
        {
                components = N*Dims;
                init(components, s, x0->const_ptr(), v0->const_ptr() );
        }

        void restart1(long double s, const Point_t&x0, const Vector_t&v0)
        {
                components = Dims;
                init(components, s, x0.const_ptr(), v0.const_ptr() );
        }

        void Accel(real, const real *x, const real *v, real *d2x_ds2)
        {
                for(int c=0; c<components/Dims; c++)
                {
                Point_t  P;
                Vector_t V;

                        for(int i=0; i<Dims; i++)
                        {
                                P[i] = x[c*Dims + i];
                                V[i] = v[c*Dims + i];
                        }

                const Vector_t&a = myAcceleration(P, V);

                        for(int j=0; j<Dims; j++)
                                d2x_ds2[j+c*Dims] = -a[j];
                }
        }

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

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

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

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

} // namespace Vecal

#endif // __GEODESIC853_HPP