TangentialDifferentialEquation.hpp

#ifndef __FISH_CEPHALUS_TangentialDifferentialEquation_HPP
#define __FISH_CEPHALUS_TangentialDifferentialEquation_HPP

#include <aerie/AnalyticFunction.hpp>
#include <ode/dop853def.hpp>
#include <field/Field.hpp>

namespace Fiber
{


struct  TangentialDifferentialEquationAnalytic
{
        typedef ::Eagle::tvector3 tvector3;

        typedef Eagle::PhysicalSpace::AnalyticFunction<tvector3> Eqn_t;

        const Eqn_t*Eqn;

        enum { N = 3 };
        typedef double real;

        double  Time;
        
        TangentialDifferentialEquationAnalytic()
        : Eqn(0)
        , Time(0.0)
        {}

        struct RealArray_t : Eagle::FixedArray<real, N> 
        {
                void init(int) {}
        };
        
        struct DOPVarsArray_t : Eagle::FixedArray<Traum::dop_vars<real>, N >
        {
                void init(int) {}
        };
        
        void DiffEqn(int nEquations, real s, const real *q, real *dq_ds) const
        {
        Eagle::point3   Q;
                for(int k=0; k<N; k++)
                        Q[k] = q[k]; 

                assert( Eqn ); 
        tvector3 DataDir = (*Eqn)(Time, Q);

                for(int i=0; i<N; i++)
                {
                        dq_ds[i] = DataDir[i];
                }
        }
};

typedef Traum::dop853<  TangentialDifferentialEquationAnalytic> LineIntegratorAnalytic;


template <class LineIntegratorAnalyticContainer>
bool    setupAnalyticIntegrators(LineIntegratorAnalyticContainer&AnalyticLineIntegrators, const RefPtr<Field>&ToIntegrateField, double Time)
{
        if (!ToIntegrateField)
                return false; 

        if (MemCore::InterfacePtr<TangentialDifferentialEquationAnalytic::Eqn_t>
                    AnalyticField = *ToIntegrateField)
        {
                for(typename LineIntegratorAnalyticContainer::iterator 
                            it = AnalyticLineIntegrators.begin(); 
                            it != AnalyticLineIntegrators.end(); it++)
                {
                TangentialDifferentialEquationAnalytic::Eqn_t&AF = *AnalyticField;
                        (*it).Eqn = &AF; 
                        // note: need to set the time at each integration point in the following code 
                        (*it).Time = Time;
                } 
                return true;
        } 
        return false;
}


} // namespace Fiber

#endif // __FISH_CEPHALUS_TangentialDifferentialEquation_HPP