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