| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4DormandPrinceRK78 implementation
27 //
28 // Dormand-Prince 8(7)13M non-FSAL, based on R
29 // P.J. Prince, J.R. Dormand, "High order embe
30 // Journal of Computational and Applied Mathem
31 // Pages 67-75, ISSN 0377-0427, DOI: 10.1016/0
32 //
33 // Created: Somnath Banerjee, Google Summer of
34 // Supervision: John Apostolakis, CERN
35 // -------------------------------------------
36
37 #include "G4DormandPrinceRK78.hh"
38 #include "G4LineSection.hh"
39
40 // Constructor
41 //
42 G4DormandPrinceRK78::G4DormandPrinceRK78(G4Equ
43 G4int
44 G4boo
45 : G4MagIntegratorStepper(EqRhs, noIntegrati
46 {
47 const G4int numberOfVariables = noIntegrat
48
49 // New Chunk of memory being created for u
50
51 // aki - for storing intermediate RHS
52 //
53 ak2 = new G4double[numberOfVariables];
54 ak3 = new G4double[numberOfVariables];
55 ak4 = new G4double[numberOfVariables];
56 ak5 = new G4double[numberOfVariables];
57 ak6 = new G4double[numberOfVariables];
58 ak7 = new G4double[numberOfVariables];
59 ak8 = new G4double[numberOfVariables];
60 ak9 = new G4double[numberOfVariables];
61 ak10 = new G4double[numberOfVariables];
62 ak11 = new G4double[numberOfVariables];
63 ak12 = new G4double[numberOfVariables];
64 ak13 = new G4double[numberOfVariables];
65
66 const G4int numStateVars = std::max(noInte
67 yTemp = new G4double[numStateVars];
68 yIn = new G4double[numStateVars] ;
69
70 fLastInitialVector = new G4double[numState
71 fLastFinalVector = new G4double[numStateVa
72
73 fLastDyDx = new G4double[numStateVars];
74
75 fMidVector = new G4double[numStateVars];
76 fMidError = new G4double[numStateVars];
77
78 if( primary )
79 {
80 fAuxStepper = new G4DormandPrinceRK78(Eq
81 }
82 }
83
84 // Destructor
85 //
86 G4DormandPrinceRK78::~G4DormandPrinceRK78()
87 {
88 // Clear all previously allocated memory f
89
90 delete [] ak2;
91 delete [] ak3;
92 delete [] ak4;
93 delete [] ak5;
94 delete [] ak6;
95 delete [] ak7;
96 delete [] ak8;
97 delete [] ak9;
98 delete [] ak10;
99 delete [] ak11;
100 delete [] ak12;
101 delete [] ak13;
102 delete [] yTemp;
103 delete [] yIn;
104
105 delete [] fLastInitialVector;
106 delete [] fLastFinalVector;
107 delete [] fLastDyDx;
108 delete [] fMidVector;
109 delete [] fMidError;
110
111 delete fAuxStepper;
112 }
113
114
115 // The following scheme and the set of coeffic
116 // Table2. RK8(7)13M (Rational approximations)
117 // P. J. Prince and J. R. Dormand, "High order
118 // Journal of Computational and Applied Math.,
119 //
120 // Stepper :
121 //
122 // Passing in the value of yInput[],the first
123 // Giving back yOut and yErr arrays for output
124 //
125 void G4DormandPrinceRK78::Stepper(const G4doub
126 const G4doub
127 G4doub
128 G4doub
129 G4doub
130 {
131 G4int i;
132
133 // The various constants defined on the ba
134 //
135 const G4double b21 = 1.0/18,
136 b31 = 1.0/48.0 ,
137 b32 = 1.0/16.0 ,
138
139 b41 = 1.0/32.0 ,
140 b42 = 0.0 ,
141 b43 = 3.0/32.0 ,
142
143 b51 = 5.0/16.0 ,
144 b52 = 0.0 ,
145 b53 = -75.0/64.0 ,
146 b54 = 75.0/64.0 ,
147
148 b61 = 3.0/80.0 ,
149 b62 = 0.0 ,
150 b63 = 0.0 ,
151 b64 = 3.0/16.0 ,
152 b65 = 3.0/20.0 ,
153
154 b71 = 29443841.0/614563906.
155 b72 = 0.0 ,
156 b73 = 0.0 ,
157 b74 = 77736538.0/692538347.
158 b75 = -28693883.0/112500000
159 b76 = 23124283.0/1800000000
160
161 b81 = 16016141.0/946692911.
162 b82 = 0.0 ,
163 b83 = 0.0 ,
164 b84 = 61564180.0/158732637.
165 b85 = 22789713.0/633445777.
166 b86 = 545815736.0/277105722
167 b87 = -180193667.0/10433075
168
169 b91 = 39632708.0/573591083.
170 b92 = 0.0 ,
171 b93 = 0.0 ,
172 b94 = -433636366.0/68370161
173 b95 = -421739975.0/26162923
174 b96 = 100302831.0/723423059
175 b97 = 790204164.0/839813087
176 b98 = 800635310.0/378307128
177
178 b101 = 246121993.0/13408477
179 b102 = 0.0 ,
180 b103 = 0.0 ,
181 b104 = -37695042795.0/15268
182 b105 = -309121744.0/1061227
183 b106 = -12992083.0/4907669
184 b107 = 6005943493.0/2108947
185 b108 = 393006217.0/13966734
186 b109 = 123872331.0/10010297
187
188 b111 = -1028468189.0/846180
189 b112 = 0.0 ,
190 b113 = 0.0 ,
191 b114 = 8478235783.0/5085128
192 b115 = 1311729495.0/1432422
193 b116 = -10304129995.0/17013
194 b117 = -48777925059.0/3047
195 b118 = 15336726248.0/103282
196 b119 = -45442868181.0/3398
197 b1110 = 3065993473.0/597172
198
199 b121 = 185892177.0/71811604
200 b122 = 0.0 ,
201 b123 = 0.0 ,
202 b124 = -3185094517.0/667107
203 b125 = -477755414.0/1098053
204 b126 = -703635378.0/2307392
205 b127 = 5731566787.0/102754
206 b128 = 5232866602.0/8500665
207 b129 = -4093664535.0/808688
208 b1210 = 3962137247.0/180595
209 b1211 = 65686358.0/48791008
210
211 b131 = 403863854.0/49106310
212 b132 = 0.0 ,
213 b133 = 0.0 ,
214 b134 = -5068492393.0/434740
215 b135 = -411421997.0/5430438
216 b136 = 652783627.0/91429660
217 b137 = 11173962825.0/925320
218 b138 = -13158990841.0/61847
219 b139 = 3936647629.0/1978049
220 b1310 = -160528059.0/685178
221 b1311 = 248638103.0/1413531
222 b1312 = 0.0 ,
223
224 c1 = 14005451.0/335480064.0
225 // c2 = 0.0 ,
226 // c3 = 0.0 ,
227 // c4 = 0.0 ,
228 // c5 = 0.0 ,
229 c6 = -59238493.0/1068277825
230 c7 = 181606767.0/758867731.
231 c8 = 561292985.0/797845732.
232 c9 = -1041891430.0/1371343
233 c10 = 760417239.0/115116529
234 c11 = 118820643.0/751138087
235 c12 = - 528747749.0/2220607
236 c13 = 1.0/4.0 ,
237
238 c_1 = 13451932.0/455176623.
239 // c_2 = 0.0 ,
240 // c_3 = 0.0 ,
241 // c_4 = 0.0 ,
242 // c_5 = 0.0 ,
243 c_6 = -808719846.0/97600014
244 c_7 = 1757004468.0/56451593
245 c_8 = 656045339.0/265891186
246 c_9 = -3867574721.0/1518517
247 c_10 = 465885868.0/32273653
248 c_11 = 53011238.0/667516719
249 c_12 = 2.0/45.0 ,
250 c_13 = 0.0 ,
251
252 dc1 = c_1 - c1 ,
253 // dc2 = c_2 - c2 ,
254 // dc3 = c_3 - c3 ,
255 // dc4 = c_4 - c4 ,
256 // dc5 = c_5 - c5 ,
257 dc6 = c_6 - c6 ,
258 dc7 = c_7 - c7 ,
259 dc8 = c_8 - c8 ,
260 dc9 = c_9 - c9 ,
261 dc10 = c_10 - c10 ,
262 dc11 = c_11 - c11 ,
263 dc12 = c_12 - c12 ,
264 dc13 = c_13 - c13 ;
265 //
266 // end of declaration !
267
268 const G4int numberOfVariables = GetNumberO
269
270 // The number of variables to be integrate
271 //
272 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
273
274 // Saving yInput because yInput and yOut c
275 //
276 for(i=0; i<numberOfVariables; ++i)
277 {
278 yIn[i]=yInput[i];
279 }
280 // RightHandSide(yIn, dydx) ; // 1st Sta
281
282 for(i=0; i<numberOfVariables; ++i)
283 {
284 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
285 }
286 RightHandSide(yTemp, ak2) ; /
287
288 for(i=0; i<numberOfVariables; ++i)
289 {
290 yTemp[i] = yIn[i] + Step*(b31*dydx[i]
291 }
292 RightHandSide(yTemp, ak3) ; /
293
294 for(i=0; i<numberOfVariables; ++i)
295 {
296 yTemp[i] = yIn[i] + Step*(b41*dydx[i]
297 }
298 RightHandSide(yTemp, ak4) ; /
299
300 for(i=0; i<numberOfVariables; ++i)
301 {
302 yTemp[i] = yIn[i] + Step*(b51*dydx[i]
303 b54*ak4[i])
304 }
305 RightHandSide(yTemp, ak5) ; /
306
307 for(i=0; i<numberOfVariables; ++i)
308 {
309 yTemp[i] = yIn[i] + Step*(b61*dydx[i]
310 b64*ak4[i] +
311 }
312 RightHandSide(yTemp, ak6) ; /
313
314 for(i=0; i<numberOfVariables; ++i)
315 {
316 yTemp[i] = yIn[i] + Step*(b71*dydx[i]
317 b74*ak4[i] +
318 }
319 RightHandSide(yTemp, ak7); /
320
321 for(i=0; i<numberOfVariables; ++i)
322 {
323 yTemp[i] = yIn[i] + Step*(b81*dydx[i]
324 b84*ak4[i] +
325 b87*ak7[i]);
326 }
327 RightHandSide(yTemp, ak8); /
328
329 for(i=0; i<numberOfVariables; ++i)
330 {
331 yTemp[i] = yIn[i] + Step*(b91*dydx[i]
332 b94*ak4[i] +
333 b97*ak7[i] +
334 }
335 RightHandSide(yTemp, ak9); /
336
337 for(i=0; i<numberOfVariables; ++i)
338 {
339 yTemp[i] = yIn[i] + Step*(b101*dydx[i]
340 b104*ak4[i]
341 b107*ak7[i]
342 }
343 RightHandSide(yTemp, ak10); /
344
345 for(i=0; i<numberOfVariables; ++i)
346 {
347 yTemp[i] = yIn[i] + Step*(b111*dydx[i]
348 b114*ak4[i]
349 b117*ak7[i]
350 b1110*ak10[i
351 }
352 RightHandSide(yTemp, ak11); /
353
354 for(i=0; i<numberOfVariables; ++i)
355 {
356 yTemp[i] = yIn[i] + Step*(b121*dydx[i]
357 b124*ak4[i]
358 b127*ak7[i]
359 b1210*ak10[i
360 }
361 RightHandSide(yTemp, ak12); /
362
363 for(i=0; i<numberOfVariables; ++i)
364 {
365 yTemp[i] = yIn[i]+Step*(b131*dydx[i] +
366 b134*ak4[i] +
367 b137*ak7[i] +
368 b1310*ak10[i]
369 }
370 RightHandSide(yTemp, ak13); /
371
372 for(i=0; i<numberOfVariables; ++i)
373 {
374 // Accumulate increments with proper w
375
376 yOut[i] = yIn[i] + Step*(c1*dydx[i] +
377 // + c4 * ak4
378 + c6*ak6[i] +
379 c7*ak7[i] + c
380 + c11*ak11[i]
381
382 // Estimate error as difference betwee
383 //
384 yErr[i] = Step*(dc1*dydx[i] + // dc2*a
385 // + dc5*ak5[i]
386 + dc6*ak6[i] + dc7*ak7[i
387 + dc9*ak9[i] + dc10*ak10
388 + dc13*ak13[i] ) ;
389
390 // Store Input and Final values, for p
391 //
392 fLastInitialVector[i] = yIn[i] ;
393 fLastFinalVector[i] = yOut[i];
394 fLastDyDx[i] = dydx[i];
395 }
396 fLastStepLength = Step;
397
398 return ;
399 }
400
401 // DistChord
402 //
403 G4double G4DormandPrinceRK78::DistChord() con
404 {
405 G4double distLine, distChord;
406 G4ThreeVector initialPoint, finalPoint, mi
407
408 // Store last initial and final points
409 // (they will be overwritten in self-Stepp
410 //
411 initialPoint = G4ThreeVector( fLastInitial
412 fLastInitialV
413 finalPoint = G4ThreeVector( fLastFinalVe
414 fLastFinalVec
415
416 // Do half a step using StepNoErr
417
418 fAuxStepper->Stepper( fLastInitialVector,
419 fMidVector, fMidErr
420
421 midPoint = G4ThreeVector( fMidVector[0], f
422
423 // Use stored values of Initial and Endpoi
424 // distance of Chord
425 //
426 if (initialPoint != finalPoint)
427 {
428 distLine = G4LineSection::Distline(mid
429 distChord = distLine;
430 }
431 else
432 {
433 distChord = (midPoint-initialPoint).ma
434 }
435 return distChord;
436 }
437