RungeKutta/src/RKOne4OrderStep.h

0001 #ifndef RKOne4OrderStep_H
0002 #define RKOne4OrderStep_H
0003
0004 #include "FWCore/Utilities/interface/Visibility.h"
0005 #include "RKDistance.h"
0006 #include "RK4OneStepTempl.h"
0007
0008 #include <utility>
0009
0010 template <typename T, int N>
0011 class dso_internal RKOne4OrderStep {
0012 public:
0013   typedef T Scalar;
0014   typedef RKSmallVector<T, N> Vector;
0015
0016   std::pair<Vector, T> operator()(Scalar startPar,
0017                                   const Vector& startState,
0018                                   const RKDerivative<T, N>& deriv,
0019                                   const RKDistance<T, N>& dist,
0020                                   Scalar step) {
0021     const Scalar huge = 1.e5;  // ad hoc protection against infinities, must be done better!
0022     const Scalar hugediff = 100.;
0023
0024     RK4OneStepTempl<T, N> solver;
0025     Vector one(solver(startPar, startState, deriv, step));
0026     if (std::abs(one[0]) > huge || std::abs(one(1)) > huge)
0027       return std::pair<Vector, Scalar>(one, hugediff);
0028
0029     Vector firstHalf(solver(startPar, startState, deriv, step / 2));
0030     Vector secondHalf(solver(startPar + step / 2, firstHalf, deriv, step / 2));
0031     Scalar diff = dist(one, secondHalf, startPar + step);
0032     return std::pair<Vector, Scalar>(secondHalf, diff);
0033   }
0034 };
0035 #endif