Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4CashKarpRKF45.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // G4CashKarpRKF45 implementation
 27 //
 28 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
 29 // order method (giving fifth-order accuracy) for the solution of an ODE.
 30 // Two different fourth order estimates are calculated; their difference
 31 // gives an error estimate. [ref. Numerical Recipes in C, 2nd Edition]
 32 // It is used to integrate the equations of the motion of a particle 
 33 // in a magnetic field.
 34 //
 35 // [ref. Numerical Recipes in C, 2nd Edition]
 36 //
 37 // Authors: J.Apostolakis, V.Grichine - 30.01.1997
 38 // -------------------------------------------------------------------
 39 
 40 #include "G4CashKarpRKF45.hh"
 41 #include "G4LineSection.hh"
 42 
 43 /////////////////////////////////////////////////////////////////////
 44 //
 45 // Constructor
 46 //
 47 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfMotion *EqRhs, 
 48                                  G4int noIntegrationVariables, 
 49                                  G4bool primary)
 50   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
 51 {
 52   const G4int numberOfVariables =
 53       std::max( noIntegrationVariables,
 54                ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
 55   // For better alignment with cache-line
 56   
 57   ak2 = new G4double[numberOfVariables] ;
 58   ak3 = new G4double[numberOfVariables] ;
 59   ak4 = new G4double[numberOfVariables] ;
 60   ak5 = new G4double[numberOfVariables] ;
 61   ak6 = new G4double[numberOfVariables] ;
 62   // ak7 = 0;
 63 
 64   // Must ensure space extra 'state' variables exists - i.e. yIn[7]
 65   const G4int numStateMax  = std::max(GetNumberOfStateVariables(), 8);  
 66   const G4int numStateVars = std::max(noIntegrationVariables,
 67                                       numStateMax );
 68                                    // GetNumberOfStateVariables() ); 
 69                                       
 70   yTemp = new G4double[numStateVars] ;
 71   yIn = new G4double[numStateVars] ;
 72 
 73   fLastInitialVector = new G4double[numStateVars] ;
 74   fLastFinalVector = new G4double[numStateVars] ;
 75   fLastDyDx = new G4double[numberOfVariables];
 76 
 77   fMidVector = new G4double[numStateVars];
 78   fMidError =  new G4double[numStateVars];
 79   if( primary )
 80   { 
 81     fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
 82   }
 83 }
 84 
 85 /////////////////////////////////////////////////////////////////////
 86 //
 87 // Destructor
 88 //
 89 G4CashKarpRKF45::~G4CashKarpRKF45()
 90 {
 91   delete [] ak2;
 92   delete [] ak3;
 93   delete [] ak4;
 94   delete [] ak5;
 95   delete [] ak6;
 96   // delete [] ak7;
 97   delete [] yTemp;
 98   delete [] yIn;
 99 
100   delete [] fLastInitialVector;
101   delete [] fLastFinalVector;
102   delete [] fLastDyDx;
103   delete [] fMidVector;
104   delete [] fMidError; 
105 
106   delete fAuxStepper;
107 }
108 
109 //////////////////////////////////////////////////////////////////////
110 //
111 // Given values for n = 6 variables yIn[0,...,n-1] 
112 // known  at x, use the fifth-order Cash-Karp Runge-
113 // Kutta-Fehlberg-4-5 method to advance the solution over an interval
114 // Step and return the incremented variables as yOut[0,...,n-1]. Also
115 // return an estimate of the local truncation error yErr[] using the
116 // embedded 4th-order method. The user supplies routine
117 // RightHandSide(y,dydx), which returns derivatives dydx for y .
118 //
119 void
120 G4CashKarpRKF45::Stepper(const G4double yInput[],
121                          const G4double dydx[],
122                                G4double Step,
123                                G4double yOut[],
124                                G4double yErr[])
125 {
126   // const G4int nvar = 6 ;
127   // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
128   G4int i;
129 
130   const G4double  b21 = 0.2 ,
131                   b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
132                   b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
133 
134                   b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
135                   b54 = 35.0/27.0 ,
136 
137                   b61 = 1631.0/55296.0 , b62 =   175.0/512.0 ,
138                   b63 =  575.0/13824.0 , b64 = 44275.0/110592.0 ,
139                   b65 =  253.0/4096.0 ,
140 
141                   c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
142                   c6 = 512.0/1771.0 , dc5 = -277.0/14336.0 ;
143 
144   const G4double dc1 = c1 - 2825.0/27648.0 ,  dc3 = c3 - 18575.0/48384.0 ,
145                  dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
146 
147   // Initialise time to t0, needed when it is not updated by the integration.
148   //        [ Note: Only for time dependent fields (usually electric) 
149   //                  is it neccessary to integrate the time.] 
150   yOut[7] = yTemp[7] = yIn[7] = yInput[7]; 
151 
152   const G4int numberOfVariables= this->GetNumberOfVariables(); 
153     // The number of variables to be integrated over
154 
155   //  Saving yInput because yInput and yOut can be aliases for same array
156 
157   for(i=0; i<numberOfVariables; ++i) 
158   {
159     yIn[i]=yInput[i];
160   }
161   // RightHandSide(yIn, dydx) ;            // 1st Step
162 
163   for(i=0; i<numberOfVariables; ++i) 
164   {
165     yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
166   }
167   RightHandSide(yTemp, ak2) ;              // 2nd Step
168 
169   for(i=0; i<numberOfVariables; ++i)
170   {
171     yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
172   }
173   RightHandSide(yTemp, ak3) ;              // 3rd Step
174 
175   for(i=0; i<numberOfVariables; ++i)
176   {
177     yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
178   }
179   RightHandSide(yTemp, ak4) ;              // 4th Step
180 
181   for(i=0; i<numberOfVariables; ++i)
182   {
183     yTemp[i] = yIn[i] + Step*(b51*dydx[i]
184                       + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]) ;
185   }
186   RightHandSide(yTemp, ak5) ;              // 5th Step
187 
188   for(i=0; i<numberOfVariables; ++i)
189   {
190     yTemp[i] = yIn[i] + Step*(b61*dydx[i]
191                       + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
192   }
193   RightHandSide(yTemp, ak6) ;              // 6th Step
194 
195   for(i=0; i<numberOfVariables; ++i)
196   {
197     // Accumulate increments with proper weights
198     //
199     yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
200 
201     // Estimate error as difference between 4th and 5th order methods
202     //
203     yErr[i] = Step*(dc1*dydx[i]
204             + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i]) ;
205 
206     // Store Input and Final values, for possible use in calculating chord
207     //
208     fLastInitialVector[i] = yIn[i] ;
209     fLastFinalVector[i]   = yOut[i];
210     fLastDyDx[i]          = dydx[i];
211   }
212   // NormaliseTangentVector( yOut ); // Not wanted
213 
214   fLastStepLength = Step;
215 
216   return;
217 } 
218 
219 ///////////////////////////////////////////////////////////////////////////////
220 //
221 void
222 G4CashKarpRKF45::StepWithEst( const G4double*,
223                               const G4double*,
224                                     G4double,
225                                     G4double*,
226                                     G4double&,
227                                     G4double&,
228                               const G4double*,
229                                     G4double*  )    
230 {
231   G4Exception("G4CashKarpRKF45::StepWithEst()", "GeomField0001",
232               FatalException, "Method no longer used.");
233   return ;
234 }
235 
236 /////////////////////////////////////////////////////////////////
237 //
238 G4double  G4CashKarpRKF45::DistChord() const
239 {
240   G4double distLine, distChord; 
241   G4ThreeVector initialPoint, finalPoint, midPoint;
242 
243   // Store last initial and final points
244   // (they will be overwritten in self-Stepper call!)
245   //
246   initialPoint = G4ThreeVector( fLastInitialVector[0], 
247                                 fLastInitialVector[1], fLastInitialVector[2]); 
248   finalPoint   = G4ThreeVector( fLastFinalVector[0],  
249                                 fLastFinalVector[1],  fLastFinalVector[2]); 
250 
251   // Do half a step using StepNoErr
252   //
253   fAuxStepper->Stepper( fLastInitialVector, fLastDyDx,
254                         0.5 * fLastStepLength, fMidVector, fMidError );
255 
256   midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);       
257 
258   // Use stored values of Initial and Endpoint + new Midpoint to evaluate
259   // distance of Chord
260   //
261   if (initialPoint != finalPoint) 
262   {
263      distLine  = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
264      distChord = distLine;
265   }
266   else
267   {
268      distChord = (midPoint-initialPoint).mag();
269   }
270   return distChord;
271 }
272