RungeKuttaFehlberg.hpp

Go to the documentation of this file.

#pragma ident "$Id: RungeKuttaFehlberg.hpp 2457 2010-08-18 14:20:12Z coandrei $"

#ifndef GPSTK_RUNGE_KUTTA_FEHLBERG_HPP
#define GPSTK_RUNGE_KUTTA_FEHLBERG_HPP


//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  Wei Yan - Chinese Academy of Sciences . 2009, 2010
//
//============================================================================


#include "Integrator.hpp"


namespace gpstk
{


   class RungeKuttaFehlberg : public Integrator
   {
   public:
         // Parametera for RKF78
      struct RKF78Param
      {
         double a[13] ;
         double b[13][12] ;
         double c1[13] ;
         double c2[13] ;
      };

   public:

      RungeKuttaFehlberg();

      virtual ~RungeKuttaFehlberg()
      { };

      void test();

      virtual Vector<double> integrateTo(const double&           t, 
                                         const Vector<double>&   y, 
                                         EquationOfMotion*       peom,
                                         const double&           tf );
      
      RungeKuttaFehlberg& setAccuracy(const double& accuracy)
      { accuracyEps = accuracy; return (*this);}
         
      RungeKuttaFehlberg& setMinStepSize(const double& step)
      { minStepSize = step; return (*this ); }
      
      RungeKuttaFehlberg& setAdaptive(const bool& adaptive = true)
      { isAdaptive = adaptive; return (*this); }


   protected:
         

      Vector<double> integrateFixedStep(const double&           t, 
                                        const Vector<double>&   y, 
                                        EquationOfMotion*       peom,
                                        const double&           tf );

      Vector<double> integrateAdaptive(const double&           t, 
                                       const Vector<double>&   y, 
                                       EquationOfMotion*       peom,
                                       const double&           tf );
      

         // takes one "quality-controlled" Runge-Kutta-Fehlberg step 
         // 0 = Success
         // 1 = Unable to allocate workspace memory
         // 2 = Stepsize underflow
      int rkfqcs(double&           x,
                 Vector<double>&   y,
                 const double&     htry,
                 const double&     accuracy,
                 EquationOfMotion* peom,
                 Vector<double>&   yscal,
                 double& hdid,
                 double& hnext);
         
         // RKF78 single step
         // 0 = Success
         // 1 = Failed to allocate memory
      int rkfs78(const double&          x,
                 const Vector<double>&  y,
                 const double&          h,
                 EquationOfMotion*      peom,
                 Vector<double>&        yout,
                 Vector<double>&        yerr);

   protected:
      
      double accuracyEps;

      double minStepSize;

      bool isAdaptive;

      static const double RKF_EPS;

      static const double RKF_MAXSTEP;

   private:
     
      const static struct RKF78Param rkf78_param;
      
         // Easy accessing RKF78 parameters
      double A(int i){ return rkf78_param.a[i]; }

      double B(int i, int j){ return rkf78_param.b[i][j]; }
      
      double C(int i){ return rkf78_param.c1[i]; }

      double C2(int i){ return rkf78_param.c2[i]; }
      
      double DC(int i){ return (rkf78_param.c1[i]-rkf78_param.c2[i]); }


   }; // End of class 'RungeKuttaFehlberg'

      // @}

}  // End of namespace 'gpstk'


#endif // GPSTK_RUNGE_KUTTA_FEHLBERG_HPP

---