Integrate4D.hpp

//
// $Id: Integrate4D.hpp,v 1.6 2007/03/28 22:31:38 georg Exp $
//
// $Log: Integrate4D.hpp,v $
// Revision 1.6  2007/03/28 22:31:38  georg
// dtors adapted for gcc4
//
// Revision 1.5  2007/02/22 17:08:42  werner
// A major revision of the vector/fixed array interface, unifying three independent
// libraries into one. The vecal (Vector Algebra), fish/vector (Multidimensional
// vector arrays) and the Eagle library (Geometric Algebra) thus share a common
// root now.
//
// Revision 1.4  2004/05/11 17:16:24  werner
// *** empty log message ***
//
// 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 __INTEGRATE4D_HPP
#define __INTEGRATE4D_HPP "Created 25.04.2004 11:02:10 by werner benger"

#include <ode/Integrator.hpp>
#include <eagle/FixedArray.hpp>
#include <ode/dop853decl.hpp>

#include "spacetime.h"

namespace Traum
{
        using namespace Eagle;

template <class real=long double, int Dims=4>
        class IntegratePathBase;

template <class Real, int N>
        struct  SecondOrderDiffEquationConstantSize
        {
                typedef Real real;

                struct RealArray_t : FixedArray<real, 2*N>             
                {
                        void init(int) {}
                };

                struct DOPVarsArray_t : FixedArray<dop_vars<real>, 2*N >
                {
                        void init(int) {}
                };

                IntegratePathBase<real, N>*parent;

                //virtual void Accel(real s, const real *x, const real *v, real *d2x_ds2) = 0;

                void DiffEqn(int nEquations, real s, const real *q, real *dq_ds)
                {
                const real *x = q;
                const real *v = q + N;
                real       *dx_ds = dq_ds;
                real       *dv_ds = dq_ds + N;

                        for(int i=0; i<N; i++)
                        {
                                dx_ds[i] = v[i];
                        }
                        parent->Accel(s, x, v, dv_ds);
                }
        };


template <>
class SPACETIME_API IntegratePathBase<long double, 4> : public Integrator<long double>
{
        typedef long double real;
        enum { n = 4 };

        typedef dop853< SecondOrderDiffEquationConstantSize<real, n> > PathIntegrator;
        PathIntegrator*myIntegrator;

public:
        IntegratePathBase<long double, 4>();
//      ~IntegratePathBase<long double, 4>();
        ~IntegratePathBase();

        real *initialize(int n, const real&s);
        override long  step_nr()      const;
        override const real&dx()      const;
        override const real&x0()      const;
        override const real&x1()      const;
        override int   size()         const;
        override const real&f(int i)  const;
        override       real f(int i, const real&x)  const;

        override success_code  advance(bool backward=false);

        virtual void Accel(real s, const real *x, const real *v, real *d2x_ds2) = 0;
};




template <class Real, int Dims>
class IntegratePath : public IntegratePathBase<Real, Dims>
{
public:
        typedef Real real;

        template <class sValue, class Value>
        void init(const sValue&s, const Value*q, const Value*dq)
        {
        real *x = initialize(2*Dims, s);
        real *v = x + Dims;
                for(int i=0; i<Dims; i++)
                {
                        x[i] =  q[i];
                        v[i] = dq[i];
                }
        }

        template <class sValue, class Value>
        void init(const sValue&s, const FixedArray<Value, Dims>&q, const FixedArray<Value, Dims>&dq)
        {
        real *x = initialize(2*Dims, s);
        real *v = x + Dims;
                for(int i=0; i<Dims; i++)
                {
                        x[i] =  q[i];
                        v[i] = dq[i];                   
                }
        }


        virtual void Accel(real s, const real *x, const real *v, real *d2x_ds2) = 0;

static  unsigned nPaths()
        {
                return Dims; 
        }

        const real&s0()  const
        {
                return this->x0(); 
        }

        const real&s1()  const
        {
                return this->x1(); 
        }

        const real&q(int i)  const
        {
                return this->f(i); 
        }

        real q(unsigned i, real s) const
        {
                return this->f(i, s ); 
        }

        const real&dq(int i)  const
        {
                return this->f(i+nPaths() ); 
        }

        real dq(unsigned i, real s) const
        {
                return this->f(i+nPaths(), s ); 
        }
};

typedef IntegratePath<long double, 4> IntegratePath4D;
//typedef IntegratePath<     double, 4> IntegratePath4Dd;
//typedef IntegratePath<     float , 4> IntegratePath4Df;


} // namespace Traum

#endif // __INTEGRATE4D_HPP