include/Geant4/G4TMagErrorStepper.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 // G4TMagErrorStepper
0027 //
0028 // Class description:
0029 //
0030 // Templated version of G4MagErrorStepper
0031 //
0032 //
0033 // Created: Josh Xie  (supported by Google Summer of Code 2014 )
0034 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0035 // Adapted from G4G4TMagErrorStepper class
0036 // --------------------------------------------------------------------
0037 #ifndef G4TMAG_ERROR_STEPPER_HH
0038 #define G4TMAG_ERROR_STEPPER_HH
0039
0040 #include "G4Types.hh"
0041 #include "G4MagIntegratorStepper.hh"
0042 #include "G4ThreeVector.hh"
0043 #include "G4LineSection.hh"
0044
0045 template <class T_Stepper, class T_Equation, unsigned int N>
0046 class G4TMagErrorStepper : public G4MagIntegratorStepper
0047 {
0048  public:  // with description
0049   G4TMagErrorStepper(T_Equation* EqRhs, G4int numberOfVariables,
0050                      G4int numStateVariables = 12)
0051     : G4MagIntegratorStepper(EqRhs, numberOfVariables, numStateVariables)
0052     , fEquation_Rhs(EqRhs)
0053   {
0054     // G4int nvar = std::max(this->GetNumberOfVariables(), 8);
0055   }
0056
0057   virtual ~G4TMagErrorStepper() { ; }
0058
0059   inline void RightHandSide(G4double y[], G4double dydx[])
0060   {
0061     fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0062   }
0063
0064   inline void Stepper(const G4double yInput[], const G4double dydx[],
0065                       G4double hstep, G4double yOutput[], G4double yError[]) override final;
0066
0067   inline G4double DistChord() const override final;
0068
0069  private:
0070   G4TMagErrorStepper(const G4TMagErrorStepper&);
0071   G4TMagErrorStepper& operator=(const G4TMagErrorStepper&);
0072   // Private copy constructor and assignment operator.
0073
0074  private:
0075   // STATE
0076   G4ThreeVector fInitialPoint, fMidPoint, fFinalPoint;
0077   // Data stored in order to find the chord
0078
0079   // Dependent Objects, owned --- part of the STATE
0080   G4double yInitial[N < 8 ? 8 : N];
0081   G4double yMiddle[N < 8 ? 8 : N];
0082   G4double dydxMid[N < 8 ? 8 : N];
0083   G4double yOneStep[N < 8 ? 8 : N];
0084   // The following arrays are used only for temporary storage
0085   // they are allocated at the class level only for efficiency -
0086   // so that calls to new and delete are not made in Stepper().
0087
0088   T_Equation* fEquation_Rhs;
0089 };
0090
0091 // ------------   Implementation -----------------------
0092
0093 template <class T_Stepper, class T_Equation, unsigned int N >
0094 void G4TMagErrorStepper<T_Stepper,T_Equation,N>::
0095 Stepper(const G4double yInput[],
0096         const G4double dydx[],
0097               G4double hstep,
0098               G4double yOutput[],
0099               G4double yError[])
0100 // The stepper for the Runge Kutta integration. The stepsize
0101 // is fixed, with the Step size given by hstep.
0102 // Integrates ODE starting values y[0 to N].
0103 // Outputs yout[] and its estimated error yerr[].
0104 {
0105   const unsigned int maxvar = GetNumberOfStateVariables();
0106
0107   //  Saving yInput because yInput and yOutput can be aliases for same array
0108   for(unsigned int i = 0; i < N; ++i)
0109      yInitial[i] = yInput[i];
0110   yInitial[7] =
0111      yInput[7];  // Copy the time in case ... even if not really needed
0112   yMiddle[7]  = yInput[7];  // Copy the time from initial value
0113   yOneStep[7] = yInput[7];  // As it contributes to final value of yOutput ?
0114   // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
0115   for(unsigned int i = N; i < maxvar; ++i)
0116      yOutput[i] = yInput[i];
0117
0118   G4double halfStep = hstep * 0.5;
0119
0120   // Do two half steps
0121
0122   static_cast<T_Stepper*>(this)->DumbStepper(yInitial, dydx, halfStep,
0123                                                yMiddle);
0124   this->RightHandSide(yMiddle, dydxMid);
0125   static_cast<T_Stepper*>(this)->DumbStepper(yMiddle, dydxMid, halfStep,
0126                                              yOutput);
0127
0128   // Store midpoint, chord calculation
0129
0130   fMidPoint = G4ThreeVector(yMiddle[0], yMiddle[1], yMiddle[2]);
0131
0132   // Do a full Step
0133   static_cast<T_Stepper*>(this)->DumbStepper(yInitial, dydx, hstep, yOneStep);
0134   for(unsigned int i = 0; i < N; ++i)
0135   {
0136     yError[i] = yOutput[i] - yOneStep[i];
0137     yOutput[i] +=
0138        yError[i] *
0139        T_Stepper::IntegratorCorrection;  // Provides accuracy increased
0140     // by 1 order via the
0141     // Richardson Extrapolation
0142   }
0143
0144   fInitialPoint = G4ThreeVector(yInitial[0], yInitial[1], yInitial[2]);
0145   fFinalPoint   = G4ThreeVector(yOutput[0], yOutput[1], yOutput[2]);
0146
0147   return;
0148 }
0149
0150
0151 template <class T_Stepper, class T_Equation, unsigned int N >
0152 inline G4double
0153 G4TMagErrorStepper<T_Stepper,T_Equation,N>::DistChord() const
0154 {
0155   // Estimate the maximum distance from the curve to the chord
0156   //
0157   //  We estimate this using the distance of the midpoint to
0158   //  chord (the line between
0159   //
0160   //  Method below is good only for angle deviations < 2 pi,
0161   //   This restriction should not a problem for the Runge cutta methods,
0162   //   which generally cannot integrate accurately for large angle deviations.
0163   G4double distLine, distChord;
0164
0165   if(fInitialPoint != fFinalPoint)
0166   {
0167     distLine = G4LineSection::Distline(fMidPoint, fInitialPoint, fFinalPoint);
0168     // This is a class method that gives distance of Mid
0169     //  from the Chord between the Initial and Final points.
0170
0171     distChord = distLine;
0172   }
0173   else
0174   {
0175      distChord = (fMidPoint - fInitialPoint).mag();
0176   }
0177
0178   return distChord;
0179 }
0180
0181 #endif /* G4TMagErrorStepper_HH */