| 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 // G4FSALDormandPrince745 implementation << 26 // DormandPrince7 - 5(4) implementation by Somnath Banerjee
>> 27 // Supervision / code review: John Apostolakis
27 // 28 //
28 // The Butcher table of the FDormand-Prince-7- << 29 // Sponsored by Google in Google Summer of Code 2015.
29 // 30 //
30 // 0 | << 31 // First version: 25 May 2015
31 // 1/5 | 1/5 <<
32 // 3/10| 3/40 9/40 <<
33 // 4/5 | 44/45 56/15 32/9 <<
34 // 8/9 | 19372/6561 25360/2187 64448/6561 2 <<
35 // 1 | 9017/3168 355/33 46732/5247 4 <<
36 // 1 | 35/384 0 500/1113 1 <<
37 // ------------------------------------------ <<
38 // 35/384 0 500/1113 1 <<
39 // 5179/57600 0 7571/16695 3 <<
40 // 32 //
41 // Created: Somnath Banerjee, Google Summer of << 33 // G4FSALDormandPrince745.cc
42 // Supervision: John Apostolakis, CERN << 34 // Geant4
43 // ------------------------------------------- << 35 //
>> 36 //
>> 37 // This is the source file of G4FSALDormandPrince745 class containing the
>> 38 // definition of the stepper() method that evaluates one step in
>> 39 // field propagation.
>> 40 // The Butcher table of the FDormand-Prince-7-4-5 method is as follows :
>> 41 //
>> 42 // 0 |
>> 43 // 1/5 | 1/5
>> 44 // 3/10| 3/40 9/40
>> 45 // 4/5 | 44/45 −56/15 32/9
>> 46 // 8/9 | 19372/6561 −25360/2187 64448/6561 −212/729
>> 47 // 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656
>> 48 // 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84
>> 49 // ----------f--------------------------------------------------------------
>> 50 // 35/384 0 500/1113 125/192 −2187/6784 11/84 0
>> 51 // 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40
>> 52
44 53
45 #include "G4FSALDormandPrince745.hh" 54 #include "G4FSALDormandPrince745.hh"
46 #include "G4LineSection.hh" 55 #include "G4LineSection.hh"
47 #include <cmath> 56 #include <cmath>
48 57
49 // Constructor << 58 //Constructor
50 // << 59 G4FSALDormandPrince745::G4FSALDormandPrince745(G4EquationOfMotion *EqRhs,
51 G4FSALDormandPrince745::G4FSALDormandPrince745 << 60 G4int noIntegrationVariables,
52 << 61 G4bool primary)
53 << 62 : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables){
54 : G4VFSALIntegrationStepper(EqRhs, noIntegra << 63
55 { << 64 const G4int numberOfVariables = noIntegrationVariables;
56 const G4int numberOfVariables = noIntegratio << 65
57 << 66 //New Chunk of memory being created for use by the stepper
58 // New Chunk of memory being created for use << 67
59 << 68 //aki - for storing intermediate RHS
60 // aki - for storing intermediate RHS << 69 ak2 = new G4double[numberOfVariables];
61 // << 70 ak3 = new G4double[numberOfVariables];
62 ak2 = new G4double[numberOfVariables]; << 71 ak4 = new G4double[numberOfVariables];
63 ak3 = new G4double[numberOfVariables]; << 72 ak5 = new G4double[numberOfVariables];
64 ak4 = new G4double[numberOfVariables]; << 73 ak6 = new G4double[numberOfVariables];
65 ak5 = new G4double[numberOfVariables]; << 74 ak7 = new G4double[numberOfVariables];
66 ak6 = new G4double[numberOfVariables]; << 75
67 ak7 = new G4double[numberOfVariables]; << 76 yTemp = new G4double[numberOfVariables] ;
68 << 77 yIn = new G4double[numberOfVariables] ;
69 // Also always allocate arrays for interpola << 78
70 // << 79 pseudoDydx_for_DistChord = new G4double[numberOfVariables];
71 ak8 = new G4double[numberOfVariables]; << 80
72 ak9 = new G4double[numberOfVariables]; << 81 fLastInitialVector = new G4double[numberOfVariables] ;
73 << 82 fLastFinalVector = new G4double[numberOfVariables] ;
74 yTemp = new G4double[numberOfVariables] ; << 83 fLastDyDx = new G4double[numberOfVariables];
75 yIn = new G4double[numberOfVariables] ; << 84
76 << 85 fMidVector = new G4double[numberOfVariables];
77 pseudoDydx_for_DistChord = new G4double[numb << 86 fMidError = new G4double[numberOfVariables];
78 << 87 if( primary )
79 fInitialDyDx = new G4double[numberOfVariable << 88 {
80 fLastInitialVector = new G4double[numberOfVa << 89 fAuxStepper = new G4FSALDormandPrince745(EqRhs, numberOfVariables,
81 fLastFinalVector = new G4double[numberOfVari << 90 !primary);
82 fLastDyDx = new G4double[numberOfVariables]; << 91 }
83 <<
84 fMidVector = new G4double[numberOfVariables] <<
85 fMidError = new G4double[numberOfVariables] <<
86 <<
87 if( primary ) <<
88 { <<
89 fAuxStepper = new G4FSALDormandPrince745(E <<
90 } <<
91 } 92 }
92 93
93 // Destructor <<
94 // <<
95 G4FSALDormandPrince745::~G4FSALDormandPrince74 <<
96 { <<
97 // Clear all previously allocated memory f <<
98 94
99 delete [] ak2; ak2 = nullptr; << 95 //Destructor
100 delete [] ak3; ak3 = nullptr; << 96 G4FSALDormandPrince745::~G4FSALDormandPrince745(){
101 delete [] ak4; ak4 = nullptr; << 97 //clear all previously allocated memory for stepper and DistChord
102 delete [] ak5; ak5 = nullptr; << 98 delete[] ak2;
103 delete [] ak6; ak6 = nullptr; << 99 delete[] ak3;
104 delete [] ak7; ak7 = nullptr; << 100 delete[] ak4;
105 delete [] ak8; ak8 = nullptr; << 101 delete[] ak5;
106 delete [] ak9; ak9 = nullptr; << 102 delete[] ak6;
107 << 103 delete[] ak7;
108 delete [] yTemp; yTemp = nullptr; << 104
109 delete [] yIn; yIn = nullptr; << 105 delete[] yTemp;
110 << 106 delete[] yIn;
111 delete [] pseudoDydx_for_DistChord; pseud << 107
112 delete [] fInitialDyDx; fInit << 108 delete[] fLastInitialVector;
113 << 109 delete[] fLastFinalVector;
114 delete [] fLastInitialVector; fLastInit << 110 delete[] fLastDyDx;
115 delete [] fLastFinalVector; fLastFina << 111 delete[] fMidVector;
116 delete [] fLastDyDx; fLastDyDx << 112 delete[] fMidError;
117 delete [] fMidVector; fMidVecto << 113
118 delete [] fMidError; fMidError << 114 delete fAuxStepper;
>> 115
>> 116 delete[] pseudoDydx_for_DistChord;
119 117
120 delete fAuxStepper; fAuxStepper = nullptr; <<
121 } 118 }
122 119
123 // Stepper << 120
124 // << 121 //Stepper :
>> 122
125 // Passing in the value of yInput[],the first 123 // Passing in the value of yInput[],the first time dydx[] and Step length
126 // Giving back yOut and yErr arrays for output 124 // Giving back yOut and yErr arrays for output and error respectively
127 // << 125
128 void G4FSALDormandPrince745::Stepper(const G4d 126 void G4FSALDormandPrince745::Stepper(const G4double yInput[],
129 const G4d << 127 const G4double dydx[],
130 G4d << 128 G4double Step,
131 G4d << 129 G4double yOut[],
132 G4d << 130 G4double yErr[],
133 G4d << 131 G4double nextDydx[]
>> 132 )
134 { 133 {
135 G4int i; 134 G4int i;
136 135
137 // The various constants defined on the ba << 136 //The various constants defined on the basis of butcher tableu
138 << 137 const G4double //G4double - only once
139 const G4double b21 = 0.2 , << 138 b21 = 0.2 ,
140 b31 = 3.0/40.0, b32 = 9.0/4 <<
141 139
142 b41 = 44.0/45.0, b42 = -56. << 140 b31 = 3.0/40.0, b32 = 9.0/40.0 ,
143 141
144 b51 = 19372.0/6561.0, b52 = << 142 b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
145 b53 = 64448.0/6561.0, b54 = <<
146 143
147 b61 = 9017.0/3168.0 , b62 = << 144 b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
148 b63 = 46732.0/5247.0 , << 145 b54 = -212.0/729.0 ,
149 b65 = -5103.0/18656.0 , <<
150 146
151 b71 = 35.0/384.0, b72 = 0., << 147 b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
152 b73 = 500.0/1113.0, b74 = 1 << 148 b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
153 b75 = -2187.0/6784.0, b76 = << 149 b65 = -5103.0/18656.0 ,
154 150
155 // c1 = 35.0/384.0, c2 = .0, << 151 b71 = 35.0/384.0, b72 = 0.,
156 // c3 = 500.0/1113.0, c4 = 125 << 152 b73 = 500.0/1113.0, b74 = 125.0/192.0,
157 // c5 = -2187.0/6784.0, c6 = 1 << 153 b75 = -2187.0/6784.0, b76 = 11.0/84.0,
158 // c7 = 0, <<
159 154
160 dc1 = b71 - 5179.0/57600.0, << 155 // c1 = 35.0/384.0, c2 = .0,
161 dc2 = b72 - .0, << 156 // c3 = 500.0/1113.0, c4 = 125.0/192.0,
162 dc3 = b73 - 7571.0/16695.0, << 157 // c5 = -2187.0/6784.0, c6 = 11.0/84.0,
163 dc4 = b74 - 393.0/640.0, << 158 // c7 = 0,
164 dc5 = b75 + 92097.0/339200. << 159
165 dc6 = b76 - 187.0/2100.0, << 160 dc1 = b71 - 5179.0/57600.0,
166 dc7 = - 1.0/40.0 ; //end o << 161 dc2 = b72 - .0,
167 << 162 dc3 = b73 - 7571.0/16695.0,
168 const G4int numberOfVariables = GetNumberO << 163 dc4 = b74 - 393.0/640.0,
169 // The number of variables to be integra << 164 dc5 = b75 + 92097.0/339200.0,
170 << 165 dc6 = b76 - 187.0/2100.0,
171 // Saving yInput because yInput and yOut c << 166 dc7 = - 1.0/40.0 ; //end of declaration
172 // << 167
173 for(i=0; i<numberOfVariables; ++i) << 168
>> 169 const G4int numberOfVariables= this->GetNumberOfVariables();
>> 170 G4double *DyDx = new G4double[numberOfVariables];
>> 171
>> 172 // The number of variables to be integrated over
>> 173 yOut[7] = yTemp[7] = yIn[7];
>> 174 // Saving yInput because yInput and yOut can be aliases for same array
>> 175
>> 176 for(i=0;i<numberOfVariables;i++)
174 { 177 {
175 yIn[i] = yInput[i]; << 178 yIn[i]=yInput[i];
176 fInitialDyDx[i] = dydx[i]; << 179 DyDx[i] = dydx[i];
177 } 180 }
178 // Ensure that time is initialised - in ca <<
179 // <<
180 yOut[7] = yTemp[7] = yInput[7]; <<
181 // RightHandSide(yIn, DyDx) ; // 1st Ste <<
182 181
183 for(i=0; i<numberOfVariables; ++i) << 182
>> 183
>> 184 // RightHandSide(yIn, DyDx) ;
>> 185 // 1st Step - Not doing, getting passed
>> 186
>> 187 for(i=0;i<numberOfVariables;i++)
184 { 188 {
185 yTemp[i] = yIn[i] + b21*Step*fInitialD << 189 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
186 } 190 }
187 RightHandSide(yTemp, ak2) ; / 191 RightHandSide(yTemp, ak2) ; // 2nd Step
188 192
189 for(i=0; i<numberOfVariables; ++i) << 193 for(i=0;i<numberOfVariables;i++)
190 { 194 {
191 yTemp[i] = yIn[i] + Step*(b31*fInitial << 195 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
192 } 196 }
193 RightHandSide(yTemp, ak3) ; / 197 RightHandSide(yTemp, ak3) ; // 3rd Step
194 198
195 for(i=0; i<numberOfVariables; ++i) << 199 for(i=0;i<numberOfVariables;i++)
196 { 200 {
197 yTemp[i] = yIn[i] + Step*(b41*fInitial << 201 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
198 + b42*ak2[i] + <<
199 } 202 }
200 RightHandSide(yTemp, ak4) ; / 203 RightHandSide(yTemp, ak4) ; // 4th Step
201 204
202 for(i=0; i<numberOfVariables; ++i) << 205 for(i=0;i<numberOfVariables;i++)
203 { 206 {
204 yTemp[i] = yIn[i] + Step*(b51*fInitial << 207 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
205 + b52*ak2[i] + << 208 b54*ak4[i]) ;
206 } 209 }
207 RightHandSide(yTemp, ak5) ; / 210 RightHandSide(yTemp, ak5) ; // 5th Step
208 211
209 for(i=0; i<numberOfVariables; ++i) << 212 for(i=0;i<numberOfVariables;i++)
210 { 213 {
211 yTemp[i] = yIn[i] + Step*(b61*fInitial << 214 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
212 + b63*ak3[i] + << 215 b64*ak4[i] + b65*ak5[i]) ;
213 } 216 }
214 RightHandSide(yTemp, ak6) ; / 217 RightHandSide(yTemp, ak6) ; // 6th Step
215 218
216 for(i=0; i<numberOfVariables; ++i) << 219 for(i=0;i<numberOfVariables;i++)
217 { 220 {
218 yOut[i] = yIn[i] + Step*(b71*fInitialD << 221 yOut[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
219 + b74*ak4[i] + << 222 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
220 } 223 }
221 RightHandSide(yOut, ak7); // 224 RightHandSide(yOut, ak7); //7th and Final step
222 225
223 for(i=0; i<numberOfVariables; ++i) << 226 for(i=0;i<numberOfVariables;i++)
224 { 227 {
225 228
226 yErr[i] = Step*(dc1*fInitialDyDx[i] + << 229 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
227 + dc4*ak4[i] + dc5*ak5[i << 230 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
>> 231
228 232
229 // Store Input and Final values, for p 233 // Store Input and Final values, for possible use in calculating chord
230 // <<
231 fLastInitialVector[i] = yIn[i] ; 234 fLastInitialVector[i] = yIn[i] ;
232 fLastFinalVector[i] = yOut[i]; 235 fLastFinalVector[i] = yOut[i];
233 fLastDyDx[i] = fInitialDyDx[i << 236 fLastDyDx[i] = DyDx[i];
234 nextDydx[i] = ak7[i]; 237 nextDydx[i] = ak7[i];
>> 238
>> 239
235 } 240 }
>> 241
236 fLastStepLength = Step; 242 fLastStepLength = Step;
237 243
238 return ; 244 return ;
239 } 245 }
240 246
241 // DistChord << 247
242 // << 248 //The following has not been tested
243 G4double G4FSALDormandPrince745::DistChord() c << 249
>> 250 //The DistChord() function fot the class - must define it here.
>> 251 G4double G4FSALDormandPrince745::DistChord() const
244 { 252 {
245 G4double distLine, distChord; 253 G4double distLine, distChord;
246 G4ThreeVector initialPoint, finalPoint, mi 254 G4ThreeVector initialPoint, finalPoint, midPoint;
247 255
248 // Store last initial and final points << 256 // Store last initial and final points (they will be overwritten in self-Stepper call!)
249 // (they will be overwritten in self-Stepp << 257 initialPoint = G4ThreeVector( fLastInitialVector[0],
250 // <<
251 initialPoint = G4ThreeVector(fLastInitialV <<
252 fLastInitialV 258 fLastInitialVector[1], fLastInitialVector[2]);
253 finalPoint = G4ThreeVector(fLastFinalVec << 259 finalPoint = G4ThreeVector( fLastFinalVector[0],
254 fLastFinalVec 260 fLastFinalVector[1], fLastFinalVector[2]);
255 261
256 // Do half a step using StepNoErr 262 // Do half a step using StepNoErr
257 263
258 fAuxStepper->Stepper( fLastInitialVector, 264 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
259 fMidVector, fMidErro << 265 fMidVector, fMidError, pseudoDydx_for_DistChord );
260 266
261 midPoint = G4ThreeVector( fMidVector[0], f << 267 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
262 268
263 // Use stored values of Initial and Endpoi 269 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
264 // distance of Chord << 270 // distance of Chord
265 // << 271
>> 272
266 if (initialPoint != finalPoint) 273 if (initialPoint != finalPoint)
267 { 274 {
268 distLine = G4LineSection::Distline( m << 275 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
269 distChord = distLine; 276 distChord = distLine;
270 } 277 }
271 else 278 else
272 { 279 {
273 distChord = (midPoint-initialPoint).ma 280 distChord = (midPoint-initialPoint).mag();
274 } 281 }
275 return distChord; 282 return distChord;
276 } 283 }
277 284
278 // interpolate <<
279 // <<
280 void G4FSALDormandPrince745::interpolate( cons <<
281 cons <<
282 <<
283 <<
284 <<
285 { <<
286 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7 <<
287 285
288 const G4int numberOfVariables = GetNumberO << 286 void G4FSALDormandPrince745::interpolate( const G4double yInput[],
>> 287 const G4double dydx[],
>> 288 G4double yOut[],
>> 289 G4double Step,
>> 290 G4double tau){
>> 291
>> 292 G4double
>> 293 bf1, bf2, bf3, bf4, bf5, bf6, bf7;
>> 294
>> 295
>> 296
>> 297
>> 298
>> 299 const G4int numberOfVariables= this->GetNumberOfVariables();
289 300
290 G4double tau0 = tau; 301 G4double tau0 = tau;
291 302
292 for(G4int i=0;i<numberOfVariables; ++i) << 303 for(int i=0;i<numberOfVariables;i++)
293 { 304 {
294 yIn[i]=yInput[i]; 305 yIn[i]=yInput[i];
295 } 306 }
296 307
297 G4double tau_2 = tau0*tau0 , << 308 G4double
298 tau_3 = tau0*tau_2, << 309 tau_2 = tau0*tau0 ,
299 tau_4 = tau_2*tau_2; << 310 tau_3 = tau0*tau_2,
300 << 311 tau_4 = tau_2*tau_2;
301 bf1 = (157015080.0*tau_4 - 13107642775.0*t << 312
302 + 34969693132.0*tau_2- 32272833064.0* << 313 bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
303 / 11282082432.0; << 314 + 11282082432.0)/11282082432.0,
304 bf2 = 0.0; << 315 bf2 = 0.0 ,
305 bf3 = - 100.0*tau*(15701508.0*tau_3 - 9141 << 316 bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2 + 2074956840.0*tau
306 + 2074956840.0*tau - 1323 << 317 - 1323431896.0)/32700410799.0,
307 bf4 = 25.0*tau*(94209048.0*tau_3- 15184142 << 318 bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2+ 2460397220.0*tau - 889289856.0)/5641041216.0 ,
308 + 2460397220.0*tau - 8892898 << 319 bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2 + 687873124.0*tau - 259006536.0)/199316789632.0 ,
309 / 5641041216.0; << 320 bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2 + 946554244.0*tau - 361440756.0)/2467955532.0 ,
310 bf5 = -2187.0*tau*(52338360.0*tau_3 - 4518 << 321 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)/ 29380423.0 ;
311 + 687873124.0*tau - 25900 << 322
312 / 199316789632.0; << 323
313 bf6 = 11.0*tau*(106151040.0*tau_3- 661884 << 324 for( int i=0; i<numberOfVariables; i++){
314 + 946554244.0*tau - 3614407 << 325 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
315 / 2467955532.0; << 326 + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
316 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 8 <<
317 / 29380423.0; <<
318 <<
319 for(G4int i=0; i<numberOfVariables; ++i) <<
320 { <<
321 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] <<
322 + bf4*ak4[i] <<
323 + bf7*ak7[i] <<
324 } 327 }
>> 328
>> 329
>> 330
325 } 331 }
326 332
327 // SetupInterpolate <<
328 // <<
329 void G4FSALDormandPrince745::SetupInterpolate( 333 void G4FSALDormandPrince745::SetupInterpolate(const G4double yInput[],
330 << 334 const G4double dydx[],
331 << 335 const G4double Step ){
332 { <<
333 // Coefficients for the additional stages <<
334 // <<
335 G4double b81 = 6245.0/62208.0 , <<
336 b82 = 0.0 , <<
337 b83 = 8875.0/103032.0 , <<
338 b84 = -125.0/1728.0 , <<
339 b85 = 801.0/13568.0 , <<
340 b86 = -13519.0/368064.0 , <<
341 b87 = 11105.0/368064.0 , <<
342 <<
343 b91 = 632855.0/4478976.0 , <<
344 b92 = 0.0 , <<
345 b93 = 4146875.0/6491016.0 , <<
346 b94 = 5490625.0/14183424.0 , <<
347 b95 = -15975.0/108544.0 , <<
348 b96 = 8295925.0/220286304.0 , <<
349 b97 = -1779595.0/62938944.0 , <<
350 b98 = -805.0/4104.0 ; <<
351 336
352 const G4int numberOfVariables = GetNumberO << 337 //Coefficients for the additional stages :
>> 338 G4double
>> 339 b81 = 6245.0/62208.0 ,
>> 340 b82 = 0.0 ,
>> 341 b83 = 8875.0/103032.0 ,
>> 342 b84 = -125.0/1728.0 ,
>> 343 b85 = 801.0/13568.0 ,
>> 344 b86 = -13519.0/368064.0 ,
>> 345 b87 = 11105.0/368064.0 ,
>> 346
>> 347 b91 = 632855.0/4478976.0 ,
>> 348 b92 = 0.0 ,
>> 349 b93 = 4146875.0/6491016.0 ,
>> 350 b94 = 5490625.0/14183424.0 ,
>> 351 b95 = -15975.0/108544.0 ,
>> 352 b96 = 8295925.0/220286304.0 ,
>> 353 b97 = -1779595.0/62938944.0 ,
>> 354 b98 = -805.0/4104.0 ;
353 355
354 // Saving yInput because yInput and yOut c << 356 const G4int numberOfVariables= this->GetNumberOfVariables();
355 // << 357
356 for(G4int i=0; i<numberOfVariables; ++i) << 358 // Saving yInput because yInput and yOut can be aliases for same array
>> 359 for(int i=0;i<numberOfVariables;i++)
357 { 360 {
358 yIn[i] = yInput[i]; << 361 yIn[i]=yInput[i];
359 } 362 }
360 363
361 yTemp[7] = yIn[7]; << 364 yTemp[7] = yIn[7];
362 365
363 // Evaluate the extra stages << 366 ak8 = new G4double[numberOfVariables];
364 // << 367 ak9 = new G4double[numberOfVariables];
365 for(G4int i=0; i<numberOfVariables; ++i) << 368
>> 369 //Evaluate the extra stages :
>> 370 for(int i=0;i<numberOfVariables;i++)
366 { 371 {
367 yTemp[i] = yIn[i] + Step*( b81*dydx[i] 372 yTemp[i] = yIn[i] + Step*( b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
368 b84*ak4[i] 373 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
369 b87*ak7[i] 374 b87*ak7[i] );
370 } 375 }
371 RightHandSide( yTemp, ak8 ); << 376 RightHandSide( yTemp, ak8 ); //8th Stage
372 377
373 for(G4int i=0; i<numberOfVariables; ++i) << 378 for(int i=0;i<numberOfVariables;i++)
374 { 379 {
375 yTemp[i] = yIn[i] + Step * ( b91*dydx[ 380 yTemp[i] = yIn[i] + Step * ( b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
376 b94*ak4[i] 381 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
377 b97*ak7[i] 382 b97*ak7[i] + b98*ak8[i] );
378 } 383 }
379 RightHandSide( yTemp, ak9 ); << 384 RightHandSide( yTemp, ak9 ); //9th Stage
>> 385
>> 386
>> 387
380 } 388 }
381 389
382 // Interpolate <<
383 // <<
384 void G4FSALDormandPrince745::Interpolate( cons <<
385 cons <<
386 cons <<
387 <<
388 <<
389 { <<
390 // Define the coefficients for the polynom <<
391 390
>> 391
>> 392 void G4FSALDormandPrince745::Interpolate( const G4double yInput[],
>> 393 const G4double dydx[],
>> 394 const G4double Step,
>> 395 G4double yOut[],
>> 396 G4double tau ){
>> 397 //Define the coefficients for the polynomials
392 G4double bi[10][5], b[10]; 398 G4double bi[10][5], b[10];
393 G4int numberOfVariables = GetNumberOfVaria << 399 G4int numberOfVariables = this->GetNumberOfVariables();
394 400
395 // COEFFICIENTS OF bi[1] 401 // COEFFICIENTS OF bi[1]
396 bi[1][0] = 1.0 , 402 bi[1][0] = 1.0 ,
397 bi[1][1] = -38039.0/7040.0 , 403 bi[1][1] = -38039.0/7040.0 ,
398 bi[1][2] = 125923.0/10560.0 , 404 bi[1][2] = 125923.0/10560.0 ,
399 bi[1][3] = -19683.0/1760.0 , 405 bi[1][3] = -19683.0/1760.0 ,
400 bi[1][4] = 3303.0/880.0 , 406 bi[1][4] = 3303.0/880.0 ,
401 // -------------------------------------- 407 // --------------------------------------------------------
402 // 408 //
403 // COEFFICIENTS OF bi[2] 409 // COEFFICIENTS OF bi[2]
404 bi[2][0] = 0.0 , 410 bi[2][0] = 0.0 ,
405 bi[2][1] = 0.0 , 411 bi[2][1] = 0.0 ,
406 bi[2][2] = 0.0 , 412 bi[2][2] = 0.0 ,
407 bi[2][3] = 0.0 , 413 bi[2][3] = 0.0 ,
408 bi[2][4] = 0.0 , 414 bi[2][4] = 0.0 ,
409 // -------------------------------------- 415 // --------------------------------------------------------
410 // 416 //
411 // COEFFICIENTS OF bi[3] 417 // COEFFICIENTS OF bi[3]
412 bi[3][0] = 0.0 , 418 bi[3][0] = 0.0 ,
413 bi[3][1] = -12500.0/4081.0 , 419 bi[3][1] = -12500.0/4081.0 ,
414 bi[3][2] = 205000.0/12243.0 , 420 bi[3][2] = 205000.0/12243.0 ,
415 bi[3][3] = -90000.0/4081.0 , 421 bi[3][3] = -90000.0/4081.0 ,
416 bi[3][4] = 36000.0/4081.0 , 422 bi[3][4] = 36000.0/4081.0 ,
417 // -------------------------------------- 423 // --------------------------------------------------------
418 // 424 //
419 // COEFFICIENTS OF bi[4] 425 // COEFFICIENTS OF bi[4]
420 bi[4][0] = 0.0 , 426 bi[4][0] = 0.0 ,
421 bi[4][1] = -3125.0/704.0 , 427 bi[4][1] = -3125.0/704.0 ,
422 bi[4][2] = 25625.0/1056.0 , 428 bi[4][2] = 25625.0/1056.0 ,
423 bi[4][3] = -5625.0/176.0 , 429 bi[4][3] = -5625.0/176.0 ,
424 bi[4][4] = 1125.0/88.0 , 430 bi[4][4] = 1125.0/88.0 ,
425 // -------------------------------------- 431 // --------------------------------------------------------
426 // 432 //
427 // COEFFICIENTS OF bi[5] 433 // COEFFICIENTS OF bi[5]
428 bi[5][0] = 0.0 , 434 bi[5][0] = 0.0 ,
429 bi[5][1] = 164025.0/74624.0 , 435 bi[5][1] = 164025.0/74624.0 ,
430 bi[5][2] = -448335.0/37312.0 , 436 bi[5][2] = -448335.0/37312.0 ,
431 bi[5][3] = 295245.0/18656.0 , 437 bi[5][3] = 295245.0/18656.0 ,
432 bi[5][4] = -59049.0/9328.0 , 438 bi[5][4] = -59049.0/9328.0 ,
433 // -------------------------------------- 439 // --------------------------------------------------------
434 // 440 //
435 // COEFFICIENTS OF bi[6] 441 // COEFFICIENTS OF bi[6]
436 bi[6][0] = 0.0 , 442 bi[6][0] = 0.0 ,
437 bi[6][1] = -25.0/28.0 , 443 bi[6][1] = -25.0/28.0 ,
438 bi[6][2] = 205.0/42.0 , 444 bi[6][2] = 205.0/42.0 ,
439 bi[6][3] = -45.0/7.0 , 445 bi[6][3] = -45.0/7.0 ,
440 bi[6][4] = 18.0/7.0 , 446 bi[6][4] = 18.0/7.0 ,
441 // -------------------------------------- 447 // --------------------------------------------------------
442 // 448 //
443 // COEFFICIENTS OF bi[7] 449 // COEFFICIENTS OF bi[7]
444 bi[7][0] = 0.0 , 450 bi[7][0] = 0.0 ,
445 bi[7][1] = -2.0/11.0 , 451 bi[7][1] = -2.0/11.0 ,
446 bi[7][2] = 73.0/55.0 , 452 bi[7][2] = 73.0/55.0 ,
447 bi[7][3] = -171.0/55.0 , 453 bi[7][3] = -171.0/55.0 ,
448 bi[7][4] = 108.0/55.0 , 454 bi[7][4] = 108.0/55.0 ,
449 // -------------------------------------- 455 // --------------------------------------------------------
450 // 456 //
451 // COEFFICIENTS OF bi[8] 457 // COEFFICIENTS OF bi[8]
452 bi[8][0] = 0.0 , 458 bi[8][0] = 0.0 ,
453 bi[8][1] = 189.0/22.0 , 459 bi[8][1] = 189.0/22.0 ,
454 bi[8][2] = -1593.0/55.0 , 460 bi[8][2] = -1593.0/55.0 ,
455 bi[8][3] = 3537.0/110.0 , 461 bi[8][3] = 3537.0/110.0 ,
456 bi[8][4] = -648.0/55.0 , 462 bi[8][4] = -648.0/55.0 ,
457 // -------------------------------------- 463 // --------------------------------------------------------
458 // 464 //
459 // COEFFICIENTS OF bi[9] 465 // COEFFICIENTS OF bi[9]
460 bi[9][0] = 0.0 , 466 bi[9][0] = 0.0 ,
461 bi[9][1] = 351.0/110.0 , 467 bi[9][1] = 351.0/110.0 ,
462 bi[9][2] = -999.0/55.0 , 468 bi[9][2] = -999.0/55.0 ,
463 bi[9][3] = 2943.0/110.0 , 469 bi[9][3] = 2943.0/110.0 ,
464 bi[9][4] = -648.0/55.0 ; 470 bi[9][4] = -648.0/55.0 ;
465 // -------------------------------------- 471 // --------------------------------------------------------
>> 472
466 473
467 for(G4int i = 0; i< numberOfVariables; ++i << 474
468 { << 475 for(G4int i = 0; i< numberOfVariables; i++)
469 yIn[i] = yInput[i]; 476 yIn[i] = yInput[i];
470 } <<
471 477
472 G4double tau0 = tau; 478 G4double tau0 = tau;
473 << 479 // Calculating the polynomials :
474 // Calculating the polynomials << 480
475 // << 481 for(int i=1; i<=9; i++){ //Here i is NOT the coordinate no. , it's stage no.
476 for(auto i=1; i<=9; ++i) // i is NOT the <<
477 { <<
478 b[i] = 0; 482 b[i] = 0;
479 tau = 1.0; 483 tau = 1.0;
480 for(auto j=0; j<=4; ++j) << 484 for(int j=0; j<=4; j++){
481 { <<
482 b[i] += bi[i][j]*tau; 485 b[i] += bi[i][j]*tau;
483 tau*=tau0; 486 tau*=tau0;
484 } 487 }
485 } 488 }
486 489
487 for(G4int i=0; i<numberOfVariables; ++i) / << 490 for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
488 { <<
489 yOut[i] = yIn[i] + Step*tau0*(b[1]*dyd 491 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
490 b[4]*ak4 492 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
491 b[7]*ak7 493 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
492 } 494 }
>> 495
493 } 496 }
>> 497
>> 498
>> 499
>> 500
>> 501
>> 502
>> 503
494 504