Geant4_10
G4CashKarpRKF45.cc
Go to the documentation of this file.
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 //
27 // $Id: G4CashKarpRKF45.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
30 // order method (giving fifth-order accuracy) for the solution of an ODE.
31 // Two different fourth order estimates are calculated; their difference
32 // gives an error estimate. [ref. Numerical Recipes in C, 2nd Edition]
33 // It is used to integrate the equations of the motion of a particle
34 // in a magnetic field.
35 //
36 // [ref. Numerical Recipes in C, 2nd Edition]
37 //
38 // -------------------------------------------------------------------
39 
40 #include "G4CashKarpRKF45.hh"
41 #include "G4LineSection.hh"
42 
44 //
45 // Constructor
46 
48  G4int noIntegrationVariables,
49  G4bool primary)
50  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
51  fLastStepLength(0.), fAuxStepper(0)
52 {
53  const G4int numberOfVariables = noIntegrationVariables;
54 
55  ak2 = new G4double[numberOfVariables] ;
56  ak3 = new G4double[numberOfVariables] ;
57  ak4 = new G4double[numberOfVariables] ;
58  ak5 = new G4double[numberOfVariables] ;
59  ak6 = new G4double[numberOfVariables] ;
60  ak7 = 0;
61  yTemp = new G4double[numberOfVariables] ;
62  yIn = new G4double[numberOfVariables] ;
63 
64  fLastInitialVector = new G4double[numberOfVariables] ;
65  fLastFinalVector = new G4double[numberOfVariables] ;
66  fLastDyDx = new G4double[numberOfVariables];
67 
68  fMidVector = new G4double[numberOfVariables];
69  fMidError = new G4double[numberOfVariables];
70  if( primary )
71  {
72  fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
73  }
74 }
75 
77 //
78 // Destructor
79 
81 {
82  delete[] ak2;
83  delete[] ak3;
84  delete[] ak4;
85  delete[] ak5;
86  delete[] ak6;
87  // delete[] ak7;
88  delete[] yTemp;
89  delete[] yIn;
90 
91  delete[] fLastInitialVector;
92  delete[] fLastFinalVector;
93  delete[] fLastDyDx;
94  delete[] fMidVector;
95  delete[] fMidError;
96 
97  delete fAuxStepper;
98 }
99 
101 //
102 // Given values for n = 6 variables yIn[0,...,n-1]
103 // known at x, use the fifth-order Cash-Karp Runge-
104 // Kutta-Fehlberg-4-5 method to advance the solution over an interval
105 // Step and return the incremented variables as yOut[0,...,n-1]. Also
106 // return an estimate of the local truncation error yErr[] using the
107 // embedded 4th-order method. The user supplies routine
108 // RightHandSide(y,dydx), which returns derivatives dydx for y .
109 
110 void
112  const G4double dydx[],
113  G4double Step,
114  G4double yOut[],
115  G4double yErr[])
116 {
117  // const G4int nvar = 6 ;
118  // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
119  G4int i;
120 
121  const G4double b21 = 0.2 ,
122  b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
123  b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
124 
125  b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
126  b54 = 35.0/27.0 ,
127 
128  b61 = 1631.0/55296.0 , b62 = 175.0/512.0 ,
129  b63 = 575.0/13824.0 , b64 = 44275.0/110592.0 ,
130  b65 = 253.0/4096.0 ,
131 
132  c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
133  c6 = 512.0/1771.0 ,
134  dc5 = -277.0/14336.0 ;
135 
136  const G4double dc1 = c1 - 2825.0/27648.0 , dc3 = c3 - 18575.0/48384.0 ,
137  dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
138 
139  // Initialise time to t0, needed when it is not updated by the integration.
140  // [ Note: Only for time dependent fields (usually electric)
141  // is it neccessary to integrate the time.]
142  yOut[7] = yTemp[7] = yIn[7];
143 
144  const G4int numberOfVariables= this->GetNumberOfVariables();
145  // The number of variables to be integrated over
146 
147  // Saving yInput because yInput and yOut can be aliases for same array
148 
149  for(i=0;i<numberOfVariables;i++)
150  {
151  yIn[i]=yInput[i];
152  }
153  // RightHandSide(yIn, dydx) ; // 1st Step
154 
155  for(i=0;i<numberOfVariables;i++)
156  {
157  yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
158  }
159  RightHandSide(yTemp, ak2) ; // 2nd Step
160 
161  for(i=0;i<numberOfVariables;i++)
162  {
163  yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
164  }
165  RightHandSide(yTemp, ak3) ; // 3rd Step
166 
167  for(i=0;i<numberOfVariables;i++)
168  {
169  yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
170  }
171  RightHandSide(yTemp, ak4) ; // 4th Step
172 
173  for(i=0;i<numberOfVariables;i++)
174  {
175  yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
176  b54*ak4[i]) ;
177  }
178  RightHandSide(yTemp, ak5) ; // 5th Step
179 
180  for(i=0;i<numberOfVariables;i++)
181  {
182  yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
183  b64*ak4[i] + b65*ak5[i]) ;
184  }
185  RightHandSide(yTemp, ak6) ; // 6th Step
186 
187  for(i=0;i<numberOfVariables;i++)
188  {
189  // Accumulate increments with proper weights
190 
191  yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
192 
193  // Estimate error as difference between 4th and
194  // 5th order methods
195 
196  yErr[i] = Step*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] +
197  dc5*ak5[i] + dc6*ak6[i]) ;
198 
199  // Store Input and Final values, for possible use in calculating chord
200  fLastInitialVector[i] = yIn[i] ;
201  fLastFinalVector[i] = yOut[i];
202  fLastDyDx[i] = dydx[i];
203  }
204  // NormaliseTangentVector( yOut ); // Not wanted
205 
206  fLastStepLength =Step;
207 
208  return ;
209 }
210 
212 
213 void
214 G4CashKarpRKF45::StepWithEst( const G4double*,
215  const G4double*,
216  G4double,
217  G4double*,
218  G4double&,
219  G4double&,
220  const G4double*,
221  G4double* )
222 {
223  G4Exception("G4CashKarpRKF45::StepWithEst()", "GeomField0001",
224  FatalException, "Method no longer used.");
225  return ;
226 }
227 
229 
231 {
232  G4double distLine, distChord;
233  G4ThreeVector initialPoint, finalPoint, midPoint;
234 
235  // Store last initial and final points (they will be overwritten in self-Stepper call!)
236  initialPoint = G4ThreeVector( fLastInitialVector[0],
237  fLastInitialVector[1], fLastInitialVector[2]);
238  finalPoint = G4ThreeVector( fLastFinalVector[0],
239  fLastFinalVector[1], fLastFinalVector[2]);
240 
241  // Do half a step using StepNoErr
242 
243  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
244  fMidVector, fMidError );
245 
246  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
247 
248  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
249  // distance of Chord
250 
251 
252  if (initialPoint != finalPoint)
253  {
254  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
255  distChord = distLine;
256  }
257  else
258  {
259  distChord = (midPoint-initialPoint).mag();
260  }
261  return distChord;
262 }
263 
264 
CLHEP::Hep3Vector G4ThreeVector
G4CashKarpRKF45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
TCanvas * c1
Definition: plotHisto.C:7
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfVariables() const
bool G4bool
Definition: G4Types.hh:79
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double DistChord() const
Definition: Step.hh:41
void RightHandSide(const double y[], double dydx[])
double G4double
Definition: G4Types.hh:76