| 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 // G4TsitourasRK45 implementation 26 // G4TsitourasRK45 implementation
27 // 27 //
28 // Author: Somnath Banerjee, Google Summer of 28 // Author: Somnath Banerjee, Google Summer of Code 2015, 11.06.2015
29 // Supervision: John Apostolakis, CERN 29 // Supervision: John Apostolakis, CERN
30 // ------------------------------------------- 30 // -------------------------------------------------------------------
31 31
32 #include "G4TsitourasRK45.hh" 32 #include "G4TsitourasRK45.hh"
33 #include "G4LineSection.hh" 33 #include "G4LineSection.hh"
34 34
35 ////////////////////////////////////////////// 35 /////////////////////////////////////////////////////////////////////
36 // 36 //
37 // Constructor 37 // Constructor
38 // 38 //
39 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfM 39 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfMotion *EqRhs,
40 G4int noInteg 40 G4int noIntegrationVariables,
41 G4bool primar 41 G4bool primary)
42 : G4MagIntegratorStepper(EqRhs, noIntegratio 42 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
43 { 43 {
44 const G4int numberOfVariables = noIntegratio 44 const G4int numberOfVariables = noIntegrationVariables;
45 45
46 ak2 = new G4double[numberOfVariables] ; 46 ak2 = new G4double[numberOfVariables] ;
47 ak3 = new G4double[numberOfVariables] ; 47 ak3 = new G4double[numberOfVariables] ;
48 ak4 = new G4double[numberOfVariables] ; 48 ak4 = new G4double[numberOfVariables] ;
49 ak5 = new G4double[numberOfVariables] ; 49 ak5 = new G4double[numberOfVariables] ;
50 ak6 = new G4double[numberOfVariables] ; 50 ak6 = new G4double[numberOfVariables] ;
51 ak7 = new G4double[numberOfVariables] ; 51 ak7 = new G4double[numberOfVariables] ;
52 ak8 = new G4double[numberOfVariables] ; 52 ak8 = new G4double[numberOfVariables] ;
53 53
54 // Must ensure space extra 'state' variables 54 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
55 // 55 //
56 const G4int numStateMax = std::max(GetNumbe 56 const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8);
57 const G4int numStateVars = std::max(noIntegr 57 const G4int numStateVars = std::max(noIntegrationVariables,
58 numState 58 numStateMax );
59 59
60 yTemp = new G4double[numStateVars] ; 60 yTemp = new G4double[numStateVars] ;
61 yIn = new G4double[numStateVars] ; 61 yIn = new G4double[numStateVars] ;
62 62
63 fLastInitialVector = new G4double[numberOfVa 63 fLastInitialVector = new G4double[numberOfVariables] ;
64 fLastFinalVector = new G4double[numberOfVari 64 fLastFinalVector = new G4double[numberOfVariables] ;
65 65
66 fLastDyDx = new G4double[numberOfVariables]; 66 fLastDyDx = new G4double[numberOfVariables];
67 67
68 fMidVector = new G4double[numberOfVariables] 68 fMidVector = new G4double[numberOfVariables];
69 fMidError = new G4double[numberOfVariables] 69 fMidError = new G4double[numberOfVariables];
70 70
71 if( primary ) 71 if( primary )
72 { 72 {
73 fAuxStepper = new G4TsitourasRK45(EqRhs, n 73 fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
74 } 74 }
75 } 75 }
76 76
77 ////////////////////////////////////////////// 77 /////////////////////////////////////////////////////////////////////
78 // 78 //
79 // Destructor 79 // Destructor
80 // 80 //
81 G4TsitourasRK45::~G4TsitourasRK45() 81 G4TsitourasRK45::~G4TsitourasRK45()
82 { 82 {
83 delete [] ak2; 83 delete [] ak2;
84 delete [] ak3; 84 delete [] ak3;
85 delete [] ak4; 85 delete [] ak4;
86 delete [] ak5; 86 delete [] ak5;
87 delete [] ak6; 87 delete [] ak6;
88 delete [] ak7; 88 delete [] ak7;
89 delete [] ak8; 89 delete [] ak8;
90 90
91 delete [] yTemp; 91 delete [] yTemp;
92 delete [] yIn; 92 delete [] yIn;
93 93
94 delete [] fLastInitialVector; 94 delete [] fLastInitialVector;
95 delete [] fLastFinalVector; 95 delete [] fLastFinalVector;
96 delete [] fLastDyDx; 96 delete [] fLastDyDx;
97 delete [] fMidVector; 97 delete [] fMidVector;
98 delete [] fMidError; 98 delete [] fMidError;
99 99
100 delete fAuxStepper; 100 delete fAuxStepper;
101 } 101 }
102 102
103 // The following coefficients have been obtain 103 // The following coefficients have been obtained from
104 // Table 1: The Coefficients of the new pair 104 // Table 1: The Coefficients of the new pair
105 // 105 //
106 // C. Tsitouras, "Runge–Kutta pairs of order 106 // C. Tsitouras, "Runge–Kutta pairs of order 5(4) satisfying only
107 // the first column simplifying 107 // the first column simplifying assumption"
108 // Computers & Mathematics with Applications, 108 // Computers & Mathematics with Applications, vol.62, no.2, pp.770-775, 2011.
109 // 109 //
110 // A corresponding matlab code was also found 110 // A corresponding matlab code was also found at:
111 // http://users.ntua.gr/tsitoura/new54.m 111 // http://users.ntua.gr/tsitoura/new54.m
112 // 112 //
113 // Doing a step 113 // Doing a step
114 // 114 //
115 void 115 void
116 G4TsitourasRK45::Stepper( const G4double yInpu 116 G4TsitourasRK45::Stepper( const G4double yInput[],
117 const G4double dydx[ 117 const G4double dydx[],
118 G4double Step, 118 G4double Step,
119 G4double yOut[ 119 G4double yOut[],
120 G4double yErr[ 120 G4double yErr[] )
121 { 121 {
122 const G4double b21 = 0.161 , 122 const G4double b21 = 0.161 ,
123 b31 = -0.008480655492356988 123 b31 = -0.00848065549235698854 ,
124 b32 = 0.335480655492356989 124 b32 = 0.335480655492356989 ,
125 125
126 b41 = 2.89715305710549343 126 b41 = 2.89715305710549343 ,
127 b42 = -6.35944848997507484 127 b42 = -6.35944848997507484 ,
128 b43 = 4.36229543286958141 , 128 b43 = 4.36229543286958141 ,
129 129
130 b51 = 5.325864828439257, 130 b51 = 5.325864828439257,
131 b52 = -11.748883564062828, 131 b52 = -11.748883564062828,
132 b53 = 7.49553934288983621 , 132 b53 = 7.49553934288983621 ,
133 b54 = -0.09249506636175525, 133 b54 = -0.09249506636175525,
134 134
135 b61 = 5.8614554429464200, 135 b61 = 5.8614554429464200,
136 b62 = -12.9209693178471093 136 b62 = -12.9209693178471093 ,
137 b63 = 8.1593678985761586 , 137 b63 = 8.1593678985761586 ,
138 b64 = -0.071584973281400997 138 b64 = -0.071584973281400997,
139 b65 = -0.028269050394068382 139 b65 = -0.0282690503940683829,
140 140
141 b71 = 0.0964607668180652295 141 b71 = 0.0964607668180652295 ,
142 b72 = 0.01, 142 b72 = 0.01,
143 b73 = 0.479889650414499575, 143 b73 = 0.479889650414499575,
144 b74 = 1.37900857410374189, 144 b74 = 1.37900857410374189,
145 b75 = -3.2900695154360807, 145 b75 = -3.2900695154360807,
146 b76 = 2.32471052409977398, 146 b76 = 2.32471052409977398,
147 147
148 // c1 = 0.001780011052226 , 148 // c1 = 0.001780011052226 ,
149 // c2 = 0.000816434459657 , 149 // c2 = 0.000816434459657 ,
150 // c3 = -0.007880878010262 , 150 // c3 = -0.007880878010262 ,
151 // c4 = 0.144711007173263 , 151 // c4 = 0.144711007173263 ,
152 // c5 = -0.582357165452555 , 152 // c5 = -0.582357165452555 ,
153 // c6 = 0.458082105929187 , 153 // c6 = 0.458082105929187 ,
154 // c7 = 1.0/66.0 ; 154 // c7 = 1.0/66.0 ;
155 155
156 dc1 = 0.0935237485818927066 156 dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
157 dc2 = 0.0086528831415663676 157 dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
158 dc3 = 0.492893099131431868 158 dc3 = 0.492893099131431868 - b73 , // + 0.007880878010262,
159 dc4 = 1.14023541226785810 - 159 dc4 = 1.14023541226785810 - b74 , // 0.144711007173263,
160 dc5 = - 2.3291801924393646 160 dc5 = - 2.3291801924393646 - b75, // + 0.582357165452555,
161 dc6 = 1.56887504931661552 - 161 dc6 = 1.56887504931661552 - b76 , // - 0.458082105929187,
162 dc7 = 0.025; //- 1.0/66.0 ; 162 dc7 = 0.025; //- 1.0/66.0 ;
163 163
164 // dc1 = -3.0/1280.0, 164 // dc1 = -3.0/1280.0,
165 // dc2 = 0.0, 165 // dc2 = 0.0,
166 // dc3 = 6561.0/632320.0, 166 // dc3 = 6561.0/632320.0,
167 // dc4 = -343.0/20800.0, 167 // dc4 = -343.0/20800.0,
168 // dc5 = 243.0/12800.0, 168 // dc5 = 243.0/12800.0,
169 // dc6 = -1.0/95.0, 169 // dc6 = -1.0/95.0,
170 // dc7 = 0.0 ; 170 // dc7 = 0.0 ;
171 171
172 const G4int numberOfVariables = GetNumberO 172 const G4int numberOfVariables = GetNumberOfVariables();
173 173
174 // The number of variables to be integrate 174 // The number of variables to be integrated over
175 // 175 //
176 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; 176 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
177 177
178 // Saving yInput because yInput and yOut 178 // Saving yInput because yInput and yOut can be aliases for same array
179 // 179 //
180 for(G4int i=0; i<numberOfVariables; ++i) 180 for(G4int i=0; i<numberOfVariables; ++i)
181 { 181 {
182 yIn[i]=yInput[i]; 182 yIn[i]=yInput[i];
183 } 183 }
184 // RightHandSide(yIn, dydx) ; // 1st Ste 184 // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
185 185
186 for(G4int i=0; i<numberOfVariables; ++i) 186 for(G4int i=0; i<numberOfVariables; ++i)
187 { 187 {
188 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; 188 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
189 } 189 }
190 RightHandSide(yTemp, ak2) ; / 190 RightHandSide(yTemp, ak2) ; // 2nd Stage
191 191
192 for(G4int i=0; i<numberOfVariables; ++i) 192 for(G4int i=0; i<numberOfVariables; ++i)
193 { 193 {
194 yTemp[i] = yIn[i] + Step*(b31*dydx[i] 194 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
195 } 195 }
196 RightHandSide(yTemp, ak3) ; / 196 RightHandSide(yTemp, ak3) ; // 3rd Stage
197 197
198 for(G4int i=0; i<numberOfVariables; ++i) 198 for(G4int i=0; i<numberOfVariables; ++i)
199 { 199 {
200 yTemp[i] = yIn[i] + Step*(b41*dydx[i] 200 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
201 } 201 }
202 RightHandSide(yTemp, ak4) ; / 202 RightHandSide(yTemp, ak4) ; // 4th Stage
203 203
204 for(G4int i=0; i<numberOfVariables; ++i) 204 for(G4int i=0; i<numberOfVariables; ++i)
205 { 205 {
206 yTemp[i] = yIn[i] + Step*(b51*dydx[i] 206 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
207 b54*ak4[i]) 207 b54*ak4[i]) ;
208 } 208 }
209 RightHandSide(yTemp, ak5) ; / 209 RightHandSide(yTemp, ak5) ; // 5th Stage
210 210
211 for(G4int i=0; i<numberOfVariables; ++i) 211 for(G4int i=0; i<numberOfVariables; ++i)
212 { 212 {
213 yTemp[i] = yIn[i] + Step*(b61*dydx[i] 213 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
214 b64*ak4[i] + 214 b64*ak4[i] + b65*ak5[i]) ;
215 } 215 }
216 RightHandSide(yTemp, ak6) ; / 216 RightHandSide(yTemp, ak6) ; // 6th Stage
217 217
218 for(G4int i=0; i<numberOfVariables; ++i) 218 for(G4int i=0; i<numberOfVariables; ++i)
219 { 219 {
220 yOut[i] = yIn[i] + Step*(b71*dydx[i] + 220 yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
221 b74*ak4[i] + 221 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
222 } 222 }
223 RightHandSide(yOut, ak7); / 223 RightHandSide(yOut, ak7); // 7th Stage
224 224
225 //Calculate the error in the step: 225 //Calculate the error in the step:
226 for(G4int i=0; i<numberOfVariables; ++i) 226 for(G4int i=0; i<numberOfVariables; ++i)
227 { 227 {
228 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[ 228 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
229 dc5*ak5[i] + dc6*ak6[ 229 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
230 230
231 // Store Input and Final values, for p 231 // Store Input and Final values, for possible use in calculating chord
232 // 232 //
233 fLastInitialVector[i] = yIn[i] ; 233 fLastInitialVector[i] = yIn[i] ;
234 fLastFinalVector[i] = yOut[i]; 234 fLastFinalVector[i] = yOut[i];
235 fLastDyDx[i] = dydx[i]; 235 fLastDyDx[i] = dydx[i];
236 } 236 }
237 237
238 fLastStepLength = Step; 238 fLastStepLength = Step;
239 239
240 return ; 240 return ;
241 } 241 }
242 242
243 void G4TsitourasRK45::SetupInterpolation() 243 void G4TsitourasRK45::SetupInterpolation()
244 // (const G4double *yInput, const G4double * 244 // (const G4double *yInput, const G4double *dydx, const G4double Step)
245 { 245 {
246 // Nothing to be done 246 // Nothing to be done
247 } 247 }
248 248
249 void G4TsitourasRK45::Interpolate(const G4doub 249 void G4TsitourasRK45::Interpolate(const G4double* yInput,
250 const G4doub 250 const G4double* dydx,
251 const G4doub 251 const G4double Step,
252 G4doub 252 G4double* yOut,
253 G4doub 253 G4double tau)
254 { 254 {
255 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7 255 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
256 // Coefficients for all the seven stages 256 // Coefficients for all the seven stages.
257 257
258 const G4int numberOfVariables = GetNumberO 258 const G4int numberOfVariables = GetNumberOfVariables();
259 259
260 G4double tau0 = tau; 260 G4double tau0 = tau;
261 261
262 for(G4int i=0; i<numberOfVariables; ++i) 262 for(G4int i=0; i<numberOfVariables; ++i)
263 { 263 {
264 yIn[i] = yInput[i]; 264 yIn[i] = yInput[i];
265 } 265 }
266 266
267 G4double tau_2 = tau0*tau0 ; 267 G4double tau_2 = tau0*tau0 ;
268 // tau_3 = tau0*tau_2, 268 // tau_3 = tau0*tau_2,
269 // tau_4 = tau_2*tau_2; 269 // tau_4 = tau_2*tau_2;
270 270
271 bf1 = -1.0530884977290216*tau*(tau - 1.329 271 bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
272 1.4364028541716351*tau + 0.713 272 1.4364028541716351*tau + 0.7139816917074209);
273 bf2 = 0.1017*tau_2*(tau_2 - 2.196656833824 273 bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
274 1.2949852507374631); 274 1.2949852507374631);
275 bf3 = 2.490627285651252793*tau_2*(tau_2 - 275 bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
276 + 1.5780 276 + 1.57803468208092486);
277 bf4 = -16.54810288924490272*(tau - 1.21712 277 bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
278 (tau - 0.6 278 (tau - 0.61620406037800089)*tau_2;
279 bf5 = 47.37952196281928122*(tau - 1.203071 279 bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
280 (tau - 0.6 280 (tau - 0.658047292653547382)*tau_2;
281 bf6 = -34.87065786149660974*(tau - 1.2)*(t 281 bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
282 0.666666666666 282 0.666666666666666667)*tau_2;
283 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2; 283 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
284 284
285 // Putting together the coefficients calcu 285 // Putting together the coefficients calculated as the respective
286 // stage coefficients 286 // stage coefficients
287 // 287 //
288 for(G4int i=0; i<numberOfVariables; ++i) 288 for(G4int i=0; i<numberOfVariables; ++i)
289 { 289 {
290 yOut[i] = yIn[i] + Step*( bf1*dydx[i] 290 yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
291 + bf4*ak4[i] 291 + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
292 + bf7*ak7[i] 292 + bf7*ak7[i] ) ;
293 } 293 }
294 } 294 }
295 295
296 G4double G4TsitourasRK45::DistChord() const 296 G4double G4TsitourasRK45::DistChord() const
297 { 297 {
298 G4double distLine, distChord; 298 G4double distLine, distChord;
299 G4ThreeVector initialPoint, finalPoint, midP 299 G4ThreeVector initialPoint, finalPoint, midPoint;
300 300
301 // Store last initial and final points (they 301 // Store last initial and final points (they will be
302 // overwritten in self-Stepper call!) 302 // overwritten in self-Stepper call!)
303 // 303 //
304 initialPoint = G4ThreeVector( fLastInitialVe 304 initialPoint = G4ThreeVector( fLastInitialVector[0],
305 fLastInitialVe 305 fLastInitialVector[1], fLastInitialVector[2]);
306 finalPoint = G4ThreeVector( fLastFinalVect 306 finalPoint = G4ThreeVector( fLastFinalVector[0],
307 fLastFinalVect 307 fLastFinalVector[1], fLastFinalVector[2]);
308 308
309 // Do half a step using StepNoErr 309 // Do half a step using StepNoErr
310 // 310 //
311 fAuxStepper->Stepper( fLastInitialVector, fL 311 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
312 fMidVector, fMidError ); 312 fMidVector, fMidError );
313 313
314 midPoint = G4ThreeVector( fMidVector[0], fMi 314 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
315 315
316 // Use stored values of Initial and Endpoint 316 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
317 // distance of Chord 317 // distance of Chord
318 // 318 //
319 if (initialPoint != finalPoint) 319 if (initialPoint != finalPoint)
320 { 320 {
321 distLine = G4LineSection::Distline( midP 321 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
322 distChord = distLine; 322 distChord = distLine;
323 } 323 }
324 else 324 else
325 { 325 {
326 distChord = (midPoint-initialPoint).mag() 326 distChord = (midPoint-initialPoint).mag();
327 } 327 }
328 return distChord; 328 return distChord;
329 } 329 }
330 330