RK4PreciseStep.cc

Go to the documentation of this file.

#include "RK4PreciseStep.h"
#include "RK4OneStep.h"
//#include "Utilities/UI/interface/SimpleConfigurable.h"
#include <iostream>

CartesianState
RK4PreciseStep::operator()( const CartesianState& start, const RKCartesianDerivative& deriv,
                double step, double eps) const
{
    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
{
    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);
}

double RK4PreciseStep::distance( const CartesianState& a, const CartesianState& b) const
{
    return (a.position() - b.position()).mag() + (a.momentum() - b.momentum()).mag() / b.momentum().mag();
}

bool RK4PreciseStep::verbose() const
{
  return true;
}