include/Geant4/G4DormandPrince745.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4DormandPrince745
0027 //
0028 // Class desription:
0029 //
0030 // An implementation of the 5th order embedded RK method from the paper:
0031 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0032 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
0033 //
0034 // DormandPrince7 - 5(4) embedded RK method
0035
0036 // Author: Somnath Banerjee (CERN, Google Summer of Code 2015), 25.05.2015
0037 // Supervision: John Apostolakis (CERN)
0038 // --------------------------------------------------------------------
0039 #ifndef G4DORMAND_PRINCE_745_HH
0040 #define G4DORMAND_PRINCE_745_HH
0041
0042 #include "G4MagIntegratorStepper.hh"
0043 #include "G4FieldUtils.hh"
0044
0045 /**
0046  * @brief G4DormandPrince745 implements the 5th order embedded Runge-Kutta
0047  * method, non-FSAL definition of the stepper() method that evaluates one step
0048  * in field propagation.
0049  */
0050
0051 class G4DormandPrince745 : public G4MagIntegratorStepper
0052 {
0053   public:
0054
0055     /**
0056      * Constructor for G4DormandPrince745.
0057      *  @param[in] equation Pointer to the provided equation of motion.
0058      *  @param[in] numberOfVariables The number of integration variables.
0059      */
0060     G4DormandPrince745(G4EquationOfMotion* equation,
0061                        G4int numberOfVariables = 6);
0062
0063     /**
0064      * Default Destructor.
0065      */
0066     ~G4DormandPrince745() override = default;
0067
0068     /**
0069      * Copy constructor and assignment operator not allowed.
0070      */
0071     G4DormandPrince745(const G4DormandPrince745&) = delete;
0072     G4DormandPrince745& operator=(const G4DormandPrince745&) = delete;
0073
0074     /**
0075      * The stepper for the Runge Kutta integration.
0076      * The stepsize is fixed, with the step size given by 'hstep'.
0077      * Integrates ODE starting values yInput[0 to 6].
0078      * Outputs yOutput[] and its estimated error yError[].
0079      *  @param[in] yInput Starting values array of integration variables.
0080      *  @param[in] dydx Derivatives array.
0081      *  @param[in] hstep The given step size.
0082      *  @param[out] yOutput Integration output.
0083      *  @param[out] yError The estimated error.
0084      */
0085     void Stepper(const G4double yInput[],
0086                  const G4double dydx[],
0087                        G4double hstep,
0088                        G4double yOutput[],
0089                        G4double yError[]) override;
0090
0091     /**
0092      * Same as the Stepper() function above, with dydx also in ouput.
0093      *  @param[in] yInput Starting values array of integration variables.
0094      *  @param[in] dydx Derivatives array.
0095      *  @param[in] hstep The given step size.
0096      *  @param[out] yOutput Integration output.
0097      *  @param[out] yError The estimated error.
0098      *  @param[out] dydxOutput dysx in output.
0099      */
0100     void Stepper(const G4double yInput[],
0101                  const G4double dydx[],
0102                        G4double hstep,
0103                        G4double yOutput[],
0104                        G4double yError[],
0105                        G4double dydxOutput[]);
0106
0107
0108     /**
0109      * Interface method for interpolation setup. Does nothing here.
0110      */
0111     inline void SetupInterpolation() {}
0112
0113     /**
0114      * Calculates the output at the tau fraction of Step.
0115      * Lower (4th) order interpolant given by Dormand and Prince.
0116      */
0117     void Interpolate4thOrder(G4double yOut[], G4double tau) const;
0118
0119     /**
0120      * Wrapper for Interpolate4thOrder() function above.
0121      */
0122     inline void Interpolate(G4double tau, G4double yOut[]) const
0123     {
0124       Interpolate4thOrder(yOut, tau);
0125     }
0126
0127     /**
0128      * Sets up the extra stages for the 5th order interpolant.
0129      */
0130     void SetupInterpolation5thOrder();
0131
0132     /**
0133      * Calculates the interpolated result 'yOut' with the coefficients.
0134      * Interpolant of 5th order given by Baker, Dormand, Gilmore and Prince.
0135      */
0136     void Interpolate5thOrder(G4double yOut[], G4double tau) const;
0137
0138     /**
0139      * Returns the distance from chord line.
0140      */
0141     G4double DistChord() const override;
0142
0143     /**
0144      * Returns the order, 4, of integration.
0145      */
0146     inline G4int IntegratorOrder() const override { return 4; }
0147
0148     /**
0149      * Returns the stepper type-ID, "kDormandPrince745".
0150      */
0151     inline G4StepperType StepperType() const override { return kDormandPrince745; }
0152
0153     /**
0154      * Methods to return the stepper name and description.
0155      */
0156     const G4String& StepperTypeName() const;
0157     const G4String& StepperDescription() const;
0158
0159     /**
0160      * Returns the field state in output.
0161      */
0162     inline const field_utils::State& GetYOut() const { return fyOut; }
0163
0164     /**
0165      * Returns a pointer to the equation of motion.
0166      */
0167     inline G4EquationOfMotion* GetSpecificEquation() { return GetEquationOfMotion(); }
0168
0169   private:
0170
0171     field_utils::State ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
0172     field_utils::State fyIn, fyOut, fdydxIn;
0173
0174     G4double fLastStepLength = -1.0;
0175 };
0176
0177 #endif