RungeKutta/src/RK4OneStep.cc

0001 #include "RK4OneStep.h"
0002 #include "RKCartesianDerivative.h"
0003
0004 #include <iostream>
0005
0006 CartesianState RK4OneStep::operator()(const CartesianState& start,
0007                                       const RKCartesianDerivative& deriv,
0008                                       double step) const {
0009   double s0 = 0;  // derivatives do't depend on absolute value of the integration variable
0010   CartesianState k1 = step * deriv(s0, start);
0011   CartesianState k2 = step * deriv(s0 + step / 2, start + k1 / 2);
0012   CartesianState k3 = step * deriv(s0 + step / 2, start + k2 / 2);
0013   CartesianState k4 = step * deriv(s0 + step, start + k3);
0014   /*
0015   std::cout << "k1 = " << k1.position() << k1.momentum() << std::endl;
0016   std::cout << "k2 = " << k2.position() << k2.momentum() << std::endl;
0017   std::cout << "k3 = " << k3.position() << k3.momentum() << std::endl;
0018   std::cout << "k4 = " << k4.position() << k4.momentum() << std::endl;
0019 */
0020   CartesianState result = start + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6;
0021   return result;
0022 }