CLHEP/GenericFunctions/AdaptiveRKStepper.hh

0001 #ifndef _AdaptiveRKStepper_h_
0002 #define _AdaptiveRKStepper_h_
0003 #include "CLHEP/GenericFunctions/RKIntegrator.hh"
0004 //
0005 // This is the default RKStepper routine, used within the RKIntegrator
0006 // when no other stepper is provided.
0007 //
0008 namespace Genfun {
0009
0010   class AdaptiveRKStepper:public RKIntegrator::RKStepper {
0011
0012   public:
0013
0014     class EEStepper;
0015
0016     // Constructor. A default stepper is created if none
0017     // is specified--it is an Embedded Runge Kutta stepper
0018     // using a Cash-Karp extended Butcher Tableau, which
0019     // is fourth-order.
0020     AdaptiveRKStepper(const EEStepper *eeStepper=NULL);
0021
0022     // Copy constructor:
0023     AdaptiveRKStepper(const AdaptiveRKStepper & right);
0024
0025     // Destructor:
0026     virtual ~AdaptiveRKStepper();
0027
0028     // Takes *multiple* steps if needed up to or exceeding the
0029     // specified time.  This called from the RKFunction and not
0030     // by users.    If a time limit is specified, only one step
0031     // is taken--this mode is used to interpolate between cached
0032     // meshpoints.
0033
0034     // Input:  data contains diffEqn information and cache:
0035     //         sData a summary of starting point information
0036     //         timeLimit (if nonzero) time limit for step
0037     // Output  dData a summary of into at the ending point.
0038     //
0039
0040     virtual void step (const RKIntegrator::RKData       * data,
0041                const RKIntegrator::RKData::Data & sdata,
0042                RKIntegrator::RKData::Data       & ddata,
0043                double                             timeLimit) const ;
0044     // Clone
0045     virtual AdaptiveRKStepper *clone() const;
0046
0047     // Accessors and modifiers to algorithmic parameters.  Roughly
0048     // speaking these are ordered according to importance: the user
0049     // will often wish to modify the tolerance and the starting
0050     // stepsize, but rarely should need to touch any of the others.
0051
0052     // The tolerance:
0053     double & tolerance();
0054     const double & tolerance() const;
0055
0056     // The starting stepsize:
0057     double & startingStepsize();
0058     const double & startingStepsize() const;
0059
0060     // The safety factor.  Step size increases are moderated by this
0061     // factor:
0062     double & safetyFactor();
0063     const double & safetyFactor() const;
0064
0065     // The minimum amount by which a step size is decreased:
0066     double & rmin();
0067     const double & rmin() const;
0068
0069     // The maximum amount by which a step size is increased:
0070     double & rmax();
0071     const double & rmax() const;
0072
0073
0074   private:
0075
0076     // It is illegal to assign an AdaptiveRKStepper:
0077     AdaptiveRKStepper & operator=(const AdaptiveRKStepper & right);
0078
0079     const        EEStepper   *eeStepper;
0080     double       T;
0081     double       sStepsize;
0082     double       S;
0083     double       Rmin;
0084     double       Rmax;
0085     mutable      double       stepsize;
0086   };
0087
0088
0089   //
0090   // An abstract base class for steppers that return an error
0091   // estimate at each step:
0092   //
0093   class AdaptiveRKStepper::EEStepper {
0094   public:
0095
0096     virtual ~EEStepper();
0097     virtual void step   (const RKIntegrator::RKData       * data,
0098              const RKIntegrator::RKData::Data & sdata,
0099              RKIntegrator::RKData::Data       & ddata,
0100              std::vector<double> & errors) const = 0;
0101     virtual EEStepper *clone() const=0;
0102     virtual unsigned int order() const=0;
0103   };
0104 }
0105
0106 #endif