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