G4CashKarpRKF45.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4CashKarpRKF45.cc 69786 2013-05-15 09:38:51Z gcosmo $
00028 //
00029 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
00030 // order method (giving fifth-order accuracy) for the solution of an ODE.
00031 // Two different fourth order estimates are calculated; their difference
00032 // gives an error estimate. [ref. Numerical Recipes in C, 2nd Edition]
00033 // It is used to integrate the equations of the motion of a particle 
00034 // in a magnetic field.
00035 //
00036 //  [ref. Numerical Recipes in C, 2nd Edition]
00037 //
00038 // -------------------------------------------------------------------
00039 
00040 #include "G4CashKarpRKF45.hh"
00041 #include "G4LineSection.hh"
00042 
00044 //
00045 // Constructor
00046 
00047 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfMotion *EqRhs, 
00048                                  G4int noIntegrationVariables, 
00049                                  G4bool primary)
00050   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
00051     fLastStepLength(0.), fAuxStepper(0)
00052 {
00053   const G4int numberOfVariables = noIntegrationVariables;
00054 
00055   ak2 = new G4double[numberOfVariables] ;  
00056   ak3 = new G4double[numberOfVariables] ; 
00057   ak4 = new G4double[numberOfVariables] ; 
00058   ak5 = new G4double[numberOfVariables] ; 
00059   ak6 = new G4double[numberOfVariables] ; 
00060   ak7 = 0;
00061   yTemp = new G4double[numberOfVariables] ; 
00062   yIn = new G4double[numberOfVariables] ;
00063 
00064   fLastInitialVector = new G4double[numberOfVariables] ;
00065   fLastFinalVector = new G4double[numberOfVariables] ;
00066   fLastDyDx = new G4double[numberOfVariables];
00067 
00068   fMidVector = new G4double[numberOfVariables];
00069   fMidError =  new G4double[numberOfVariables];
00070   if( primary )
00071   { 
00072     fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
00073   }
00074 }
00075 
00077 //
00078 // Destructor
00079 
00080 G4CashKarpRKF45::~G4CashKarpRKF45()
00081 {
00082   delete[] ak2;
00083   delete[] ak3;
00084   delete[] ak4;
00085   delete[] ak5;
00086   delete[] ak6;
00087   // delete[] ak7;
00088   delete[] yTemp;
00089   delete[] yIn;
00090 
00091   delete[] fLastInitialVector;
00092   delete[] fLastFinalVector;
00093   delete[] fLastDyDx;
00094   delete[] fMidVector;
00095   delete[] fMidError; 
00096 
00097   delete fAuxStepper;
00098 }
00099 
00101 //
00102 // Given values for n = 6 variables yIn[0,...,n-1] 
00103 // known  at x, use the fifth-order Cash-Karp Runge-
00104 // Kutta-Fehlberg-4-5 method to advance the solution over an interval
00105 // Step and return the incremented variables as yOut[0,...,n-1]. Also
00106 // return an estimate of the local truncation error yErr[] using the
00107 // embedded 4th-order method. The user supplies routine
00108 // RightHandSide(y,dydx), which returns derivatives dydx for y .
00109 
00110 void
00111 G4CashKarpRKF45::Stepper(const G4double yInput[],
00112                          const G4double dydx[],
00113                                G4double Step,
00114                                G4double yOut[],
00115                                G4double yErr[])
00116 {
00117   // const G4int nvar = 6 ;
00118   // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
00119  G4int i;
00120 
00121  const G4double  b21 = 0.2 ,
00122                  b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
00123                  b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
00124 
00125                  b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
00126                  b54 = 35.0/27.0 ,
00127 
00128                  b61 = 1631.0/55296.0 , b62 =   175.0/512.0 ,
00129                  b63 =  575.0/13824.0 , b64 = 44275.0/110592.0 ,
00130                  b65 =  253.0/4096.0 ,
00131 
00132                  c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
00133                  c6 = 512.0/1771.0 ,
00134                                           dc5 = -277.0/14336.0 ;
00135 
00136  const G4double dc1 = c1 - 2825.0/27648.0 ,  dc3 = c3 - 18575.0/48384.0 ,
00137     dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
00138 
00139  // Initialise time to t0, needed when it is not updated by the integration.
00140  //        [ Note: Only for time dependent fields (usually electric) 
00141  //                  is it neccessary to integrate the time.] 
00142  yOut[7] = yTemp[7]   = yIn[7]; 
00143 
00144  const G4int numberOfVariables= this->GetNumberOfVariables(); 
00145  // The number of variables to be integrated over
00146 
00147    //  Saving yInput because yInput and yOut can be aliases for same array
00148 
00149    for(i=0;i<numberOfVariables;i++) 
00150    {
00151      yIn[i]=yInput[i];
00152    }
00153  // RightHandSide(yIn, dydx) ;              // 1st Step
00154 
00155  for(i=0;i<numberOfVariables;i++) 
00156  {
00157    yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
00158  }
00159  RightHandSide(yTemp, ak2) ;              // 2nd Step
00160 
00161  for(i=0;i<numberOfVariables;i++)
00162  {
00163     yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
00164  }
00165  RightHandSide(yTemp, ak3) ;              // 3rd Step
00166 
00167  for(i=0;i<numberOfVariables;i++)
00168  {
00169     yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
00170  }
00171  RightHandSide(yTemp, ak4) ;              // 4th Step
00172 
00173  for(i=0;i<numberOfVariables;i++)
00174  {
00175     yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
00176                       b54*ak4[i]) ;
00177  }
00178  RightHandSide(yTemp, ak5) ;              // 5th Step
00179 
00180  for(i=0;i<numberOfVariables;i++)
00181  {
00182     yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
00183                       b64*ak4[i] + b65*ak5[i]) ;
00184  }
00185  RightHandSide(yTemp, ak6) ;              // 6th Step
00186 
00187  for(i=0;i<numberOfVariables;i++)
00188  {
00189     // Accumulate increments with proper weights
00190 
00191     yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
00192 
00193     // Estimate error as difference between 4th and
00194     // 5th order methods
00195 
00196     yErr[i] = Step*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] +
00197               dc5*ak5[i] + dc6*ak6[i]) ;
00198 
00199     // Store Input and Final values, for possible use in calculating chord
00200     fLastInitialVector[i] = yIn[i] ;
00201     fLastFinalVector[i]   = yOut[i];
00202     fLastDyDx[i]          = dydx[i];
00203  }
00204  // NormaliseTangentVector( yOut ); // Not wanted
00205 
00206  fLastStepLength =Step;
00207 
00208  return ;
00209 } 
00210 
00212 
00213 void
00214 G4CashKarpRKF45::StepWithEst( const G4double*,
00215                               const G4double*,
00216                                     G4double,
00217                                     G4double*,
00218                                     G4double&,
00219                                     G4double&,
00220                               const G4double*,
00221                                     G4double*  )    
00222 {
00223   G4Exception("G4CashKarpRKF45::StepWithEst()", "GeomField0001",
00224               FatalException, "Method no longer used.");
00225   return ;
00226 }
00227 
00229 
00230 G4double  G4CashKarpRKF45::DistChord() const
00231 {
00232   G4double distLine, distChord; 
00233   G4ThreeVector initialPoint, finalPoint, midPoint;
00234 
00235   // Store last initial and final points (they will be overwritten in self-Stepper call!)
00236   initialPoint = G4ThreeVector( fLastInitialVector[0], 
00237                                 fLastInitialVector[1], fLastInitialVector[2]); 
00238   finalPoint   = G4ThreeVector( fLastFinalVector[0],  
00239                                 fLastFinalVector[1],  fLastFinalVector[2]); 
00240 
00241   // Do half a step using StepNoErr
00242 
00243   fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 
00244            fMidVector,   fMidError );
00245 
00246   midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);       
00247 
00248   // Use stored values of Initial and Endpoint + new Midpoint to evaluate
00249   //  distance of Chord
00250 
00251 
00252   if (initialPoint != finalPoint) 
00253   {
00254      distLine  = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
00255      distChord = distLine;
00256   }
00257   else
00258   {
00259      distChord = (midPoint-initialPoint).mag();
00260   }
00261   return distChord;
00262 }
00263 
00264 

Generated on Mon May 27 17:47:51 2013 for Geant4 by  doxygen 1.4.7