#include <RK4PreciseStep.h>
Public Member Functions | |
double | distance (const CartesianState &a, const CartesianState &b) const |
CartesianState | operator() (const CartesianState &start, const RKCartesianDerivative &deriv, double step, double eps) const |
std::pair< CartesianState, double > | stepWithAccuracy (const CartesianState &start, const RKCartesianDerivative &deriv, double step) const |
Private Member Functions | |
bool | verbose () const |
Definition at line 10 of file RK4PreciseStep.h.
double RK4PreciseStep::distance | ( | const CartesianState & | a, |
const CartesianState & | b | ||
) | const |
Definition at line 63 of file RK4PreciseStep.cc.
References mag(), CartesianState::momentum(), and CartesianState::position().
Referenced by stepWithAccuracy().
CartesianState RK4PreciseStep::operator() | ( | const CartesianState & | start, |
const RKCartesianDerivative & | deriv, | ||
double | step, | ||
double | eps | ||
) | const |
Definition at line 7 of file RK4PreciseStep.cc.
References gather_cfg::cout, max(), min, CartesianState::position(), funct::pow(), dqm_diff::start, launcher::step, stepWithAccuracy(), and verbose().
{ const double Safety = 0.9; double remainigStep = step; double stepSize = step; CartesianState currentStart = start; int nsteps = 0; std::pair<CartesianState, double> tryStep; do { tryStep = stepWithAccuracy( currentStart, deriv, stepSize); nsteps++; if (tryStep.second <eps) { if (remainigStep - stepSize < eps/2) { if (verbose()) std::cout << "Accuracy reached, and full step taken in " << nsteps << " steps" << std::endl; return tryStep.first; // we are there } else { remainigStep -= stepSize; // increase step size double factor = std::min( Safety * pow( fabs(eps/tryStep.second),0.2), 4.); stepSize = std::min( stepSize*factor, remainigStep); currentStart = tryStep.first; if (verbose()) std::cout << "Accuracy reached, but " << remainigStep << " remain after " << nsteps << " steps. Step size increased by " << factor << " to " << stepSize << std::endl; } } else { // decrease step size double factor = std::max( Safety * pow( fabs(eps/tryStep.second),0.25), 0.1); stepSize *= factor; if (verbose()) std::cout << "Accuracy not yet reached: delta = " << tryStep.second << ", step reduced by " << factor << " to " << stepSize << ", (R,z)= " << currentStart.position().perp() << ", " << currentStart.position().z() << std::endl; } } while (remainigStep > eps/2); return tryStep.first; }
std::pair< CartesianState, double > RK4PreciseStep::stepWithAccuracy | ( | const CartesianState & | start, |
const RKCartesianDerivative & | deriv, | ||
double | step | ||
) | const |
Definition at line 52 of file RK4PreciseStep.cc.
References diffTreeTool::diff, and distance().
Referenced by operator()().
{ RK4OneStep solver; CartesianState one(solver(start, deriv, step)); CartesianState firstHalf(solver(start, deriv, step/2)); CartesianState secondHalf(solver(firstHalf, deriv, step/2)); double diff = distance(one, secondHalf); return std::pair<CartesianState, double>(secondHalf,diff); }
bool RK4PreciseStep::verbose | ( | ) | const [private] |
Definition at line 68 of file RK4PreciseStep.cc.
Referenced by operator()().
{ // static bool verb = SimpleConfigurable<bool>(false,"RK4PreciseStep:verbose").value(); static bool verb = true; return verb; }