source: trunk/source/geometry/magneticfield/src/G4CashKarpRKF45.cc@ 1036

Last change on this file since 1036 was 921, checked in by garnier, 17 years ago

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