| 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 // G4FSALBogackiShampine45 implementation 26 // G4FSALBogackiShampine45 implementation
27 // 27 //
28 // The Butcher table of the Bogacki-Shampine-8 28 // The Butcher table of the Bogacki-Shampine-8-4-5 method is as follows:
29 // 29 //
30 // 0 | 30 // 0 |
31 // 1/6 | 1/6 31 // 1/6 | 1/6
32 // 2/9 | 2/27 4/27 32 // 2/9 | 2/27 4/27
33 // 3/7 | 183/1372 -162/343 1053/1372 33 // 3/7 | 183/1372 -162/343 1053/1372
34 // 2/3 | 68/297 -4/11 42/143 196 34 // 2/3 | 68/297 -4/11 42/143 1960/3861
35 // 3/4 | 597/22528 81/352 63099/5857 35 // 3/4 | 597/22528 81/352 63099/585728 58653/366080 4617/20480
36 // 1 | 174197/959244 -30942/79937 8152137/19 36 // 1 | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219
37 // 1 | 587/8064 0 4440339/15 37 // 1 | 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080
38 // ------------------------------------------- 38 // -------------------------------------------------------------------------------------------------------------------
39 // 587/8064 0 4440339/15 39 // 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 0
40 // 2479/34992 0 123/416 40 // 2479/34992 0 123/416 612941/3411720 43/1440 2272/6561 79937/1113912 3293/556956
41 // 41 //
42 // Created: Somnath Banerjee, Google Summer of 42 // Created: Somnath Banerjee, Google Summer of Code 2015, 26 May 2015
43 // Supervision: John Apostolakis, CERN 43 // Supervision: John Apostolakis, CERN
44 // ------------------------------------------- 44 // --------------------------------------------------------------------
45 45
46 // Plan is that this source file / class will 46 // Plan is that this source file / class will be merged with the updated
47 // BogackiShampine45 class, which contains imp 47 // BogackiShampine45 class, which contains improvements (May 2016)
48 48
49 #include <cassert> 49 #include <cassert>
50 50
51 #include "G4FSALBogackiShampine45.hh" 51 #include "G4FSALBogackiShampine45.hh"
52 #include "G4LineSection.hh" 52 #include "G4LineSection.hh"
53 53
54 G4bool G4FSALBogackiShampine45::fPreparedCon 54 G4bool G4FSALBogackiShampine45::fPreparedConstants = false;
55 G4double G4FSALBogackiShampine45::bi[12][7]; 55 G4double G4FSALBogackiShampine45::bi[12][7];
56 56
57 // Constructor 57 // Constructor
58 // 58 //
59 G4FSALBogackiShampine45::G4FSALBogackiShampine 59 G4FSALBogackiShampine45::G4FSALBogackiShampine45(G4EquationOfMotion* EqRhs,
60 60 G4int noIntegrationVariables,
61 61 G4bool primary)
62 : G4VFSALIntegrationStepper(EqRhs, noIntegr 62 : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
63 { 63 {
64 const G4int numberOfVariables = noIntegrat 64 const G4int numberOfVariables = noIntegrationVariables;
65 65
66 // New Chunk of memory being created for u 66 // New Chunk of memory being created for use by the stepper
67 67
68 // aki - for storing intermediate RHS 68 // aki - for storing intermediate RHS
69 // 69 //
70 ak2 = new G4double[numberOfVariables]; 70 ak2 = new G4double[numberOfVariables];
71 ak3 = new G4double[numberOfVariables]; 71 ak3 = new G4double[numberOfVariables];
72 ak4 = new G4double[numberOfVariables]; 72 ak4 = new G4double[numberOfVariables];
73 ak5 = new G4double[numberOfVariables]; 73 ak5 = new G4double[numberOfVariables];
74 ak6 = new G4double[numberOfVariables]; 74 ak6 = new G4double[numberOfVariables];
75 ak7 = new G4double[numberOfVariables]; 75 ak7 = new G4double[numberOfVariables];
76 ak8 = new G4double[numberOfVariables]; 76 ak8 = new G4double[numberOfVariables];
77 77
78 ak9 = new G4double[numberOfVariables]; 78 ak9 = new G4double[numberOfVariables];
79 ak10 = new G4double[numberOfVariables]; 79 ak10 = new G4double[numberOfVariables];
80 ak11 = new G4double[numberOfVariables]; 80 ak11 = new G4double[numberOfVariables];
81 DyDx = new G4double[numberOfVariables]; 81 DyDx = new G4double[numberOfVariables];
82 82
83 assert ( GetNumberOfStateVariables() >= 8 83 assert ( GetNumberOfStateVariables() >= 8 );
84 const G4int numStateVars = std::max(noInte 84 const G4int numStateVars = std::max(noIntegrationVariables,
85 GetNum 85 GetNumberOfStateVariables() );
86 86
87 // Must ensure space extra 'state' variabl 87 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
88 // 88 //
89 yTemp = new G4double[numStateVars]; 89 yTemp = new G4double[numStateVars];
90 yIn = new G4double[numStateVars] ; 90 yIn = new G4double[numStateVars] ;
91 91
92 fLastInitialVector = new G4double[numState 92 fLastInitialVector = new G4double[numStateVars] ;
93 fLastFinalVector = new G4double[numStateVa 93 fLastFinalVector = new G4double[numStateVars] ;
94 fLastDyDx = new G4double[numberOfVariables 94 fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
95 95
96 fMidVector = new G4double[numStateVars]; 96 fMidVector = new G4double[numStateVars];
97 fMidError = new G4double[numStateVars]; 97 fMidError = new G4double[numStateVars];
98 98
99 pseudoDydx_for_DistChord = new G4double[nu 99 pseudoDydx_for_DistChord = new G4double[numberOfVariables];
100 100
101 fMidVector = new G4double[numberOfVariable 101 fMidVector = new G4double[numberOfVariables];
102 fMidError = new G4double[numberOfVariable 102 fMidError = new G4double[numberOfVariables];
103 if( primary ) 103 if( primary )
104 { 104 {
105 fAuxStepper = new G4FSALBogackiShampine4 105 fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
106 106 !primary);
107 } 107 }
108 if( !fPreparedConstants ) 108 if( !fPreparedConstants )
109 { 109 {
110 PrepareConstants(); 110 PrepareConstants();
111 } 111 }
112 } 112 }
113 113
114 // Destructor 114 // Destructor
115 // 115 //
116 G4FSALBogackiShampine45::~G4FSALBogackiShampin 116 G4FSALBogackiShampine45::~G4FSALBogackiShampine45()
117 { 117 {
118 // Clear all previously allocated memory f 118 // Clear all previously allocated memory for stepper and DistChord
119 119
120 delete [] ak2; 120 delete [] ak2;
121 delete [] ak3; 121 delete [] ak3;
122 delete [] ak4; 122 delete [] ak4;
123 delete [] ak5; 123 delete [] ak5;
124 delete [] ak6; 124 delete [] ak6;
125 delete [] ak7; 125 delete [] ak7;
126 delete [] ak8; 126 delete [] ak8;
127 delete [] ak9; 127 delete [] ak9;
128 delete [] ak10; 128 delete [] ak10;
129 delete [] ak11; 129 delete [] ak11;
130 delete [] DyDx; 130 delete [] DyDx;
131 delete [] yTemp; 131 delete [] yTemp;
132 delete [] yIn; 132 delete [] yIn;
133 133
134 delete [] fLastInitialVector; 134 delete [] fLastInitialVector;
135 delete [] fLastFinalVector; 135 delete [] fLastFinalVector;
136 delete [] fLastDyDx; 136 delete [] fLastDyDx;
137 delete [] fMidVector; 137 delete [] fMidVector;
138 delete [] fMidError; 138 delete [] fMidError;
139 139
140 delete fAuxStepper; 140 delete fAuxStepper;
141 141
142 delete [] pseudoDydx_for_DistChord; 142 delete [] pseudoDydx_for_DistChord;
143 } 143 }
144 144
145 // Stepper 145 // Stepper
146 // 146 //
147 // Passing in the value of yInput[],the first 147 // Passing in the value of yInput[],the first time dydx[] and Step length
148 // Giving back yOut and yErr arrays for output 148 // Giving back yOut and yErr arrays for output and error respectively
149 // 149 //
150 void G4FSALBogackiShampine45::Stepper(const G4 150 void G4FSALBogackiShampine45::Stepper(const G4double yInput[],
151 const G4 151 const G4double dydx[],
152 G4 152 G4double Step,
153 G4 153 G4double yOut[],
154 G4 154 G4double yErr[],
155 G4 155 G4double nextDydx[])
156 { 156 {
157 G4int i; 157 G4int i;
158 158
159 // The various constants defined on the ba 159 // The various constants defined on the basis of butcher tableu
160 160
161 const G4double b21 = 1.0/6.0 , 161 const G4double b21 = 1.0/6.0 ,
162 b31 = 2.0/27.0 , b32 = 4.0/ 162 b31 = 2.0/27.0 , b32 = 4.0/27.0,
163 163
164 b41 = 183.0/1372.0 , b42 = 164 b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
165 165
166 b51 = 68.0/297.0, b52 = -4. 166 b51 = 68.0/297.0, b52 = -4.0/11.0,
167 b53 = 42.0/143.0, b54 = 196 167 b53 = 42.0/143.0, b54 = 1960.0/3861.0,
168 168
169 b61 = 597.0/22528.0, b62 = 169 b61 = 597.0/22528.0, b62 = 81.0/352.0,
170 b63 = 63099.0/585728.0, b64 170 b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
171 b65 = 4617.0/20480.0, 171 b65 = 4617.0/20480.0,
172 172
173 b71 = 174197.0/959244.0, b7 173 b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
174 b73 = 8152137.0/19744439.0, 174 b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
175 b75 = -29421.0/29068.0, b7 175 b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
176 176
177 b81 = 587.0/8064.0, b82 = 177 b81 = 587.0/8064.0, b82 = 0.0,
178 b83 = 4440339.0/15491840.0, 178 b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
179 b85 = 387.0/44800.0, b86 = 179 b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
180 b87 = 7267.0/94080.0, 180 b87 = 7267.0/94080.0,
181 181
182 182
183 // c1 = 2479.0/34992.0, 183 // c1 = 2479.0/34992.0,
184 // c2 = 0.0, 184 // c2 = 0.0,
185 // c3 = 123.0/416.0, 185 // c3 = 123.0/416.0,
186 // c4 = 612941.0/3411720.0, 186 // c4 = 612941.0/3411720.0,
187 // c5 = 43.0/1440.0, 187 // c5 = 43.0/1440.0,
188 // c6 = 2272.0/6561.0, 188 // c6 = 2272.0/6561.0,
189 // c7 = 79937.0/1113912.0, 189 // c7 = 79937.0/1113912.0,
190 // c8 = 3293.0/556956.0, 190 // c8 = 3293.0/556956.0,
191 191
192 // For the embedded higher order method on 192 // For the embedded higher order method only the difference of values
193 // taken and is used directly later instea 193 // taken and is used directly later instead of defining the last row
194 // of butcher table in a separate set of v 194 // of butcher table in a separate set of variables and taking the
195 // difference there 195 // difference there
196 196
197 dc1 = b81 - 2479.0/34992.0 197 dc1 = b81 - 2479.0/34992.0 ,
198 dc2 = 0.0, 198 dc2 = 0.0,
199 dc3 = b83 - 123.0/416.0 , 199 dc3 = b83 - 123.0/416.0 ,
200 dc4 = b84 - 612941.0/341172 200 dc4 = b84 - 612941.0/3411720.0,
201 dc5 = b85 - 43.0/1440.0, 201 dc5 = b85 - 43.0/1440.0,
202 dc6 = b86 - 2272.0/6561.0, 202 dc6 = b86 - 2272.0/6561.0,
203 dc7 = b87 - 79937.0/1113912 203 dc7 = b87 - 79937.0/1113912.0,
204 dc8 = -3293.0/556956.0; / 204 dc8 = -3293.0/556956.0; // end of declaration
205 205
206 const G4int numberOfVariables = GetNumberO 206 const G4int numberOfVariables = GetNumberOfVariables();
207 207
208 // The number of variables to be integrate 208 // The number of variables to be integrated over
209 // 209 //
210 yOut[7] = yTemp[7] = yIn[7]; 210 yOut[7] = yTemp[7] = yIn[7];
211 211
212 // Saving yInput because yInput and yOut 212 // Saving yInput because yInput and yOut can be aliases for same array
213 // 213 //
214 for(i=0; i<numberOfVariables; ++i) 214 for(i=0; i<numberOfVariables; ++i)
215 { 215 {
216 yIn[i]=yInput[i]; 216 yIn[i]=yInput[i];
217 DyDx[i] = dydx[i]; 217 DyDx[i] = dydx[i];
218 } 218 }
219 // RightHandSide(yIn, dydx) ; // 1st Ste 219 // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
220 220
221 for(i=0; i<numberOfVariables; ++i) 221 for(i=0; i<numberOfVariables; ++i)
222 { 222 {
223 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 223 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
224 } 224 }
225 RightHandSide(yTemp, ak2) ; / 225 RightHandSide(yTemp, ak2) ; // 2nd Step
226 226
227 for(i=0; i<numberOfVariables; ++i) 227 for(i=0; i<numberOfVariables; ++i)
228 { 228 {
229 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] 229 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
230 } 230 }
231 RightHandSide(yTemp, ak3) ; / 231 RightHandSide(yTemp, ak3) ; // 3rd Step
232 232
233 for(i=0; i<numberOfVariables; ++i) 233 for(i=0; i<numberOfVariables; ++i)
234 { 234 {
235 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] 235 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
236 } 236 }
237 RightHandSide(yTemp, ak4) ; / 237 RightHandSide(yTemp, ak4) ; // 4th Step
238 238
239 for(i=0; i<numberOfVariables; ++i) 239 for(i=0; i<numberOfVariables; ++i)
240 { 240 {
241 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] 241 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
242 b54*ak4[i]) 242 b54*ak4[i]) ;
243 } 243 }
244 RightHandSide(yTemp, ak5) ; / 244 RightHandSide(yTemp, ak5) ; // 5th Step
245 245
246 for(i=0; i<numberOfVariables; ++i) 246 for(i=0; i<numberOfVariables; ++i)
247 { 247 {
248 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] 248 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
249 b64*ak4[i] + 249 b64*ak4[i] + b65*ak5[i]) ;
250 } 250 }
251 RightHandSide(yTemp, ak6) ; / 251 RightHandSide(yTemp, ak6) ; // 6th Step
252 252
253 for(i=0; i<numberOfVariables; ++i) 253 for(i=0; i<numberOfVariables; ++i)
254 { 254 {
255 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] 255 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
256 b74*ak4[i] + 256 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
257 } 257 }
258 RightHandSide(yTemp, ak7); / 258 RightHandSide(yTemp, ak7); // 7th Step
259 259
260 for(i=0; i<numberOfVariables; ++i) 260 for(i=0; i<numberOfVariables; ++i)
261 { 261 {
262 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + 262 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
263 b84*ak4[i] + 263 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
264 b87*ak7[i]); 264 b87*ak7[i]);
265 } 265 }
266 RightHandSide(yOut, ak8); / 266 RightHandSide(yOut, ak8); // 8th Step - Final one Using FSAL
267 267
268 268
269 for(i=0; i<numberOfVariables; ++i) 269 for(i=0; i<numberOfVariables; ++i)
270 { 270 {
271 271
272 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[ 272 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
273 dc5*ak5[i] + dc6*ak6[i 273 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
274 274
275 275
276 // FSAL stepper : Must pass the last D 276 // FSAL stepper : Must pass the last DyDx for the next step, here ak8
277 // 277 //
278 nextDydx[i] = ak8[i]; 278 nextDydx[i] = ak8[i];
279 279
280 // Store Input and Final values, for p 280 // Store Input and Final values, for possible use in calculating chord
281 // 281 //
282 fLastInitialVector[i] = yIn[i] ; 282 fLastInitialVector[i] = yIn[i] ;
283 fLastFinalVector[i] = yOut[i]; 283 fLastFinalVector[i] = yOut[i];
284 fLastDyDx[i] = DyDx[i]; 284 fLastDyDx[i] = DyDx[i];
285 } 285 }
286 fLastStepLength = Step; 286 fLastStepLength = Step;
287 287
288 return; 288 return;
289 } 289 }
290 290
291 // DistChord 291 // DistChord
292 // 292 //
293 G4double G4FSALBogackiShampine45::DistChord() 293 G4double G4FSALBogackiShampine45::DistChord() const
294 { 294 {
295 G4double distLine, distChord; 295 G4double distLine, distChord;
296 G4ThreeVector initialPoint, finalPoint, mi 296 G4ThreeVector initialPoint, finalPoint, midPoint;
297 297
298 // Store last initial and final points 298 // Store last initial and final points
299 // (they will be overwritten in self-Stepp 299 // (they will be overwritten in self-Stepper call!)
300 // 300 //
301 initialPoint = G4ThreeVector( fLastInitial 301 initialPoint = G4ThreeVector( fLastInitialVector[0],
302 fLastInitialV 302 fLastInitialVector[1], fLastInitialVector[2]);
303 finalPoint = G4ThreeVector( fLastFinalVe 303 finalPoint = G4ThreeVector( fLastFinalVector[0],
304 fLastFinalVec 304 fLastFinalVector[1], fLastFinalVector[2]);
305 305
306 // Do half a step using StepNoErr 306 // Do half a step using StepNoErr
307 307
308 fAuxStepper->Stepper( fLastInitialVector, 308 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
309 fMidVector, fMidErro 309 fMidVector, fMidError, pseudoDydx_for_DistChord );
310 310
311 midPoint = G4ThreeVector( fMidVector[0], f 311 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
312 312
313 // Use stored values of Initial and Endpoi 313 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
314 // distance of Chord 314 // distance of Chord
315 // 315 //
316 if (initialPoint != finalPoint) 316 if (initialPoint != finalPoint)
317 { 317 {
318 distLine = G4LineSection::Distline(mid 318 distLine = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
319 distChord = distLine; 319 distChord = distLine;
320 } 320 }
321 else 321 else
322 { 322 {
323 distChord = (midPoint-initialPoint).ma 323 distChord = (midPoint-initialPoint).mag();
324 } 324 }
325 return distChord; 325 return distChord;
326 } 326 }
327 327
328 // PrepareConstants 328 // PrepareConstants
329 // 329 //
330 void G4FSALBogackiShampine45::PrepareConstants 330 void G4FSALBogackiShampine45::PrepareConstants()
331 { 331 {
332 // -------------------------------------- 332 // --------------------------------------------------------
333 // COEFFICIENTS FOR INTERPOLANT bi WITH 333 // COEFFICIENTS FOR INTERPOLANT bi WITH 11 STAGES
334 // -------------------------------------- 334 // --------------------------------------------------------
335 335
336 // Initialise all values of G4double bi[12 336 // Initialise all values of G4double bi[12][7]
337 // 337 //
338 for(auto i=1; i<12; ++i) 338 for(auto i=1; i<12; ++i)
339 { 339 {
340 for(auto j=1; j<7; ++j) 340 for(auto j=1; j<7; ++j)
341 { 341 {
342 bi[i][j] = 0.0 ; 342 bi[i][j] = 0.0 ;
343 } 343 }
344 } 344 }
345 345
346 bi[1][6] = -12134338393.0/1050809760.0 , 346 bi[1][6] = -12134338393.0/1050809760.0 ,
347 bi[1][5] = -1620741229.0/50038560.0 , 347 bi[1][5] = -1620741229.0/50038560.0 ,
348 bi[1][4] = -2048058893.0/59875200.0 , 348 bi[1][4] = -2048058893.0/59875200.0 ,
349 bi[1][3] = -87098480009.0/5254048800.0 , 349 bi[1][3] = -87098480009.0/5254048800.0 ,
350 bi[1][2] = -11513270273.0/3502699200.0 , 350 bi[1][2] = -11513270273.0/3502699200.0 ,
351 // 351 //
352 bi[3][6] = -33197340367.0/1218433216.0 , 352 bi[3][6] = -33197340367.0/1218433216.0 ,
353 bi[3][5] = -539868024987.0/6092166080.0 , 353 bi[3][5] = -539868024987.0/6092166080.0 ,
354 bi[3][4] = -39991188681.0/374902528.0 , 354 bi[3][4] = -39991188681.0/374902528.0 ,
355 bi[3][3] = -69509738227.0/1218433216.0 , 355 bi[3][3] = -69509738227.0/1218433216.0 ,
356 bi[3][2] = -29327744613.0/2436866432.0 , 356 bi[3][2] = -29327744613.0/2436866432.0 ,
357 // 357 //
358 bi[4][6] = -284800997201.0/19905339168.0 , 358 bi[4][6] = -284800997201.0/19905339168.0 ,
359 bi[4][5] = -7896875450471.0/165877826400.0 359 bi[4][5] = -7896875450471.0/165877826400.0 ,
360 bi[4][4] = -333945812879.0/5671036800.0 , 360 bi[4][4] = -333945812879.0/5671036800.0 ,
361 bi[4][3] = -16209923456237.0/497633479200. 361 bi[4][3] = -16209923456237.0/497633479200.0 ,
362 bi[4][2] = -2382590741699.0/331755652800.0 362 bi[4][2] = -2382590741699.0/331755652800.0 ,
363 // 363 //
364 bi[5][6] = -540919.0/741312.0 , 364 bi[5][6] = -540919.0/741312.0 ,
365 bi[5][5] = -103626067.0/43243200.0 , 365 bi[5][5] = -103626067.0/43243200.0 ,
366 bi[5][4] = -633779.0/211200.0 , 366 bi[5][4] = -633779.0/211200.0 ,
367 bi[5][3] = -32406787.0/18532800.0 , 367 bi[5][3] = -32406787.0/18532800.0 ,
368 bi[5][2] = -36591193.0/86486400.0 , 368 bi[5][2] = -36591193.0/86486400.0 ,
369 // 369 //
370 bi[6][6] = 7157998304.0/374350977.0 , 370 bi[6][6] = 7157998304.0/374350977.0 ,
371 bi[6][5] = 30405842464.0/623918295.0 , 371 bi[6][5] = 30405842464.0/623918295.0 ,
372 bi[6][4] = 183022264.0/5332635.0 , 372 bi[6][4] = 183022264.0/5332635.0 ,
373 bi[6][3] = -3357024032.0/1871754885.0 , 373 bi[6][3] = -3357024032.0/1871754885.0 ,
374 bi[6][2] = -611586736.0/89131185.0 , 374 bi[6][2] = -611586736.0/89131185.0 ,
375 // 375 //
376 bi[7][6] = -138073.0/9408.0 , 376 bi[7][6] = -138073.0/9408.0 ,
377 bi[7][5] = -719433.0/15680.0 , 377 bi[7][5] = -719433.0/15680.0 ,
378 bi[7][4] = -1620541.0/31360.0 , 378 bi[7][4] = -1620541.0/31360.0 ,
379 bi[7][3] = -385151.0/15680.0 , 379 bi[7][3] = -385151.0/15680.0 ,
380 bi[7][2] = -65403.0/15680.0 , 380 bi[7][2] = -65403.0/15680.0 ,
381 // 381 //
382 bi[8][6] = 1245.0/64.0 , 382 bi[8][6] = 1245.0/64.0 ,
383 bi[8][5] = 3991.0/64.0 , 383 bi[8][5] = 3991.0/64.0 ,
384 bi[8][4] = 4715.0/64.0 , 384 bi[8][4] = 4715.0/64.0 ,
385 bi[8][3] = 2501.0/64.0 , 385 bi[8][3] = 2501.0/64.0 ,
386 bi[8][2] = 149.0/16.0 , 386 bi[8][2] = 149.0/16.0 ,
387 bi[8][1] = 1.0 , 387 bi[8][1] = 1.0 ,
388 // 388 //
389 bi[9][6] = 55.0/3.0 , 389 bi[9][6] = 55.0/3.0 ,
390 bi[9][5] = 71.0 , 390 bi[9][5] = 71.0 ,
391 bi[9][4] = 103.0 , 391 bi[9][4] = 103.0 ,
392 bi[9][3] = 199.0/3.0 , 392 bi[9][3] = 199.0/3.0 ,
393 bi[9][2] = 16.0 , 393 bi[9][2] = 16.0 ,
394 // 394 //
395 bi[10][6] = -1774004627.0/75810735.0 , 395 bi[10][6] = -1774004627.0/75810735.0 ,
396 bi[10][5] = -1774004627.0/25270245.0 , 396 bi[10][5] = -1774004627.0/25270245.0 ,
397 bi[10][4] = -26477681.0/359975.0 , 397 bi[10][4] = -26477681.0/359975.0 ,
398 bi[10][3] = -11411880511.0/379053675.0 , 398 bi[10][3] = -11411880511.0/379053675.0 ,
399 bi[10][2] = -423642896.0/126351225.0 , 399 bi[10][2] = -423642896.0/126351225.0 ,
400 // 400 //
401 bi[11][6] = 35.0 , 401 bi[11][6] = 35.0 ,
402 bi[11][5] = 105.0 , 402 bi[11][5] = 105.0 ,
403 bi[11][4] = 117.0 , 403 bi[11][4] = 117.0 ,
404 bi[11][3] = 59.0 , 404 bi[11][3] = 59.0 ,
405 bi[11][2] = 12.0 ; 405 bi[11][2] = 12.0 ;
406 } 406 }
407 407
408 // ------------------------------------------- 408 // ---------------------------------------------------------------------------------------
409 409
410 void G4FSALBogackiShampine45::interpolate( con 410 void G4FSALBogackiShampine45::interpolate( const G4double yInput[],
411 con 411 const G4double dydx[],
412 412 G4double yOut[],
413 413 G4double Step,
414 414 G4double tau )
415 { 415 {
416 const G4double a91 = 455.0/6144.0 , 416 const G4double a91 = 455.0/6144.0 ,
417 a92 = 0.0 , 417 a92 = 0.0 ,
418 a93 = 10256301.0/35409920.0 418 a93 = 10256301.0/35409920.0 ,
419 a94 = 2307361.0/17971200.0 419 a94 = 2307361.0/17971200.0 ,
420 a95 = -387.0/102400.0 , 420 a95 = -387.0/102400.0 ,
421 a96 = 73.0/5130.0 , 421 a96 = 73.0/5130.0 ,
422 a97 = -7267.0/215040.0 , 422 a97 = -7267.0/215040.0 ,
423 a98 = 1.0/32.0 , 423 a98 = 1.0/32.0 ,
424 424
425 a101 = -837888343715.0/1317 425 a101 = -837888343715.0/13176988637184.0 ,
426 a102 = 30409415.0/52955362. 426 a102 = 30409415.0/52955362.0 ,
427 a103 = -48321525963.0/75916 427 a103 = -48321525963.0/759168069632.0 ,
428 a104 = 8530738453321.0/1976 428 a104 = 8530738453321.0/197654829557760.0 ,
429 a105 = 1361640523001.0/1626 429 a105 = 1361640523001.0/1626788720640.0 ,
430 a106 = -13143060689.0/38604 430 a106 = -13143060689.0/38604458898.0 ,
431 a107 = 18700221969.0/379584 431 a107 = 18700221969.0/379584034816.0 ,
432 a108 = -5831595.0/847285792 432 a108 = -5831595.0/847285792.0 ,
433 a109 = -5183640.0/26477681. 433 a109 = -5183640.0/26477681.0 ,
434 434
435 a111 = 98719073263.0/155196 435 a111 = 98719073263.0/1551965184000.0 ,
436 a112 = 1307.0/123552.0 , 436 a112 = 1307.0/123552.0 ,
437 a113 = 4632066559387.0/7018 437 a113 = 4632066559387.0/70181753241600.0 ,
438 a114 = 7828594302389.0/3821 438 a114 = 7828594302389.0/382182512025600.0 ,
439 a115 = 40763687.0/110702592 439 a115 = 40763687.0/11070259200.0 ,
440 a116 = 34872732407.0/224610 440 a116 = 34872732407.0/224610586200.0 ,
441 a117 = -2561897.0/30105600. 441 a117 = -2561897.0/30105600.0 ,
442 a118 = 1.0/10.0 , 442 a118 = 1.0/10.0 ,
443 a119 = -1.0/10.0 , 443 a119 = -1.0/10.0 ,
444 a1110 = -1403317093.0/11371 444 a1110 = -1403317093.0/11371610250.0 ;
445 445
446 const G4int numberOfVariables = GetNumberO 446 const G4int numberOfVariables = GetNumberOfVariables();
447 447
448 // Saving yInput because yInput and yOut c 448 // Saving yInput because yInput and yOut can be aliases for same array
449 // 449 //
450 for(G4int i=0; i<numberOfVariables; ++i) 450 for(G4int i=0; i<numberOfVariables; ++i)
451 { 451 {
452 yIn[i]=yInput[i]; 452 yIn[i]=yInput[i];
453 } 453 }
454 454
455 // The number of variables to be integrate 455 // The number of variables to be integrated over
456 // 456 //
457 yOut[7] = yTemp[7] = yIn[7]; 457 yOut[7] = yTemp[7] = yIn[7];
458 458
459 // Calculating extra stages 459 // Calculating extra stages
460 // 460 //
461 for(G4int i=0; i<numberOfVariables; ++i) 461 for(G4int i=0; i<numberOfVariables; ++i)
462 { 462 {
463 yTemp[i] = yIn[i] + Step*(a91*dydx[i] 463 yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
464 a94*ak4[i] + 464 a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
465 a97*ak7[i] + 465 a97*ak7[i] + a98*ak8[i] );
466 } 466 }
467 467
468 RightHandSide(yTemp, ak9); 468 RightHandSide(yTemp, ak9);
469 469
470 for(G4int i=0; i<numberOfVariables; ++i) 470 for(G4int i=0; i<numberOfVariables; ++i)
471 { 471 {
472 yTemp[i] = yIn[i] + Step*(a101*dydx[i] 472 yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
473 a104*ak4[i] 473 a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
474 a107*ak7[i] 474 a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
475 } 475 }
476 476
477 RightHandSide(yTemp, ak10); 477 RightHandSide(yTemp, ak10);
478 478
479 for(G4int i=0; i<numberOfVariables; ++i) 479 for(G4int i=0; i<numberOfVariables; ++i)
480 { 480 {
481 yTemp[i] = yIn[i] + Step*(a111*dydx[i] 481 yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
482 a114*ak4[i] 482 a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
483 a117*ak7[i] 483 a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
484 a1110*ak10[i 484 a1110*ak10[i] );
485 } 485 }
486 486
487 RightHandSide(yTemp, ak11); 487 RightHandSide(yTemp, ak11);
488 488
489 G4double tau0 = tau; 489 G4double tau0 = tau;
490 490
491 // Calculating the polynomials 491 // Calculating the polynomials
492 // 492 //
493 for(auto i=1; i<=11; ++i) // i is NOT the 493 for(auto i=1; i<=11; ++i) // i is NOT the coordinate no., it's stage no.
494 { 494 {
495 b[i] = 0.0; 495 b[i] = 0.0;
496 tau = tau0; 496 tau = tau0;
497 for(auto j=1; j<=6; ++j) 497 for(auto j=1; j<=6; ++j)
498 { 498 {
499 b[i] += bi[i][j]*tau; 499 b[i] += bi[i][j]*tau;
500 tau*=tau0; 500 tau*=tau0;
501 } 501 }
502 } 502 }
503 503
504 for(G4int i=0; i<numberOfVariables; ++i) 504 for(G4int i=0; i<numberOfVariables; ++i)
505 { 505 {
506 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] 506 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
507 b[4]*ak4[i] + 507 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
508 b[7]*ak7[i] + 508 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
509 b[10]*ak10[i] 509 b[10]*ak10[i] + b[11]*ak11[i] );
510 } 510 }
511 } 511 }
512 512