| 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 // G4DormandPrinceRK56 implementation << 26 // FDormand-Prince RK 6(5) FSAL implementation by Somnath Banerjee
>> 27 // Supervision / code review: John Apostolakis
27 // 28 //
28 // Created: Somnath Banerjee, Google Summer of << 29 // Sponsored by Google in Google Summer of Code 2015.
29 // Supervision: John Apostolakis, CERN << 30 //
30 // ------------------------------------------- << 31 // First version: 26 June 2015
>> 32 //
>> 33 // G4DormandPrince745.cc
>> 34 // Geant4
>> 35 //
>> 36 // History
>> 37 // -----------------------------
>> 38 // Created by Somnath on 26 June 2015
>> 39 //
>> 40 //
>> 41 ///////////////////////////////////////////////////////////////////////////////
31 42
32 #include "G4DormandPrinceRK56.hh" 43 #include "G4DormandPrinceRK56.hh"
33 #include "G4LineSection.hh" 44 #include "G4LineSection.hh"
34 45
35 // Constructor << 46 //Constructor
36 // << 47 G4DormandPrinceRK56::G4DormandPrinceRK56(G4EquationOfMotion *EqRhs,
37 G4DormandPrinceRK56::G4DormandPrinceRK56(G4Equ << 48 G4int noIntegrationVariables,
38 G4int << 49 G4bool primary)
39 G4boo << 50 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
40 : G4MagIntegratorStepper(EqRhs, noIntegratio << 51 fLastStepLength(-1.0), fAuxStepper(nullptr)
41 { 52 {
42 const G4int numberOfVariables = noIntegrat 53 const G4int numberOfVariables = noIntegrationVariables;
43 54
44 // New Chunk of memory being created for u << 55 //New Chunk of memory being created for use by the stepper
45 56
46 // aki - for storing intermediate RHS << 57 //aki - for storing intermediate RHS
47 // <<
48 ak2 = new G4double[numberOfVariables]; 58 ak2 = new G4double[numberOfVariables];
49 ak3 = new G4double[numberOfVariables]; 59 ak3 = new G4double[numberOfVariables];
50 ak4 = new G4double[numberOfVariables]; 60 ak4 = new G4double[numberOfVariables];
51 ak5 = new G4double[numberOfVariables]; 61 ak5 = new G4double[numberOfVariables];
52 ak6 = new G4double[numberOfVariables]; 62 ak6 = new G4double[numberOfVariables];
53 ak7 = new G4double[numberOfVariables]; 63 ak7 = new G4double[numberOfVariables];
54 ak8 = new G4double[numberOfVariables]; 64 ak8 = new G4double[numberOfVariables];
55 ak9 = new G4double[numberOfVariables]; 65 ak9 = new G4double[numberOfVariables];
56 66
57 // Memory for Additional stages 67 // Memory for Additional stages
58 // <<
59 ak10 = new G4double[numberOfVariables]; 68 ak10 = new G4double[numberOfVariables];
60 ak11 = new G4double[numberOfVariables]; 69 ak11 = new G4double[numberOfVariables];
61 ak12 = new G4double[numberOfVariables]; 70 ak12 = new G4double[numberOfVariables];
62 ak10_low = new G4double[numberOfVariables] 71 ak10_low = new G4double[numberOfVariables];
63 72
64 const G4int numStateVars = std::max(noInte 73 const G4int numStateVars = std::max(noIntegrationVariables, 8);
65 yTemp = new G4double[numStateVars]; 74 yTemp = new G4double[numStateVars];
66 yIn = new G4double[numStateVars] ; 75 yIn = new G4double[numStateVars] ;
67 76
68 fLastInitialVector = new G4double[numState 77 fLastInitialVector = new G4double[numStateVars] ;
69 fLastFinalVector = new G4double[numStateVa 78 fLastFinalVector = new G4double[numStateVars] ;
70 79
71 fLastDyDx = new G4double[numStateVars]; 80 fLastDyDx = new G4double[numStateVars];
72 81
73 fMidVector = new G4double[numStateVars]; 82 fMidVector = new G4double[numStateVars];
74 fMidError = new G4double[numStateVars]; << 83 fMidError = new G4double[numStateVars];
75 84
76 if( primary ) 85 if( primary )
77 { 86 {
78 fAuxStepper = new G4DormandPrinceRK56(Eq << 87 fAuxStepper = new G4DormandPrinceRK56(EqRhs, numberOfVariables,
>> 88 !primary);
79 } 89 }
80 } 90 }
81 91
82 // Destructor <<
83 // <<
84 G4DormandPrinceRK56::~G4DormandPrinceRK56() <<
85 { <<
86 // clear all previously allocated memory f <<
87 92
88 delete [] ak2; << 93 //Destructor
89 delete [] ak3; << 94 G4DormandPrinceRK56::~G4DormandPrinceRK56(){
90 delete [] ak4; << 95 //clear all previously allocated memory for stepper and DistChord
91 delete [] ak5; << 96 delete[] ak2;
92 delete [] ak6; << 97 delete[] ak3;
93 delete [] ak7; << 98 delete[] ak4;
94 delete [] ak8; << 99 delete[] ak5;
95 delete [] ak9; << 100 delete[] ak6;
96 << 101 delete[] ak7;
97 delete [] ak10; << 102 delete[] ak8;
98 delete [] ak10_low; << 103 delete[] ak9;
99 delete [] ak11; << 104
100 delete [] ak12; << 105 delete[] ak10;
101 << 106 delete[] ak10_low;
102 delete [] yTemp; << 107 delete[] ak11;
103 delete [] yIn; << 108 delete[] ak12;
104 << 109
105 delete [] fLastInitialVector; << 110 delete[] yTemp;
106 delete [] fLastFinalVector; << 111 delete[] yIn;
107 delete [] fLastDyDx; << 112
108 delete [] fMidVector; << 113 delete[] fLastInitialVector;
109 delete [] fMidError; << 114 delete[] fLastFinalVector;
>> 115 delete[] fLastDyDx;
>> 116 delete[] fMidVector;
>> 117 delete[] fMidError;
110 118
111 delete fAuxStepper; 119 delete fAuxStepper;
>> 120
112 } 121 }
113 122
114 // Stepper << 123
115 // << 124 //Stepper :
>> 125
116 // Passing in the value of yInput[],the first 126 // Passing in the value of yInput[],the first time dydx[] and Step length
117 // Giving back yOut and yErr arrays for output 127 // Giving back yOut and yErr arrays for output and error respectively
118 // << 128
119 void G4DormandPrinceRK56::Stepper(const G4doub 129 void G4DormandPrinceRK56::Stepper(const G4double yInput[],
120 const G4doub << 130 const G4double dydx[],
121 G4doub << 131 G4double Step,
122 G4doub << 132 G4double yOut[],
123 G4doub << 133 G4double yErr[] )
124 // G4double 134 // G4double nextDydx[] ) -- Output:
125 // endpoint << 135 // endpoint DyDx ( for future FSAL version )
126 { 136 {
127 G4int i; 137 G4int i;
128 << 138
129 // The various constants defined on the ba << 139 //The various constants defined on the basis of butcher tableu
>> 140 const G4double //G4double - only once
>> 141
130 // Old Coefficients from 142 // Old Coefficients from
131 // P.J.Prince and J.R.Dormand, "High order << 143 // Table 1. RK6(5)8M
132 // Journal of Computational and Applied Ma << 144 //---Ref---
133 // << 145 //[P. J. Prince and J. R. Dormand, “High order embedded Runge-Kutta formulae,”
134 const G4double b21 = 1.0/10.0 , << 146 // Journal of Computational and Applied Mathematics, vol. 7, no. 1, pp. 67–75,
135 b31 = -2.0/81.0 , << 147 // Dec. 1980.
136 b32 = 20.0/81.0 , << 148 //----------------
137 << 149
138 b41 = 615.0/1372.0 , << 150 b21 = 1.0/10.0 ,
139 b42 = -270.0/343.0 , << 151
140 b43 = 1053.0/1372.0 , << 152 b31 = -2.0/81.0 ,
141 << 153 b32 = 20.0/81.0 ,
142 b51 = 3243.0/5500.0 , << 154
143 b52 = -54.0/55.0 , << 155 b41 = 615.0/1372.0 ,
144 b53 = 50949.0/71500.0 , << 156 b42 = -270.0/343.0 ,
145 b54 = 4998.0/17875.0 , << 157 b43 = 1053.0/1372.0 ,
146 << 158
147 b61 = -26492.0/37125.0 , << 159 b51 = 3243.0/5500.0 ,
148 b62 = 72.0/55.0 , << 160 b52 = -54.0/55.0 ,
149 b63 = 2808.0/23375.0 , << 161 b53 = 50949.0/71500.0 ,
150 b64 = -24206.0/37125.0 , << 162 b54 = 4998.0/17875.0 ,
151 b65 = 338.0/459.0 , << 163
152 << 164 b61 = -26492.0/37125.0 ,
153 b71 = 5561.0/2376.0 , << 165 b62 = 72.0/55.0 ,
154 b72 = -35.0/11.0 , << 166 b63 = 2808.0/23375.0 ,
155 b73 = -24117.0/31603.0 , << 167 b64 = -24206.0/37125.0 ,
156 b74 = 899983.0/200772.0 , << 168 b65 = 338.0/459.0 ,
157 b75 = -5225.0/1836.0 , << 169
158 b76 = 3925.0/4056.0 , << 170 b71 = 5561.0/2376.0 ,
159 << 171 b72 = -35.0/11.0 ,
160 b81 = 465467.0/266112.0 , << 172 b73 = -24117.0/31603.0 ,
161 b82 = -2945.0/1232.0 , << 173 b74 = 899983.0/200772.0 ,
162 b83 = -5610201.0/14158144.0 << 174 b75 = -5225.0/1836.0 ,
163 b84 = 10513573.0/3212352.0 << 175 b76 = 3925.0/4056.0 ,
164 b85 = -424325.0/205632.0 , << 176
165 b86 = 376225.0/454272.0 , << 177 b81 = 465467.0/266112.0 ,
166 b87 = 0.0 , << 178 b82 = -2945.0/1232.0 ,
167 << 179 b83 = -5610201.0/14158144.0 ,
168 c1 = 61.0/864.0 , << 180 b84 = 10513573.0/3212352.0 ,
169 c2 = 0.0 , << 181 b85 = -424325.0/205632.0 ,
170 c3 = 98415.0/321776.0 , << 182 b86 = 376225.0/454272.0 ,
171 c4 = 16807.0/146016.0 , << 183 b87 = 0.0 ,
172 c5 = 1375.0/7344.0 , << 184
173 c6 = 1375.0/5408.0 , << 185 c1 = 61.0/864.0 ,
174 c7 = -37.0/1120.0 , << 186 c2 = 0.0 ,
175 c8 = 1.0/10.0 , << 187 c3 = 98415.0/321776.0 ,
176 << 188 c4 = 16807.0/146016.0 ,
177 b91 = 61.0/864.0 , << 189 c5 = 1375.0/7344.0 ,
178 b92 = 0.0 , << 190 c6 = 1375.0/5408.0 ,
179 b93 = 98415.0/321776.0 , << 191 c7 = -37.0/1120.0 ,
180 b94 = 16807.0/146016.0 , << 192 c8 = 1.0/10.0 ,
181 b95 = 1375.0/7344.0 , << 193
182 b96 = 1375.0/5408.0 , << 194 b91 = 61.0/864.0 ,
183 b97 = -37.0/1120.0 , << 195 b92 = 0.0 ,
184 b98 = 1.0/10.0 , << 196 b93 = 98415.0/321776.0 ,
185 << 197 b94 = 16807.0/146016.0 ,
186 dc1 = c1 - 821.0/10800.0 << 198 b95 = 1375.0/7344.0 ,
187 dc2 = c2 - 0.0 , << 199 b96 = 1375.0/5408.0 ,
188 dc3 = c3 - 19683.0/71825, << 200 b97 = -37.0/1120.0 ,
189 dc4 = c4 - 175273.0/912600 << 201 b98 = 1.0/10.0 ,
190 dc5 = c5 - 395.0/3672.0 , << 202
191 dc6 = c6 - 785.0/2704.0 , << 203 dc1 = c1 - 821.0/10800.0 ,
192 dc7 = c7 - 3.0/50.0 , << 204 dc2 = c2 - 0.0 ,
193 dc8 = c8 - 0.0 , << 205 dc3 = c3 - 19683.0/71825,
194 dc9 = 0.0; << 206 dc4 = c4 - 175273.0/912600.0 ,
>> 207 dc5 = c5 - 395.0/3672.0 ,
>> 208 dc6 = c6 - 785.0/2704.0 ,
>> 209 dc7 = c7 - 3.0/50.0 ,
>> 210 dc8 = c8 - 0.0 ,
>> 211 dc9 = 0.0;
195 212
196 213
197 // New Coefficients obtained from 214 // New Coefficients obtained from
198 // Table 3 RK6(5)9FM with corrected coefficien << 215 // Table 3 RK6(5)9FM with corrected coefficients
199 // << 216 //---Ref---
200 // T. S. Baker, J. R. Dormand, J. P. Gilmor 217 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
201 // "Continuous approximation with embedded << 218 // “Continuous approximation with embedded Runge-Kutta methods,”
202 // Applied Numerical Mathematics, vol. 22, << 219 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51–62, 1996.
203 // << 220 //------------------------------------
>> 221
204 // b21 = 1.0/9.0 , 222 // b21 = 1.0/9.0 ,
205 // 223 //
206 // b31 = 1.0/24.0 , 224 // b31 = 1.0/24.0 ,
207 // b32 = 1.0/8.0 , 225 // b32 = 1.0/8.0 ,
208 // 226 //
209 // b41 = 1.0/16.0 , 227 // b41 = 1.0/16.0 ,
210 // b42 = 0.0 , 228 // b42 = 0.0 ,
211 // b43 = 3.0/16.0 , 229 // b43 = 3.0/16.0 ,
212 // 230 //
213 // b51 = 280.0/729.0 , 231 // b51 = 280.0/729.0 ,
214 // b52 = 0.0 , 232 // b52 = 0.0 ,
215 // b53 = -325.0/243.0 , 233 // b53 = -325.0/243.0 ,
216 // b54 = 1100.0/729.0 , 234 // b54 = 1100.0/729.0 ,
217 // 235 //
218 // b61 = 6127.0/14680.0 , 236 // b61 = 6127.0/14680.0 ,
219 // b62 = 0.0 , 237 // b62 = 0.0 ,
220 // b63 = -1077.0/734.0 , 238 // b63 = -1077.0/734.0 ,
221 // b64 = 6494.0/4037.0 , 239 // b64 = 6494.0/4037.0 ,
222 // b65 = -9477.0/161480.0 , 240 // b65 = -9477.0/161480.0 ,
223 // 241 //
224 // b71 = -13426273320.0/14809773769.0 , 242 // b71 = -13426273320.0/14809773769.0 ,
225 // b72 = 0.0 , 243 // b72 = 0.0 ,
226 // b73 = 4192558704.0/2115681967.0 , 244 // b73 = 4192558704.0/2115681967.0 ,
227 // b74 = 14334750144.0/14809773769.0 , 245 // b74 = 14334750144.0/14809773769.0 ,
228 // b75 = 117092732328.0/14809773769.0 , 246 // b75 = 117092732328.0/14809773769.0 ,
229 // b76 = -361966176.0/40353607.0 , 247 // b76 = -361966176.0/40353607.0 ,
230 // 248 //
231 // b81 = -2340689.0/1901060.0 , 249 // b81 = -2340689.0/1901060.0 ,
232 // b82 = 0.0 , 250 // b82 = 0.0 ,
233 // b83 = 31647.0/13579.0 , 251 // b83 = 31647.0/13579.0 ,
234 // b84 = 253549596.0/149518369.0 , 252 // b84 = 253549596.0/149518369.0 ,
235 // b85 = 10559024082.0/977620105.0 , 253 // b85 = 10559024082.0/977620105.0 ,
236 // b86 = -152952.0/12173.0 , 254 // b86 = -152952.0/12173.0 ,
237 // b87 = -5764801.0/186010396.0 , 255 // b87 = -5764801.0/186010396.0 ,
238 // 256 //
239 // b91 = 203.0/2880.0 , 257 // b91 = 203.0/2880.0 ,
240 // b92 = 0.0 , 258 // b92 = 0.0 ,
241 // b93 = 0.0 , 259 // b93 = 0.0 ,
242 // b94 = 30208.0/70785.0 , 260 // b94 = 30208.0/70785.0 ,
243 // b95 = 177147.0/164560.0 , 261 // b95 = 177147.0/164560.0 ,
244 // b96 = -536.0/705.0 , 262 // b96 = -536.0/705.0 ,
245 // b97 = 1977326743.0/3619661760.0 , 263 // b97 = 1977326743.0/3619661760.0 ,
246 // b98 = -259.0/720.0 , 264 // b98 = -259.0/720.0 ,
247 // 265 //
248 // 266 //
249 // dc1 = 36567.0/458800.0 - b91, 267 // dc1 = 36567.0/458800.0 - b91,
250 // dc2 = 0.0 - b92, 268 // dc2 = 0.0 - b92,
251 // dc3 = 0.0 - b93, 269 // dc3 = 0.0 - b93,
252 // dc4 = 9925984.0/27063465.0 - b94, 270 // dc4 = 9925984.0/27063465.0 - b94,
253 // dc5 = 85382667.0/117968950.0 - b95, 271 // dc5 = 85382667.0/117968950.0 - b95,
254 // dc6 = - 310378.0/808635.0 - b96 , 272 // dc6 = - 310378.0/808635.0 - b96 ,
255 // dc7 = 262119736669.0/345979336560.0 - b 273 // dc7 = 262119736669.0/345979336560.0 - b97,
256 // dc8 = - 1.0/2.0 - b98 , 274 // dc8 = - 1.0/2.0 - b98 ,
257 // dc9 = -101.0/2294.0 ; 275 // dc9 = -101.0/2294.0 ;
>> 276
258 277
259 // end of declaration << 278 //end of declaration
>> 279
260 280
261 const G4int numberOfVariables = GetNumberO << 281 const G4int numberOfVariables= this->GetNumberOfVariables();
262 282
263 // The number of variables to be integrate 283 // The number of variables to be integrated over
264 // << 284 yOut[7] = yTemp[7] = yIn[7];
265 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; <<
266 <<
267 // Saving yInput because yInput and yOut 285 // Saving yInput because yInput and yOut can be aliases for same array
268 // << 286
269 for(i=0; i<numberOfVariables; ++i) << 287 for(i=0;i<numberOfVariables;i++)
270 { 288 {
271 yIn[i]=yInput[i]; 289 yIn[i]=yInput[i];
272 } 290 }
273 // RightHandSide(yIn, dydx) ; // 1st Stage <<
274 291
275 for(i=0; i<numberOfVariables; ++i) << 292
>> 293
>> 294 // RightHandSide(yIn, dydx) ;
>> 295 // 1st Step - Not doing, getting passed
>> 296
>> 297 for(i=0;i<numberOfVariables;i++)
276 { 298 {
277 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; 299 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
278 } 300 }
279 RightHandSide(yTemp, ak2) ; / 301 RightHandSide(yTemp, ak2) ; // 2nd Stage
280 302
281 for(i=0; i<numberOfVariables; ++i) << 303 for(i=0;i<numberOfVariables;i++)
282 { 304 {
283 yTemp[i] = yIn[i] + Step*(b31*dydx[i] 305 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
284 } 306 }
285 RightHandSide(yTemp, ak3) ; / 307 RightHandSide(yTemp, ak3) ; // 3rd Stage
286 308
287 for(i=0; i<numberOfVariables; ++i) << 309 for(i=0;i<numberOfVariables;i++)
288 { 310 {
289 yTemp[i] = yIn[i] + Step*(b41*dydx[i] 311 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
290 } 312 }
291 RightHandSide(yTemp, ak4) ; / 313 RightHandSide(yTemp, ak4) ; // 4th Stage
292 314
293 for(i=0; i<numberOfVariables; ++i) << 315 for(i=0;i<numberOfVariables;i++)
294 { 316 {
295 yTemp[i] = yIn[i] + Step*(b51*dydx[i] 317 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
296 b54*ak4[i]) 318 b54*ak4[i]) ;
297 } 319 }
298 RightHandSide(yTemp, ak5) ; / 320 RightHandSide(yTemp, ak5) ; // 5th Stage
299 321
300 for(i=0; i<numberOfVariables; ++i) << 322 for(i=0;i<numberOfVariables;i++)
301 { 323 {
302 yTemp[i] = yIn[i] + Step*(b61*dydx[i] 324 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
303 b64*ak4[i] + 325 b64*ak4[i] + b65*ak5[i]) ;
304 } 326 }
305 RightHandSide(yTemp, ak6) ; / 327 RightHandSide(yTemp, ak6) ; // 6th Stage
306 328
307 for(i=0; i<numberOfVariables; ++i) << 329 for(i=0;i<numberOfVariables;i++)
308 { 330 {
309 yTemp[i] = yIn[i] + Step*(b71*dydx[i] 331 yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
310 b74*ak4[i] + 332 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
311 } 333 }
312 RightHandSide(yTemp, ak7); / << 334 RightHandSide(yTemp, ak7); //7th Stage
313 335
314 for(i=0; i<numberOfVariables; ++i) << 336 for(i=0;i<numberOfVariables;i++)
315 { 337 {
316 yTemp[i] = yIn[i] + Step*(b81*dydx[i] 338 yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
317 b84*ak4[i] + 339 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
318 b87*ak7[i]); 340 b87*ak7[i]);
319 } 341 }
320 RightHandSide(yTemp, ak8); / << 342 RightHandSide(yTemp, ak8); //8th Stage
321 343
322 for(i=0; i<numberOfVariables; ++i) << 344 for(i=0;i<numberOfVariables;i++)
323 { 345 {
324 yOut[i] = yIn[i] + Step*(b91*dydx[i] + 346 yOut[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
325 b94*ak4[i] + 347 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
326 b97*ak7[i] + 348 b97*ak7[i] + b98*ak8[i] );
327 } 349 }
328 RightHandSide(yOut, ak9); / << 350 RightHandSide(yOut, ak9); //9th Stage
329 351
330 352
331 for(i=0; i<numberOfVariables; ++i) << 353 for(i=0;i<numberOfVariables;i++)
332 { 354 {
>> 355
333 // Estimate error as difference betwee 356 // Estimate error as difference between 5th and
334 // 6th order methods 357 // 6th order methods
335 // << 358
336 yErr[i] = Step*( dc1*dydx[i] + dc2*ak 359 yErr[i] = Step*( dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
337 + dc5*ak5[i] + dc6*ak6 360 + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
338 + dc9*ak9[i] ) ; 361 + dc9*ak9[i] ) ;
339 362
340 // Saving 'estimated' derivative at en << 363 // - Saving 'estimated' derivative at end-point
341 // nextDydx[i] = ak9[i]; 364 // nextDydx[i] = ak9[i];
342 365
343 // Store Input and Final values, for p 366 // Store Input and Final values, for possible use in calculating chord
344 // <<
345 fLastInitialVector[i] = yIn[i] ; 367 fLastInitialVector[i] = yIn[i] ;
346 fLastFinalVector[i] = yOut[i]; 368 fLastFinalVector[i] = yOut[i];
347 fLastDyDx[i] = dydx[i]; 369 fLastDyDx[i] = dydx[i];
>> 370
348 } 371 }
349 372
350 fLastStepLength = Step; 373 fLastStepLength = Step;
351 374
352 return ; 375 return ;
353 } 376 }
354 377
355 // DistChord << 378
356 // << 379 //The following has not been tested
>> 380
>> 381 //The DistChord() function fot the class - must define it here.
357 G4double G4DormandPrinceRK56::DistChord() con 382 G4double G4DormandPrinceRK56::DistChord() const
358 { 383 {
359 G4double distLine, distChord; 384 G4double distLine, distChord;
360 G4ThreeVector initialPoint, finalPoint, mi 385 G4ThreeVector initialPoint, finalPoint, midPoint;
361 386
362 // Store last initial and final points << 387 // Store last initial and final points (they will be overwritten in self-Stepper call!)
363 // (they will be overwritten in self-Stepp <<
364 // <<
365 initialPoint = G4ThreeVector( fLastInitial 388 initialPoint = G4ThreeVector( fLastInitialVector[0],
366 fLastInitialV 389 fLastInitialVector[1], fLastInitialVector[2]);
367 finalPoint = G4ThreeVector( fLastFinalVe 390 finalPoint = G4ThreeVector( fLastFinalVector[0],
368 fLastFinalVec 391 fLastFinalVector[1], fLastFinalVector[2]);
369 392
370 // Do half a step using StepNoErr 393 // Do half a step using StepNoErr
371 394
372 fAuxStepper->Stepper( fLastInitialVector, 395 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
373 fMidVector, fMidErr 396 fMidVector, fMidError );
374 397
375 midPoint = G4ThreeVector( fMidVector[0], f 398 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
376 399
377 // Use stored values of Initial and Endpoi 400 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
378 // distance of Chord << 401 // distance of Chord
379 // << 402
>> 403
380 if (initialPoint != finalPoint) 404 if (initialPoint != finalPoint)
381 { 405 {
382 distLine = G4LineSection::Distline( m << 406 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
383 distChord = distLine; 407 distChord = distLine;
384 } 408 }
385 else 409 else
386 { 410 {
387 distChord = (midPoint-initialPoint).ma 411 distChord = (midPoint-initialPoint).mag();
388 } 412 }
389 return distChord; 413 return distChord;
390 } 414 }
391 415
>> 416
392 // The following interpolation scheme has been 417 // The following interpolation scheme has been obtained from
393 // Table 5. The RK6(5)9FM process and associat 418 // Table 5. The RK6(5)9FM process and associated dense formula
394 // << 419 //---Ref---
395 // J. R. Dormand, M. A. Lockyer, N. E. McGorri << 420 // J. R. Dormand, M. A. Lockyer, N. E. McGorrigan, and P. J. Prince,
396 // "Global error estimation with runge-kutta t << 421 // “Global error estimation with runge-kutta triples,”
397 // Computers & Mathematics with Applications, << 422 // Computers & Mathematics with Applications, vol. 18, no. 9, pp. 835–846, 1989.
398 // << 423 //-----------------------------
>> 424
399 // Fifth order interpolant with one extra func 425 // Fifth order interpolant with one extra function evaluation per step
400 // << 426
401 void G4DormandPrinceRK56::SetupInterpolate_low 427 void G4DormandPrinceRK56::SetupInterpolate_low( const G4double yInput[],
402 << 428 const G4double dydx[],
403 << 429 const G4double Step ){
404 { << 430 const G4int numberOfVariables= this->GetNumberOfVariables();
405 const G4int numberOfVariables= this->GetNu << 431
406 << 432 G4double
407 G4double b_101 = 33797.0/460800.0 , << 433 b_101 = 33797.0/460800.0 ,
408 b_102 = 0. , << 434 b_102 = 0. ,
409 b_103 = 0. , << 435 b_103 = 0. ,
410 b_104 = 27757.0/70785.0 , << 436 b_104 = 27757.0/70785.0 ,
411 b_105 = 7923501.0/26329600.0 , << 437 b_105 = 7923501.0/26329600.0 ,
412 b_106 = -927.0/3760.0 , << 438 b_106 = -927.0/3760.0 ,
413 b_107 = -3314760575.0/23165835264 << 439 b_107 = -3314760575.0/23165835264.0 ,
414 b_108 = 2479.0/23040.0 , << 440 b_108 = 2479.0/23040.0 ,
415 b_109 = 1.0/64.0 ; << 441 b_109 = 1.0/64.0 ;
416 442
417 for(G4int i=0; i<numberOfVariables; ++i) << 443 for(int i=0;i<numberOfVariables;i++)
418 { 444 {
419 yIn[i]=yInput[i]; << 445 yIn[i]=yInput[i];
420 } 446 }
421 447
422 448
423 for(G4int i=0; i<numberOfVariables; ++i) << 449 for(int i=0;i<numberOfVariables;i++)
424 { 450 {
425 yTemp[i] = yIn[i] + Step*(b_101*dydx[i] << 451 yTemp[i] = yIn[i] + Step*(b_101*dydx[i] + b_102*ak2[i] + b_103*ak3[i] +
426 b_104*ak4[i] + << 452 b_104*ak4[i] + b_105*ak5[i] + b_106*ak6[i] +
427 b_107*ak7[i] + << 453 b_107*ak7[i] + b_108*ak8[i] + b_109*ak9[i]);
428 } 454 }
429 RightHandSide(yTemp, ak10_low); / << 455 RightHandSide(yTemp, ak10_low); //10th Stage
430 } 456 }
431 457
432 void G4DormandPrinceRK56::Interpolate_low( con 458 void G4DormandPrinceRK56::Interpolate_low( const G4double yInput[],
433 con << 459 const G4double dydx[],
434 con << 460 const G4double Step,
435 << 461 G4double yOut[],
436 << 462 G4double tau ){
437 { << 463 {
438 G4double bf1, bf4, bf5, bf6, bf7, bf8, bf9, <<
439 464
440 G4double tau0 = tau; << 465 G4double
441 const G4int numberOfVariables= this->GetNumb << 466 bf1, bf4, bf5, bf6, bf7, bf8, bf9, bf10;
>> 467
>> 468 G4double tau0 = tau;
>> 469 const G4int numberOfVariables= this->GetNumberOfVariables();
442 470
443 for(G4int i=0; i<numberOfVariables; ++i) << 471 for(int i=0;i<numberOfVariables;i++)
444 { << 472 {
445 yIn[i]=yInput[i]; << 473 yIn[i]=yInput[i];
446 } << 474 }
447 << 475
448 G4double tau_2 = tau0*tau0 , << 476 G4double
449 tau_3 = tau0*tau_2, << 477 tau_2 = tau0*tau0 ,
450 tau_4 = tau_2*tau_2; << 478 tau_3 = tau0*tau_2,
451 << 479 tau_4 = tau_2*tau_2;
452 // bf2 = bf3 = 0.0 << 480
453 bf1 = (66480.0*tau_4-206243.0*tau_3+237786.0 << 481 //bf2 = bf3 = 0
454 / 28800.0 ; << 482 bf1 = (66480.0*tau_4 - 206243.0*tau_3 + 237786.0*tau_2 - 124793.0*tau + 28800.0)/28800.0 ,
455 bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*ta << 483 bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*tau_2 + 119706.0*tau -40973.0)/70785.0 ,
456 / 70785.0 ; << 484 bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*tau_2 + 34786.0*tau - 9293.0)/1645600.0 ,
457 bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*t << 485 bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 23786.0*tau - 6533.0)/705.0 ,
458 / 1645600.0 ; << 486 bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797.0*tau_2 + 17574.0*tau - 1927.0)/7239323520.0 ,
459 bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 2 << 487 bf8 = 37.0*tau*(336.0*tau_3 - 661.0*tau_2 + 342.0*tau -31.0)/1440.0 ,
460 / 705.0 ; << 488 bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau +3.0)/4.0 ,
461 bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797. << 489 bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0*tau - 1.0) ;
462 / 7239323520.0 ; << 490
463 bf8 = 37.0*tau*(336.0*tau_3 - 661.0*tau_2 + << 491 for( int i=0; i<numberOfVariables; i++){
464 / 1440.0 ; << 492 yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i] + bf4*ak4[i] + bf5*ak5[i] +
465 bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau + << 493 bf6*ak6[i] + bf7*ak7[i] + bf8*ak8[i] +
466 / 4.0 ; << 494 bf9*ak9[i] + bf10*ak10_low[i] ) ;
467 bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0* << 495 }
468 << 496
469 for( G4int i=0; i<numberOfVariables; ++i) << 497
470 { << 498
471 yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i] << 499 }
472 bf6*ak6[i] + << 500
473 bf9*ak9[i] + <<
474 } <<
475 } 501 }
476 502
477 // The following scheme and set of coefficient << 503 //The following scheme and set of coefficients have been obtained from
478 // Table 2. Sixth order dense formula based on << 504 //Table 2. Sixth order dense formula based on linear optimisation for RK6(5)9FM
479 // RK6(5)9FM with extra stages C1O= 1/2, C11 = << 505 //with extra stages C1O= 1/2, C11 =1/6, c12= 5/12
480 // << 506
481 // T. S. Baker, J. R. Dormand, J. P. Gilmore, << 507 //---Ref---
482 // "Continuous approximation with embedded Run << 508 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
483 // Applied Numerical Mathematics, vol. 22, no. << 509 // “Continuous approximation with embedded Runge-Kutta methods,”
484 // << 510 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51–62, 1996.
485 // --- Sixth order interpolant with 3 addition << 511 //--------------------
486 // << 512
487 // Function for calculating the additional sta << 513
488 // << 514 // --- Sixth order interpolant with 3 additional stages per step ---
>> 515
>> 516 //Function for calculating the additional stages :
489 void G4DormandPrinceRK56::SetupInterpolate_hig 517 void G4DormandPrinceRK56::SetupInterpolate_high( const G4double yInput[],
490 << 518 const G4double dydx[],
491 << 519 const G4double Step ){
492 { <<
493 // Coefficients for the additional stages <<
494 // <<
495 G4double b101 = 33797.0/460800.0 , <<
496 b102 = 0.0 , <<
497 b103 = 0.0 , <<
498 b104 = 27757.0/70785.0 , <<
499 b105 = 7923501.0/26329600.0 , <<
500 b106 = -927.0/3760.0 , <<
501 b107 = -3314760575.0/23165835264. <<
502 b108 = 2479.0/23040.0 , <<
503 b109 = 1.0/64.0 , <<
504 <<
505 b111 = 5843.0/76800.0 , <<
506 b112 = 0.0 , <<
507 b113 = 0.0 , <<
508 b114 = 464.0/2673.0 , <<
509 b115 = 353997.0/1196800.0 , <<
510 b116 = -15068.0/57105.0 , <<
511 b117 = -282475249.0/3644974080.0 <<
512 b118 = 8678831.0/156245760.0 , <<
513 b119 = 116113.0/11718432.0 , <<
514 b1110 = -25.0/243.0 , <<
515 <<
516 b121 = 15088049.0/199065600.0 , <<
517 b122 = 0.0 , <<
518 b123 = 0.0 , <<
519 b124 = 2.0/5.0 , <<
520 b125 = 92222037.0/268083200.0 , <<
521 b126 = -433420501.0/1528586640.0 <<
522 b127 = -11549242677007.0/83630285 <<
523 b128 = 2725085557.0/26167173120. <<
524 b129 = 235429367.0/16354483200.0 <<
525 b1210 = -90924917.0/1040739840.0 <<
526 b1211 = -271149.0/21414400.0 ; <<
527 520
528 const G4int numberOfVariables = GetNumberO << 521 //Coefficients for the additional stages :
529 522
530 // Saving yInput because yInput and yOut c << 523 G4double
531 // << 524 b101 = 33797.0/460800.0 ,
532 for(G4int i=0; i<numberOfVariables; ++i) << 525 b102 = 0.0 ,
>> 526 b103 = 0.0 ,
>> 527 b104 = 27757.0/70785.0 ,
>> 528 b105 = 7923501.0/26329600.0 ,
>> 529 b106 = -927.0/3760.0 ,
>> 530 b107 = -3314760575.0/23165835264.0 ,
>> 531 b108 = 2479.0/23040.0 ,
>> 532 b109 = 1.0/64.0 ,
>> 533
>> 534 b111 = 5843.0/76800.0 ,
>> 535 b112 = 0.0 ,
>> 536 b113 = 0.0 ,
>> 537 b114 = 464.0/2673.0 ,
>> 538 b115 = 353997.0/1196800.0 ,
>> 539 b116 = -15068.0/57105.0 ,
>> 540 b117 = -282475249.0/3644974080.0 ,
>> 541 b118 = 8678831.0/156245760.0 ,
>> 542 b119 = 116113.0/11718432.0 ,
>> 543 b1110 = -25.0/243.0 ,
>> 544
>> 545 b121 = 15088049.0/199065600.0 ,
>> 546 b122 = 0.0 ,
>> 547 b123 = 0.0 ,
>> 548 b124 = 2.0/5.0 ,
>> 549 b125 = 92222037.0/268083200.0 ,
>> 550 b126 = -433420501.0/1528586640.0 ,
>> 551 b127 = -11549242677007.0/83630285291520.0 ,
>> 552 b128 = 2725085557.0/26167173120.0 ,
>> 553 b129 = 235429367.0/16354483200.0 ,
>> 554 b1210 = -90924917.0/1040739840.0 ,
>> 555 b1211 = -271149.0/21414400.0 ;
>> 556
>> 557
>> 558 const G4int numberOfVariables= this->GetNumberOfVariables();
>> 559
>> 560 // Saving yInput because yInput and yOut can be aliases for same array
>> 561 for(int i=0;i<numberOfVariables;i++)
533 { 562 {
534 yIn[i]=yInput[i]; << 563 yIn[i]=yInput[i];
535 } 564 }
>> 565
536 yTemp[7] = yIn[7]; 566 yTemp[7] = yIn[7];
>> 567
537 568
538 // Evaluate the extra stages << 569
539 // << 570 //Evaluate the extra stages :
540 for(G4int i=0; i<numberOfVariables; ++i) << 571
>> 572 for(int i=0;i<numberOfVariables;i++)
541 { 573 {
542 yTemp[i] = yIn[i] + Step*(b101*dydx[i] 574 yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
543 b104*ak4[i] 575 b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
544 b107*ak7[i] 576 b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
545 } 577 }
546 RightHandSide(yTemp, ak10); << 578 RightHandSide(yTemp, ak10); //10th Stage
547 579
548 for(G4int i=0; i<numberOfVariables; ++i) << 580 for(int i=0;i<numberOfVariables;i++)
549 { 581 {
550 yTemp[i] = yIn[i] + Step*(b111*dydx[i] 582 yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
551 b114*ak4[i] 583 b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
552 b117*ak7[i] 584 b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
553 b1110*ak10[i 585 b1110*ak10[i]);
554 } 586 }
555 RightHandSide(yTemp, ak11); << 587 RightHandSide(yTemp, ak11); //11th Stage
556 588
557 for(G4int i=0; i<numberOfVariables; ++i) << 589 for(int i=0;i<numberOfVariables;i++)
558 { 590 {
559 yTemp[i] = yIn[i] + Step*(b121*dydx[i] 591 yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
560 b124*ak4[i] 592 b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
561 b127*ak7[i] 593 b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
562 b1210*ak10[i 594 b1210*ak10[i] + b1211*ak11[i]);
563 } 595 }
564 RightHandSide(yTemp, ak12); << 596 RightHandSide(yTemp, ak12); //12th Stage
>> 597
565 } 598 }
566 599
567 // Function to interpolate to tau(passed in) f << 600
568 // << 601
>> 602 //Function to interpolate to tau(passed in) fraction of the step
569 void G4DormandPrinceRK56::Interpolate_high( co 603 void G4DormandPrinceRK56::Interpolate_high( const G4double yInput[],
570 co << 604 const G4double dydx[],
571 co << 605 const G4double Step,
572 << 606 G4double yOut[],
573 << 607 G4double tau )
574 { 608 {
575 // Define the coefficients for the polynom << 609
576 // << 610 //Define the coefficients for the polynomials
577 G4double bi[13][6], b[13]; 611 G4double bi[13][6], b[13];
578 G4int numberOfVariables = GetNumberOfVaria << 612 G4int numberOfVariables = this->GetNumberOfVariables();
579 613
580 << 614
581 // COEFFICIENTS OF bi[ 1] 615 // COEFFICIENTS OF bi[ 1]
582 bi[1][0] = 1.0 , 616 bi[1][0] = 1.0 ,
583 bi[1][1] = -18487.0/2880.0 , 617 bi[1][1] = -18487.0/2880.0 ,
584 bi[1][2] = 139189.0/7200.0 , 618 bi[1][2] = 139189.0/7200.0 ,
585 bi[1][3] = -53923.0/1800.0 , 619 bi[1][3] = -53923.0/1800.0 ,
586 bi[1][4] = 13811.0/600, 620 bi[1][4] = 13811.0/600,
587 bi[1][5] = -2071.0/300, 621 bi[1][5] = -2071.0/300,
588 // -------------------------------------- 622 // --------------------------------------------------------
589 // 623 //
590 // COEFFICIENTS OF bi[2] 624 // COEFFICIENTS OF bi[2]
591 bi[2][0] = 0.0 , 625 bi[2][0] = 0.0 ,
592 bi[2][1] = 0.0 , 626 bi[2][1] = 0.0 ,
593 bi[2][2] = 0.0 , 627 bi[2][2] = 0.0 ,
594 bi[2][3] = 0.0 , 628 bi[2][3] = 0.0 ,
595 bi[2][4] = 0.0 , 629 bi[2][4] = 0.0 ,
596 bi[2][5] = 0.0 , 630 bi[2][5] = 0.0 ,
597 // -------------------------------------- 631 // --------------------------------------------------------
598 // 632 //
599 // COEFFICIENTS OF bi[3] 633 // COEFFICIENTS OF bi[3]
600 bi[3][0] = 0.0 , 634 bi[3][0] = 0.0 ,
601 bi[3][1] = 0.0 , 635 bi[3][1] = 0.0 ,
602 bi[3][2] = 0.0 , 636 bi[3][2] = 0.0 ,
603 bi[3][3] = 0.0 , 637 bi[3][3] = 0.0 ,
604 bi[3][4] = 0.0 , 638 bi[3][4] = 0.0 ,
605 bi[3][5] = 0.0 , 639 bi[3][5] = 0.0 ,
606 // -------------------------------------- 640 // --------------------------------------------------------
607 // 641 //
608 // COEFFICIENTS OF bi[4] 642 // COEFFICIENTS OF bi[4]
609 bi[4][0] = 0.0 , 643 bi[4][0] = 0.0 ,
610 bi[4][1] = -30208.0/14157.0 , 644 bi[4][1] = -30208.0/14157.0 ,
611 bi[4][2] = 1147904.0/70785.0 , 645 bi[4][2] = 1147904.0/70785.0 ,
612 bi[4][3] = -241664.0/5445.0 , 646 bi[4][3] = -241664.0/5445.0 ,
613 bi[4][4] = 241664.0/4719.0 , 647 bi[4][4] = 241664.0/4719.0 ,
614 bi[4][5] = -483328.0/23595.0 , 648 bi[4][5] = -483328.0/23595.0 ,
615 // -------------------------------------- 649 // --------------------------------------------------------
616 // 650 //
617 // COEFFICIENTS OF bi[5] 651 // COEFFICIENTS OF bi[5]
618 bi[5][0] = 0.0 , 652 bi[5][0] = 0.0 ,
619 bi[5][1] = -177147.0/32912.0 , 653 bi[5][1] = -177147.0/32912.0 ,
620 bi[5][2] = 3365793.0/82280.0 , 654 bi[5][2] = 3365793.0/82280.0 ,
621 bi[5][3] = -2302911.0/20570.0 , 655 bi[5][3] = -2302911.0/20570.0 ,
622 bi[5][4] = 531441.0/4114.0 , 656 bi[5][4] = 531441.0/4114.0 ,
623 bi[5][5] = -531441.0/10285.0 , 657 bi[5][5] = -531441.0/10285.0 ,
624 // -------------------------------------- 658 // --------------------------------------------------------
625 // 659 //
626 // COEFFICIENTS OF bi[6] 660 // COEFFICIENTS OF bi[6]
627 bi[6][0] = 0.0 , 661 bi[6][0] = 0.0 ,
628 bi[6][1] = 536.0/141.0 , 662 bi[6][1] = 536.0/141.0 ,
629 bi[6][2] = -20368.0/705.0 , 663 bi[6][2] = -20368.0/705.0 ,
630 bi[6][3] = 55744.0/705.0 , 664 bi[6][3] = 55744.0/705.0 ,
631 bi[6][4] = -4288.0/47.0 , 665 bi[6][4] = -4288.0/47.0 ,
632 bi[6][5] = 8576.0/235, 666 bi[6][5] = 8576.0/235,
633 // -------------------------------------- 667 // --------------------------------------------------------
634 // 668 //
635 // COEFFICIENTS OF bi[7] 669 // COEFFICIENTS OF bi[7]
636 bi[7][0] = 0.0 , 670 bi[7][0] = 0.0 ,
637 bi[7][1] = -1977326743.0/723932352.0 , 671 bi[7][1] = -1977326743.0/723932352.0 ,
638 bi[7][2] = 37569208117.0/1809830880.0 , 672 bi[7][2] = 37569208117.0/1809830880.0 ,
639 bi[7][3] = -1977326743.0/34804440.0 , 673 bi[7][3] = -1977326743.0/34804440.0 ,
640 bi[7][4] = 1977326743.0/30163848.0 , 674 bi[7][4] = 1977326743.0/30163848.0 ,
641 bi[7][5] = -1977326743.0/75409620.0 , 675 bi[7][5] = -1977326743.0/75409620.0 ,
642 // -------------------------------------- 676 // --------------------------------------------------------
643 // 677 //
644 // COEFFICIENTS OF bi[8] 678 // COEFFICIENTS OF bi[8]
645 bi[8][0] = 0.0 , 679 bi[8][0] = 0.0 ,
646 bi[8][1] = 259.0/144.0 , 680 bi[8][1] = 259.0/144.0 ,
647 bi[8][2] = -4921.0/360.0 , 681 bi[8][2] = -4921.0/360.0 ,
648 bi[8][3] = 3367.0/90.0 , 682 bi[8][3] = 3367.0/90.0 ,
649 bi[8][4] = -259.0/6.0 , 683 bi[8][4] = -259.0/6.0 ,
650 bi[8][5] = 259.0/15.0 , 684 bi[8][5] = 259.0/15.0 ,
651 // -------------------------------------- 685 // --------------------------------------------------------
652 // 686 //
653 // COEFFICIENTS OF bi[9] 687 // COEFFICIENTS OF bi[9]
654 bi[9][0] = 0.0 , 688 bi[9][0] = 0.0 ,
655 bi[9][1] = 62.0/105.0 , 689 bi[9][1] = 62.0/105.0 ,
656 bi[9][2] = -2381.0/525.0 , 690 bi[9][2] = -2381.0/525.0 ,
657 bi[9][3] = 949.0/75.0 , 691 bi[9][3] = 949.0/75.0 ,
658 bi[9][4] = -2636.0/175.0 , 692 bi[9][4] = -2636.0/175.0 ,
659 bi[9][5] = 1112.0/175.0 , 693 bi[9][5] = 1112.0/175.0 ,
660 // -------------------------------------- 694 // --------------------------------------------------------
661 // 695 //
662 // COEFFICIENTS OF bi[10] 696 // COEFFICIENTS OF bi[10]
663 bi[10][0] = 0.0 , 697 bi[10][0] = 0.0 ,
664 bi[10][1] = 43.0/3.0 , 698 bi[10][1] = 43.0/3.0 ,
665 bi[10][2] = -1534.0/15.0 , 699 bi[10][2] = -1534.0/15.0 ,
666 bi[10][3] = 3767.0/15.0 , 700 bi[10][3] = 3767.0/15.0 ,
667 bi[10][4] = -1264.0/5.0 , 701 bi[10][4] = -1264.0/5.0 ,
668 bi[10][5] = 448.0/5.0 , 702 bi[10][5] = 448.0/5.0 ,
669 // -------------------------------------- 703 // --------------------------------------------------------
670 // 704 //
671 // COEFFICIENTS OF bi[11] 705 // COEFFICIENTS OF bi[11]
672 bi[11][0] = 0.0 , 706 bi[11][0] = 0.0 ,
673 bi[11][1] = 63.0/5.0 , 707 bi[11][1] = 63.0/5.0 ,
674 bi[11][2] = -1494.0/25.0 , 708 bi[11][2] = -1494.0/25.0 ,
675 bi[11][3] = 2907.0/25.0 , 709 bi[11][3] = 2907.0/25.0 ,
676 bi[11][4] = -2592.0/25.0 , 710 bi[11][4] = -2592.0/25.0 ,
677 bi[11][5] = 864.0/25.0 , 711 bi[11][5] = 864.0/25.0 ,
678 // -------------------------------------- 712 // --------------------------------------------------------
679 // 713 //
680 // COEFFICIENTS OF bi[12] 714 // COEFFICIENTS OF bi[12]
681 bi[12][0] = 0.0 , 715 bi[12][0] = 0.0 ,
682 bi[12][1] = -576.0/35.0 , 716 bi[12][1] = -576.0/35.0 ,
683 bi[12][2] = 19584.0/175.0 , 717 bi[12][2] = 19584.0/175.0 ,
684 bi[12][3] = -6336.0/25.0 , 718 bi[12][3] = -6336.0/25.0 ,
685 bi[12][4] = 41472.0/175.0 , 719 bi[12][4] = 41472.0/175.0 ,
686 bi[12][5] = -13824.0/175.0 ; 720 bi[12][5] = -13824.0/175.0 ;
687 // -------------------------------------- 721 // --------------------------------------------------------
688 722
689 for(G4int i = 0; i< numberOfVariables; ++i << 723
690 { << 724 for(G4int i = 0; i< numberOfVariables; i++)
691 yIn[i] = yInput[i]; 725 yIn[i] = yInput[i];
692 } << 726
693 <<
694 G4double tau0 = tau; 727 G4double tau0 = tau;
695 << 728 // Calculating the polynomials (coefficents for the respective stages) :
696 // Calculating the polynomials (coefficent << 729
697 // << 730 for(int i=1; i<=12; i++){ //Here i is NOT the coordinate no. , it's stage no.
698 for(auto i=1; i<=12; ++i) // i is NOT the <<
699 { <<
700 b[i] = 0; 731 b[i] = 0;
701 tau = 1.0; 732 tau = 1.0;
702 for(auto j=0; j<=5; ++j) << 733 for(int j=0; j<=5; j++){
703 { <<
704 b[i] += bi[i][j]*tau; 734 b[i] += bi[i][j]*tau;
705 tau*=tau0; 735 tau*=tau0;
706 } 736 }
707 } 737 }
708 738
709 // Calculating the interpolation at the fr << 739 // Calculating the interpolation at the fraction tau of the step using the polynomial
710 // the polynomial coefficients and the res << 740 // coefficients and the respective stages
711 // << 741
712 for(G4int i=0; i<numberOfVariables; ++i) / << 742 for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
713 { << 743 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
714 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[ << 744 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
715 b[4]*ak4[i << 745 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
716 b[7]*ak7[i <<
717 b[10]*ak10[i] 746 b[10]*ak10[i] + b[11]*ak11[i] + b[12]*ak12[i]);
718 } 747 }
719 } 748 }
>> 749
>> 750 //-----Verified--------- - hackabot
720 751