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