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();
//      ~IntegratePathBase<long double, 4>();
        ~IntegratePathBase();
 
        real *initialize(int n, const real&s);
        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) override;
 
        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