| 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 // G4DormandPrince745 implementation << 26 // $Id: G4DormandPrince745.cc 107470 2017-11-15 07:14:28Z gcosmo $
27 // 27 //
28 // DormandPrince7 - 5(4) non-FSAL << 28 // Class description:
29 // definition of the stepper() method that eva <<
30 // field propagation. <<
31 // The coefficients and the algorithm have bee <<
32 // 29 //
33 // J. R. Dormand and P. J. Prince, "A family o << 30 // DormandPrince7 - 5(4) non-FSAL
34 // Journal of computational and applied Math., <<
35 // 31 //
36 // The Butcher table of the Dormand-Prince-7-4 << 32 // This is the source file of G4DormandPrince745 class containing the
>> 33 // definition of the stepper() method that evaluates one step in
>> 34 // field propagation.
>> 35 // The coefficients and the algorithm have been adapted from
>> 36 //
>> 37 // Table 2 : Coefficients of RK5(4)7M
>> 38 // ---Ref---
>> 39 // J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,”
>> 40 // Journal of computational and applied …, vol. 6, no. 1, pp. 19–26, 1980.
>> 41 // ------------------
>> 42 //
>> 43 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
37 // 44 //
38 // 0 | 45 // 0 |
39 // 1/5 | 1/5 46 // 1/5 | 1/5
40 // 3/10| 3/40 9/40 << 47 // 3/10| 3/40 9/40
41 // 4/5 | 44/45 56/15 32/9 << 48 // 4/5 | 44/45 −56/15 32/9
42 // 8/9 | 19372/6561 25360/2187 64448/6561 << 49 // 8/9 | 19372/6561 −25360/2187 64448/6561 −212/729
43 // 1 | 9017/3168 355/33 46732/5247 << 50 // 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656
44 // 1 | 35/384 0 500/1113 << 51 // 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84
45 // ---------------------------------------- 52 // ------------------------------------------------------------------------
46 // 35/384 0 500/1113 << 53 // 35/384 0 500/1113 125/192 −2187/6784 11/84 0
47 // 5179/57600 0 7571/16695 << 54 // 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40
>> 55 //
>> 56 //
>> 57 // Implementation by Somnath Banerjee - GSoC 2015
>> 58 // Work supported by Google as part of Google Summer of Code 2015.
>> 59 // Supervision / code review: John Apostolakis
>> 60 //
>> 61 // First version: 25 May 2015 - Somnath Banerjee
>> 62 //
>> 63 // Note: Current version includes 3 versions of 'DistChord' method.
>> 64 // Default is hard-coded interpolation.
48 // 65 //
49 // Created: Somnath Banerjee, Google Summer of <<
50 // Supervision: John Apostolakis, CERN <<
51 // ------------------------------------------- <<
52 <<
53 #include "G4DormandPrince745.hh" 66 #include "G4DormandPrince745.hh"
54 #include "G4LineSection.hh" 67 #include "G4LineSection.hh"
>> 68 #include <cmath>
55 69
56 #include <cstring> << 70 //Constructor
57 << 71 G4DormandPrince745::G4DormandPrince745(G4EquationOfMotion *EqRhs,
58 using namespace field_utils; << 72 G4int noIntegrationVariables,
59 << 73 G4bool primary)
60 // Name of this steppers << 74 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
61 const G4String& G4DormandPrince745::StepperTyp << 75 fAuxStepper(0)
62 { 76 {
63 static G4String _stepperType("G4DormandPrinc << 77 const G4int numberOfVariables = // noIntegrationVariables;
64 return _stepperType; << 78 std::max( noIntegrationVariables,
>> 79 ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
>> 80 // For better alignment with cache-line
>> 81
>> 82 //New Chunk of memory being created for use by the stepper
>> 83
>> 84 //ak_i - for storing intermediate RHS
>> 85 ak2 = new G4double[numberOfVariables];
>> 86 ak3 = new G4double[numberOfVariables];
>> 87 ak4 = new G4double[numberOfVariables];
>> 88 ak5 = new G4double[numberOfVariables];
>> 89 ak6 = new G4double[numberOfVariables];
>> 90 ak7 = new G4double[numberOfVariables];
>> 91 // Also always allocate arrays for interpolation stages
>> 92 ak8 = new G4double[numberOfVariables];
>> 93 ak9 = new G4double[numberOfVariables];
>> 94
>> 95 // Must ensure space for extra 'state' variables exists - i.e. yIn[7]
>> 96 const G4int numStateVars =
>> 97 std::max(noIntegrationVariables,
>> 98 std::max( GetNumberOfStateVariables(), 8)
>> 99 );
>> 100 yTemp = new G4double[numStateVars];
>> 101 yIn = new G4double[numStateVars];
>> 102
>> 103 fLastInitialVector = new G4double[numStateVars] ;
>> 104 fLastFinalVector = new G4double[numStateVars] ;
>> 105
>> 106 // fLastDyDx = new G4double[numberOfVariables];
>> 107
>> 108 fMidVector = new G4double[numStateVars];
>> 109 fMidError = new G4double[numStateVars];
>> 110
>> 111 yTemp = new G4double[numberOfVariables] ;
>> 112 yIn = new G4double[numberOfVariables] ;
>> 113
>> 114 fLastInitialVector = new G4double[numberOfVariables] ;
>> 115 fLastFinalVector = new G4double[numberOfVariables] ;
>> 116 fInitialDyDx = new G4double[numberOfVariables];
>> 117
>> 118 fMidVector = new G4double[numberOfVariables];
>> 119 fMidError = new G4double[numberOfVariables];
>> 120 fAuxStepper = nullptr;
>> 121 if( primary )
>> 122 {
>> 123 fAuxStepper = new G4DormandPrince745(EqRhs, numberOfVariables,
>> 124 !primary);
>> 125 }
>> 126 fLastStepLength = -1.0;
65 } 127 }
66 128
67 // Description of this steppers - plus details << 129 //Destructor
68 const G4String& G4DormandPrince745::StepperDes << 130 G4DormandPrince745::~G4DormandPrince745()
69 { 131 {
70 static G4String _stepperDescription( << 132 //clear all previously allocated memory for stepper and DistChord
71 "Embedeed 5th order Runge-Kutta stepper - << 133 delete[] ak2;
72 return _stepperDescription; << 134 delete[] ak3;
>> 135 delete[] ak4;
>> 136 delete[] ak5;
>> 137 delete[] ak6;
>> 138 delete[] ak7;
>> 139 // Used only for interpolation
>> 140 delete[] ak8;
>> 141 delete[] ak9;
>> 142
>> 143 delete[] yTemp;
>> 144 delete[] yIn;
>> 145
>> 146 delete[] fLastInitialVector;
>> 147 delete[] fLastFinalVector;
>> 148 delete[] fInitialDyDx;
>> 149 delete[] fMidVector;
>> 150 delete[] fMidError;
>> 151
>> 152 delete fAuxStepper;
73 } 153 }
74 154
75 G4DormandPrince745::G4DormandPrince745(G4Equat <<
76 G4int n <<
77 : G4MagIntegratorStepper(equation, noInteg <<
78 { <<
79 } <<
80 155
81 void G4DormandPrince745::Stepper(const G4doubl << 156 // The coefficients and the algorithm have been adapted from
82 const G4doubl << 157 // Table 2 : Coefficients of RK5(4)7M
83 G4doubl << 158 // ---Ref---
84 G4doubl << 159 // J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,”
85 G4doubl << 160 // Journal of computational and applied …, vol. 6, no. 1, pp. 19–26, 1980.
86 G4doubl << 161 // ------------------
87 { <<
88 Stepper(yInput, dydx, hstep, yOutput, yError <<
89 field_utils::copy(dydxOutput, ak7); <<
90 } <<
91 162
92 // Stepper << 163 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
93 // 164 //
>> 165 // 0 |
>> 166 // 1/5 | 1/5
>> 167 // 3/10| 3/40 9/40
>> 168 // 4/5 | 44/45 −56/15 32/9
>> 169 // 8/9 | 19372/6561 −25360/2187 64448/6561 −212/729
>> 170 // 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656
>> 171 // 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84
>> 172 // ------------------------------------------------------------------------
>> 173 // 35/384 0 500/1113 125/192 −2187/6784 11/84 0
>> 174 // 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40
>> 175
>> 176
>> 177 //Stepper :
>> 178
94 // Passing in the value of yInput[],the first 179 // Passing in the value of yInput[],the first time dydx[] and Step length
95 // Giving back yOut and yErr arrays for output 180 // Giving back yOut and yErr arrays for output and error respectively
96 // << 181
97 void G4DormandPrince745::Stepper(const G4doubl 182 void G4DormandPrince745::Stepper(const G4double yInput[],
98 const G4doubl << 183 const G4double DyDx[],
99 G4doubl << 184 G4double Step,
100 G4doubl << 185 G4double yOut[],
101 G4doubl << 186 G4double yErr[] )
102 { 187 {
103 // The various constants defined on the ba << 188 G4int i;
104 // << 189
105 const G4double b21 = 0.2, << 190 //The various constants defined on the basis of butcher tableu
106 b31 = 3.0 / 40.0, b32 = 9.0 << 191 const G4double //G4double - only once
107 b41 = 44.0 / 45.0, b42 = -5 << 192 b21 = 0.2 ,
108 << 193
109 b51 = 19372.0 / 6561.0, b52 = -25360.0 << 194 b31 = 3.0/40.0, b32 = 9.0/40.0 ,
110 b54 = -212.0 / 729.0, << 195
111 << 196 b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
112 b61 = 9017.0 / 3168.0 , b62 = -355.0 << 197
113 b63 = 46732.0 / 5247.0, b64 = 49.0 / 1 << 198 b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
114 b65 = -5103.0 / 18656.0, << 199 b54 = -212.0/729.0 ,
115 << 200
116 b71 = 35.0 / 384.0, b72 = 0., << 201 b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
117 b73 = 500.0 / 1113.0, b74 = 125.0 / 19 << 202 b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
118 b75 = -2187.0 / 6784.0, b76 = 11.0 / 8 << 203 b65 = -5103.0/18656.0 ,
>> 204
>> 205 b71 = 35.0/384.0, b72 = 0.,
>> 206 b73 = 500.0/1113.0, b74 = 125.0/192.0,
>> 207 b75 = -2187.0/6784.0, b76 = 11.0/84.0,
119 208
120 //Sum of columns, sum(bij) = ei 209 //Sum of columns, sum(bij) = ei
121 // e1 = 0. , << 210 // e1 = 0. ,
122 // e2 = 1.0/5.0 , << 211 // e2 = 1.0/5.0 ,
123 // e3 = 3.0/10.0 , << 212 // e3 = 3.0/10.0 ,
124 // e4 = 4.0/5.0 , << 213 // e4 = 4.0/5.0 ,
125 // e5 = 8.0/9.0 , << 214 // e5 = 8.0/9.0 ,
126 // e6 = 1.0 , << 215 // e6 = 1.0 ,
127 // e7 = 1.0 , << 216 // e7 = 1.0 ,
128 217
129 // Difference between the higher and the l << 218 // Difference between the higher and the lower order method coeff. :
130 // b7j are the coefficients of higher orde 219 // b7j are the coefficients of higher order
131 220
132 dc1 = -(b71 - 5179.0 / 57600.0), << 221 dc1 = -( b71 - 5179.0/57600.0),
133 dc2 = -(b72 - .0), << 222 dc2 = -( b72 - .0),
134 dc3 = -(b73 - 7571.0 / 16695.0), << 223 dc3 = -( b73 - 7571.0/16695.0),
135 dc4 = -(b74 - 393.0 / 640.0), << 224 dc4 = -( b74 - 393.0/640.0),
136 dc5 = -(b75 + 92097.0 / 339200.0), << 225 dc5 = -( b75 + 92097.0/339200.0),
137 dc6 = -(b76 - 187.0 / 2100.0), << 226 dc6 = -( b76 - 187.0/2100.0),
138 dc7 = -(- 1.0 / 40.0); << 227 dc7 = -( - 1.0/40.0 ); //end of declaration
139 228
140 const G4int numberOfVariables = GetNumberO << 229 const G4int numberOfVariables= this->GetNumberOfVariables();
141 State yTemp = {0., 0., 0., 0., 0., 0., 0., <<
142 230
143 // The number of variables to be integrate 231 // The number of variables to be integrated over
144 // <<
145 yOut[7] = yTemp[7] = yInput[7]; 232 yOut[7] = yTemp[7] = yInput[7];
146 <<
147 // Saving yInput because yInput and yOut 233 // Saving yInput because yInput and yOut can be aliases for same array
148 // << 234
149 for(G4int i = 0; i < numberOfVariables; ++ << 235 for(i=0;i<numberOfVariables;i++)
150 { 236 {
151 fyIn[i] = yInput[i]; << 237 yIn[i]=yInput[i];
152 } 238 }
153 // RightHandSide(yIn, dydx); // Not don <<
154 239
155 for(G4int i = 0; i < numberOfVariables; ++ << 240 // RightHandSide(yIn, DyDx) ;
>> 241 // 1st Step - Not doing, getting passed
>> 242
>> 243 for(i=0;i<numberOfVariables;i++)
156 { 244 {
157 yTemp[i] = fyIn[i] + b21 * hstep * dyd << 245 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
158 } 246 }
159 RightHandSide(yTemp, ak2); // << 247 RightHandSide(yTemp, ak2) ; // 2nd stage
160 248
161 for(G4int i = 0; i < numberOfVariables; ++ << 249 for(i=0;i<numberOfVariables;i++)
162 { 250 {
163 yTemp[i] = fyIn[i] + hstep * (b31 * dy << 251 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
164 } 252 }
165 RightHandSide(yTemp, ak3); // << 253 RightHandSide(yTemp, ak3) ; // 3rd stage
166 254
167 for(G4int i = 0; i < numberOfVariables; ++ << 255 for(i=0;i<numberOfVariables;i++)
168 { 256 {
169 yTemp[i] = fyIn[i] + hstep * ( << 257 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
170 b41 * dydx[i] + b42 * ak2[i] + b43 <<
171 } 258 }
172 RightHandSide(yTemp, ak4); // << 259 RightHandSide(yTemp, ak4) ; // 4th stage
173 260
174 for(G4int i = 0; i < numberOfVariables; ++ << 261 for(i=0;i<numberOfVariables;i++)
175 { 262 {
176 yTemp[i] = fyIn[i] + hstep * ( << 263 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
177 b51 * dydx[i] + b52 * ak2[i] + b53 << 264 b54*ak4[i]) ;
178 } 265 }
179 RightHandSide(yTemp, ak5); // << 266 RightHandSide(yTemp, ak5) ; // 5th stage
180 267
181 for(G4int i = 0; i < numberOfVariables; ++ << 268 for(i=0;i<numberOfVariables;i++)
182 { 269 {
183 yTemp[i] = fyIn[i] + hstep * ( << 270 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
184 b61 * dydx[i] + b62 * ak2[i] + << 271 b64*ak4[i] + b65*ak5[i]) ;
185 b63 * ak3[i] + b64 * ak4[i] + b65 <<
186 } 272 }
187 RightHandSide(yTemp, ak6); // << 273 RightHandSide(yTemp, ak6) ; // 6th stage
188 274
189 for(G4int i = 0; i < numberOfVariables; ++ << 275 for(i=0;i<numberOfVariables;i++)
190 { 276 {
191 yOut[i] = fyIn[i] + hstep * ( << 277 yOut[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
192 b71 * dydx[i] + b72 * ak2[i] + b73 << 278 b74*ak4[i] + b75*ak5[i] + b76*ak6[i] );
193 b74 * ak4[i] + b75 * ak5[i] + b76 <<
194 } 279 }
195 RightHandSide(yOut, ak7); // << 280 RightHandSide(yOut, ak7); //7th and Final stage
196 281
197 for(G4int i = 0; i < numberOfVariables; ++ << 282 for(i=0;i<numberOfVariables;i++)
198 { 283 {
199 yErr[i] = hstep * ( << 284 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
200 dc1 * dydx[i] + dc2 * ak2[i] + << 285 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) + 1.5e-18 ;
201 dc3 * ak3[i] + dc4 * ak4[i] + <<
202 dc5 * ak5[i] + dc6 * ak6[i] + dc7 <<
203 ) + 1.5e-18; <<
204 286
205 // Store Input and Final values, for p 287 // Store Input and Final values, for possible use in calculating chord
206 // << 288 fLastInitialVector[i] = yIn[i] ;
207 fyOut[i] = yOut[i]; << 289 fLastFinalVector[i] = yOut[i];
208 fdydxIn[i] = dydx[i]; << 290 fInitialDyDx[i] = DyDx[i];
209 } 291 }
210 292
211 fLastStepLength = hstep; << 293 fLastStepLength = Step;
>> 294
>> 295 return ;
212 } 296 }
213 297
214 G4double G4DormandPrince745::DistChord() const <<
215 { <<
216 // Coefficients were taken from Some Pract <<
217 // by Lawrence F. Shampine, page 149, c* <<
218 // <<
219 const G4double hf1 = 6025192743.0 / 300855 <<
220 hf3 = 51252292925.0 / 65400 <<
221 hf4 = - 2691868925.0 / 4512 <<
222 hf5 = 187940372067.0 / 1594 <<
223 hf6 = - 1776094331.0 / 1974 <<
224 hf7 = 11237099.0 / 23504338 <<
225 298
226 G4ThreeVector mid; << 299 // Calculate DistChord given start, mid and end-point of step
>> 300 G4double G4DormandPrince745::DistLine( G4double yStart[], G4double yMid[], G4double yEnd[] ) const
>> 301 {
>> 302 G4double distLine, distChord;
>> 303 G4ThreeVector initialPoint, finalPoint, midPoint;
>> 304
>> 305 initialPoint = G4ThreeVector( yStart[0], yStart[1], yStart[2]);
>> 306 finalPoint = G4ThreeVector( yEnd[0], yEnd[1], yEnd[2]);
>> 307 midPoint = G4ThreeVector( yMid[0], yMid[1], yMid[2]);
227 308
228 for(G4int i = 0; i < 3; ++i) << 309 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
>> 310 // distance of Chord
>> 311 if (initialPoint != finalPoint)
>> 312 {
>> 313 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
>> 314 distChord = distLine;
>> 315 }
>> 316 else
229 { 317 {
230 mid[i] = fyIn[i] + 0.5 * fLastStepLeng << 318 distChord = (midPoint-initialPoint).mag();
231 hf1 * fdydxIn[i] + hf3 * ak3[i] + <<
232 hf4 * ak4[i] + hf5 * ak5[i] + hf6 <<
233 } 319 }
>> 320 return distChord;
>> 321 }
>> 322
>> 323 // (New) DistChord function using interpolation
>> 324 G4double G4DormandPrince745::DistChord2() const
>> 325 {
>> 326 // Copy the values of stages from this (original) into the Aux Stepper
>> 327 *fAuxStepper = *this;
234 328
235 const G4ThreeVector begin = makeVector(fyI << 329 //Preparing for the interpolation
236 const G4ThreeVector end = makeVector(fyOut << 330 fAuxStepper->SetupInterpolation(); // (fLastInitialVector, fInitialDyDx, fLastStepLength);
>> 331 //Interpolate to half step
>> 332 fAuxStepper->Interpolate( /*fLastInitialVector, fInitialDyDx, fLastStepLength,*/ 0.5, fAuxStepper->fMidVector);
237 333
238 return G4LineSection::Distline(mid, begin, << 334 return DistLine( fLastInitialVector, fAuxStepper->fMidVector, fLastFinalVector);
239 } 335 }
240 336
241 // The lower (4th) order interpolant given by << 337 G4double G4DormandPrince745::DistChord() const
242 // J. R. Dormand and P. J. Prince, "Run <<
243 // Computers & Mathematics with Applica <<
244 // pp. 1007-1017, 1986. <<
245 // <<
246 void G4DormandPrince745:: <<
247 Interpolate4thOrder(G4double yOut[], G4double <<
248 { 338 {
249 const G4int numberOfVariables = GetNumberO << 339 //Coefficients for halfway interpolation
250 << 340 const G4double
251 const G4double tau2 = tau * tau, << 341 hf1 = 5783653.0/57600000.0 ,
252 tau3 = tau * tau2, << 342 hf2 = 0. ,
253 tau4 = tau2 * tau2; << 343 hf3 = 466123.0/1192500.0 ,
>> 344 hf4 = -41347.0/1920000.0 ,
>> 345 hf5 = 16122321.0/339200000.0 ,
>> 346 hf6 = -7117.0/20000.0,
>> 347 hf7 = 183.0/10000.0 ;
>> 348
>> 349 for(int i=0; i<3; i++){
>> 350 fMidVector[i] = fLastInitialVector[i] + fLastStepLength*(
>> 351 hf1*fInitialDyDx[i] + hf2*ak2[i] + hf3*ak3[i] + hf4*ak4[i] +
>> 352 hf5*ak5[i] + hf6*ak6[i] + hf7*ak7[i] );
>> 353 }
>> 354
>> 355 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
>> 356 // distance of Chord
254 357
255 const G4double bf1 = 1.0 / 11282082432.0 * << 358 return DistLine( fLastInitialVector, fMidVector, fLastFinalVector);
256 157015080.0 * tau4 - 13107642775.0 * t << 359 }
257 32272833064.0 * tau + 11282082432.0); <<
258 360
259 const G4double bf3 = - 100.0 / 32700410799 << 361 //The original DistChord() function for the class
260 15701508.0 * tau3 - 914128567.0 * tau2 << 362 G4double G4DormandPrince745::DistChord3() const
261 1323431896.0); << 363 {
262 << 364 // Do half a step using StepNoErr
263 const G4double bf4 = 25.0 / 5641041216.0 * << 365 fAuxStepper->Stepper( fLastInitialVector, fInitialDyDx, 0.5 * fLastStepLength,
264 94209048.0 * tau3 - 1518414297.0 * tau << 366 fAuxStepper->fMidVector, fAuxStepper->fMidError) ;
265 889289856.0); << 367 return DistLine( fLastInitialVector, fAuxStepper->fMidVector, fLastFinalVector);
266 << 368 }
267 const G4double bf5 = - 2187.0 / 1993167896 <<
268 52338360.0 * tau3 - 451824525.0 * tau2 <<
269 259006536.0); <<
270 369
271 const G4double bf6 = 11.0 / 2467955532.0 * << 370 // The lower (4th) order interpolant given by Dormand and prince
272 106151040.0 * tau3 - 661884105.0 * tau << 371 // "An RK 5(4) triple"
273 946554244.0 * tau - 361440756.0); << 372 //---Ref---
274 << 373 // J. R. Dormand and P. J. Prince, “Runge-Kutta triples,”
275 const G4double bf7 = 1.0 / 29380423.0 * ta << 374 // Computers & Mathematics with Applications, vol. 12, no. 9,
276 8293050.0 * tau2 - 82437520.0 * tau + << 375 // pp. 1007–1017, 1986.
>> 376 //---------------------------
277 377
278 for(G4int i = 0; i < numberOfVariables; ++ << 378 void G4DormandPrince745::SetupInterpolation_low() // const G4double *yInput, const G4double *dydx, const G4double Step)
279 { << 379 {
280 yOut[i] = fyIn[i] + fLastStepLength * << 380 //Nothing to be done
281 bf1 * fdydxIn[i] + bf3 * ak3[i] + << 381 }
282 bf5 * ak5[i] + bf6 * ak6[i] + bf7 << 382
>> 383 void G4DormandPrince745::Interpolate_low( /* const G4double yInput[],
>> 384 const G4double dydx[],
>> 385 const G4double Step, */
>> 386 G4double yOut[],
>> 387 G4double tau )
>> 388 {
>> 389 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
>> 390 // Coefficients for all the seven stages.
>> 391 G4double Step = fLastStepLength;
>> 392 const G4double *dydx= fInitialDyDx;
>> 393
>> 394 const G4int numberOfVariables= this->GetNumberOfVariables();
>> 395
>> 396 // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
>> 397
>> 398 const G4double
>> 399 tau_2 = tau * tau,
>> 400 tau_3 = tau * tau_2,
>> 401 tau_4 = tau_2 * tau_2;
>> 402
>> 403 bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
>> 404 + 11282082432.0)/11282082432.0,
>> 405 bf2 = 0.0 ,
>> 406 bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2 + 2074956840.0*tau
>> 407 - 1323431896.0)/32700410799.0,
>> 408 bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2+ 2460397220.0*tau - 889289856.0)/5641041216.0 ,
>> 409 bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2 + 687873124.0*tau - 259006536.0)/199316789632.0 ,
>> 410 bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2 + 946554244.0*tau - 361440756.0)/2467955532.0 ,
>> 411 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)/ 29380423.0 ;
>> 412
>> 413 //Putting together the coefficients calculated as the respective stage coefficients
>> 414 for( int i=0; i<numberOfVariables; i++){
>> 415 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
>> 416 + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
283 } 417 }
284 } 418 }
285 419
>> 420
286 // Following interpolant of order 5 was given 421 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
287 // T. S. Baker, J. R. Dormand, J. P. Gi << 422 //---Ref---
288 // "Continuous approximation with embed << 423 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
289 // Applied Numerical Mathematics, vol. << 424 // “Continuous approximation with embedded Runge-Kutta methods,”
290 // << 425 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51–62, 1996.
291 // Calculating the extra stages for the interp << 426 //---------------------
292 // << 427
293 void G4DormandPrince745::SetupInterpolation5th << 428 // Calculating the extra stages for the interpolant :
294 { << 429 void G4DormandPrince745::SetupInterpolation_high( /* const G4double yInput[],
295 // Coefficients for the additional stages << 430 const G4double dydx[],
296 // << 431 const G4double Step */ ){
297 const G4double b81 = 6245.0 / 62208.0, << 432
298 b82 = 0.0, << 433 //Coefficients for the additional stages :
299 b83 = 8875.0 / 103032.0, << 434 const G4double
300 b84 = -125.0 / 1728.0, << 435 b81 = 6245.0/62208.0 ,
301 b85 = 801.0 / 13568.0, << 436 b82 = 0.0 ,
302 b86 = -13519.0 / 368064.0, << 437 b83 = 8875.0/103032.0 ,
303 b87 = 11105.0 / 368064.0, << 438 b84 = -125.0/1728.0 ,
304 << 439 b85 = 801.0/13568.0 ,
305 b91 = 632855.0 / 4478976.0, << 440 b86 = -13519.0/368064.0 ,
306 b92 = 0.0, << 441 b87 = 11105.0/368064.0 ,
307 b93 = 4146875.0 / 6491016.0 << 442
308 b94 = 5490625.0 /14183424.0 << 443 b91 = 632855.0/4478976.0 ,
309 b95 = -15975.0 / 108544.0, << 444 b92 = 0.0 ,
310 b96 = 8295925.0 / 220286304 << 445 b93 = 4146875.0/6491016.0 ,
311 b97 = -1779595.0 / 62938944 << 446 b94 = 5490625.0/14183424.0 ,
312 b98 = -805.0 / 4104.0; << 447 b95 = -15975.0/108544.0 ,
>> 448 b96 = 8295925.0/220286304.0 ,
>> 449 b97 = -1779595.0/62938944.0 ,
>> 450 b98 = -805.0/4104.0 ;
>> 451
>> 452 const G4int numberOfVariables= this->GetNumberOfVariables();
>> 453 const G4double *dydx = fInitialDyDx;
>> 454 const G4double Step = fLastStepLength;
313 455
314 const G4int numberOfVariables = GetNumberO << 456 // Saving yInput because yInput and yOut can be aliases for same array
315 State yTemp = {0., 0., 0., 0., 0., 0., 0., << 457 // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
>> 458 // yTemp[7] = yIn[7];
316 459
317 // Evaluate the extra stages << 460 //Evaluate the extra stages :
318 // << 461 for(int i=0;i<numberOfVariables;i++)
319 for(G4int i = 0; i < numberOfVariables; ++ <<
320 { 462 {
321 yTemp[i] = fyIn[i] + fLastStepLength * << 463 yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
322 b81 * fdydxIn[i] + b82 * ak2[i] + << 464 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
323 b84 * ak4[i] + b85 * ak5[i] + b86 << 465 b87*ak7[i]);
324 b87 * ak7[i] <<
325 ); <<
326 } 466 }
327 RightHandSide(yTemp, ak8); // 8th << 467 RightHandSide(yTemp, ak8); //8th Stage
328 468
329 for(G4int i = 0; i < numberOfVariables; ++ << 469 for(int i=0;i<numberOfVariables;i++)
330 { 470 {
331 yTemp[i] = fyIn[i] + fLastStepLength * << 471 yTemp[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
332 b91 * fdydxIn[i] + b92 * ak2[i] + << 472 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
333 b94 * ak4[i] + b95 * ak5[i] + b96 << 473 b97*ak7[i] + b98*ak8[i] );
334 b97 * ak7[i] + b98 * ak8[i] <<
335 ); <<
336 } 474 }
337 RightHandSide(yTemp, ak9); // 9th << 475 RightHandSide(yTemp, ak9); //9th Stage
338 } 476 }
339 477
>> 478
340 // Calculating the interpolated result yOut wi 479 // Calculating the interpolated result yOut with the coefficients
341 // << 480 void G4DormandPrince745::Interpolate_high( /* const G4double yInput[],
342 void G4DormandPrince745:: << 481 const G4double dydx[],
343 Interpolate5thOrder(G4double yOut[], G4double << 482 const G4double Step, */
344 { << 483 G4double yOut[],
345 // Define the coefficients for the polynom << 484 G4double tau ){
346 // << 485 //Define the coefficients for the polynomials
347 G4double bi[10][5]; << 486 G4double bi[10][5], b[10];
>> 487 const G4int numberOfVariables = this->GetNumberOfVariables();
>> 488 const G4double *dydx = fInitialDyDx;
>> 489 // const G4double fullStep = fLastStepLength;
>> 490
>> 491 // If given requestedStep in argument:
>> 492 // G4double tau = requestedStep / fLastStepLength;
348 493
349 // COEFFICIENTS OF bi[1] 494 // COEFFICIENTS OF bi[1]
350 bi[1][0] = 1.0, << 495 bi[1][0] = 1.0 ,
351 bi[1][1] = -38039.0 / 7040.0, << 496 bi[1][1] = -38039.0/7040.0 ,
352 bi[1][2] = 125923.0 / 10560.0, << 497 bi[1][2] = 125923.0/10560.0 ,
353 bi[1][3] = -19683.0 / 1760.0, << 498 bi[1][3] = -19683.0/1760.0 ,
354 bi[1][4] = 3303.0 / 880.0, << 499 bi[1][4] = 3303.0/880.0 ,
355 // -------------------------------------- 500 // --------------------------------------------------------
356 // 501 //
357 // COEFFICIENTS OF bi[2] 502 // COEFFICIENTS OF bi[2]
358 bi[2][0] = 0.0, << 503 bi[2][0] = 0.0 ,
359 bi[2][1] = 0.0, << 504 bi[2][1] = 0.0 ,
360 bi[2][2] = 0.0, << 505 bi[2][2] = 0.0 ,
361 bi[2][3] = 0.0, << 506 bi[2][3] = 0.0 ,
362 bi[2][4] = 0.0, << 507 bi[2][4] = 0.0 ,
363 // -------------------------------------- 508 // --------------------------------------------------------
364 // 509 //
365 // COEFFICIENTS OF bi[3] 510 // COEFFICIENTS OF bi[3]
366 bi[3][0] = 0.0, << 511 bi[3][0] = 0.0 ,
367 bi[3][1] = -12500.0 / 4081.0, << 512 bi[3][1] = -12500.0/4081.0 ,
368 bi[3][2] = 205000.0 / 12243.0, << 513 bi[3][2] = 205000.0/12243.0 ,
369 bi[3][3] = -90000.0 / 4081.0, << 514 bi[3][3] = -90000.0/4081.0 ,
370 bi[3][4] = 36000.0 / 4081.0, << 515 bi[3][4] = 36000.0/4081.0 ,
371 // -------------------------------------- 516 // --------------------------------------------------------
372 // 517 //
373 // COEFFICIENTS OF bi[4] 518 // COEFFICIENTS OF bi[4]
374 bi[4][0] = 0.0, << 519 bi[4][0] = 0.0 ,
375 bi[4][1] = -3125.0 / 704.0, << 520 bi[4][1] = -3125.0/704.0 ,
376 bi[4][2] = 25625.0 / 1056.0, << 521 bi[4][2] = 25625.0/1056.0 ,
377 bi[4][3] = -5625.0 / 176.0, << 522 bi[4][3] = -5625.0/176.0 ,
378 bi[4][4] = 1125.0 / 88.0, << 523 bi[4][4] = 1125.0/88.0 ,
379 // -------------------------------------- 524 // --------------------------------------------------------
380 // 525 //
381 // COEFFICIENTS OF bi[5] 526 // COEFFICIENTS OF bi[5]
382 bi[5][0] = 0.0, << 527 bi[5][0] = 0.0 ,
383 bi[5][1] = 164025.0 / 74624.0, << 528 bi[5][1] = 164025.0/74624.0 ,
384 bi[5][2] = -448335.0 / 37312.0, << 529 bi[5][2] = -448335.0/37312.0 ,
385 bi[5][3] = 295245.0 / 18656.0, << 530 bi[5][3] = 295245.0/18656.0 ,
386 bi[5][4] = -59049.0 / 9328.0, << 531 bi[5][4] = -59049.0/9328.0 ,
387 // -------------------------------------- 532 // --------------------------------------------------------
388 // 533 //
389 // COEFFICIENTS OF bi[6] 534 // COEFFICIENTS OF bi[6]
390 bi[6][0] = 0.0, << 535 bi[6][0] = 0.0 ,
391 bi[6][1] = -25.0 / 28.0, << 536 bi[6][1] = -25.0/28.0 ,
392 bi[6][2] = 205.0 / 42.0, << 537 bi[6][2] = 205.0/42.0 ,
393 bi[6][3] = -45.0 / 7.0, << 538 bi[6][3] = -45.0/7.0 ,
394 bi[6][4] = 18.0 / 7.0, << 539 bi[6][4] = 18.0/7.0 ,
395 // -------------------------------------- 540 // --------------------------------------------------------
396 // 541 //
397 // COEFFICIENTS OF bi[7] 542 // COEFFICIENTS OF bi[7]
398 bi[7][0] = 0.0, << 543 bi[7][0] = 0.0 ,
399 bi[7][1] = -2.0 / 11.0, << 544 bi[7][1] = -2.0/11.0 ,
400 bi[7][2] = 73.0 / 55.0, << 545 bi[7][2] = 73.0/55.0 ,
401 bi[7][3] = -171.0 / 55.0, << 546 bi[7][3] = -171.0/55.0 ,
402 bi[7][4] = 108.0 / 55.0, << 547 bi[7][4] = 108.0/55.0 ,
403 // -------------------------------------- 548 // --------------------------------------------------------
404 // 549 //
405 // COEFFICIENTS OF bi[8] 550 // COEFFICIENTS OF bi[8]
406 bi[8][0] = 0.0, << 551 bi[8][0] = 0.0 ,
407 bi[8][1] = 189.0 / 22.0, << 552 bi[8][1] = 189.0/22.0 ,
408 bi[8][2] = -1593.0 / 55.0, << 553 bi[8][2] = -1593.0/55.0 ,
409 bi[8][3] = 3537.0 / 110.0, << 554 bi[8][3] = 3537.0/110.0 ,
410 bi[8][4] = -648.0 / 55.0, << 555 bi[8][4] = -648.0/55.0 ,
411 // -------------------------------------- 556 // --------------------------------------------------------
412 // 557 //
413 // COEFFICIENTS OF bi[9] 558 // COEFFICIENTS OF bi[9]
414 bi[9][0] = 0.0, << 559 bi[9][0] = 0.0 ,
415 bi[9][1] = 351.0 / 110.0, << 560 bi[9][1] = 351.0/110.0 ,
416 bi[9][2] = -999.0 / 55.0, << 561 bi[9][2] = -999.0/55.0 ,
417 bi[9][3] = 2943.0 / 110.0, << 562 bi[9][3] = 2943.0/110.0 ,
418 bi[9][4] = -648.0 / 55.0; << 563 bi[9][4] = -648.0/55.0 ;
419 // -------------------------------------- 564 // --------------------------------------------------------
420 << 565
421 // Calculating the polynomials << 566 // for(G4int i = 0; i< numberOfVariables; i++) { yIn[i] = yInput[i]; }
422 << 567
423 G4double b[10]; << 568 // Calculating the polynomials :
424 std::memset(b, 0.0, sizeof(b)); << 569 #if 1
425 << 570 for(int iStage=1; iStage<=9; iStage++){
426 G4double tauPower = 1.0; << 571 b[iStage] = 0;
427 for(G4int j = 0; j <= 4; ++j) << 572 }
428 { << 573
429 for(G4int iStage = 1; iStage <= 9; ++iS << 574 for(int j=0; j<=4; j++){
430 { << 575 G4double tauPower = 1.0;
431 b[iStage] += bi[iStage][j] * tauPowe << 576 for(int iStage=1; iStage<=9; iStage++){
>> 577 b[iStage] += bi[iStage][j]*tauPower;
432 } 578 }
433 tauPower *= tau; 579 tauPower *= tau;
434 } 580 }
>> 581 #else
>> 582 G4double tau0 = tau;
435 583
436 const G4int numberOfVariables = GetNumberO << 584 for(int i=1; i<=9; i++){ //Here i is NOT the coordinate no. , it's stage no.
437 const G4double stepLen = fLastStepLength * << 585 b[i] = 0;
438 for(G4int i = 0; i < numberOfVariables; ++ << 586 tau = 1.0;
439 { << 587 for(int j=0; j<=4; j++){
440 yOut[i] = fyIn[i] + stepLen * ( << 588 b[i] += bi[i][j]*tau;
441 b[1] * fdydxIn[i] + b[2] * ak2[i] << 589 tau*=tau0;
442 b[4] * ak4[i] + b[5] * ak5[i] + b[ << 590 }
443 b[7] * ak7[i] + b[8] * ak8[i] + b[ << 591 }
444 ); << 592 #endif
>> 593
>> 594 G4double stepLen = fLastStepLength * tau;
>> 595 for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
>> 596 yOut[i] = yIn[i] + stepLen *(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
>> 597 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
>> 598 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
445 } 599 }
>> 600
>> 601 }
>> 602
>> 603 //G4DormandPrince745::G4DormandPrince745(G4DormandPrince745& DP_Obj){
>> 604 //
>> 605 //}
>> 606
>> 607 // Overloaded = operator
>> 608 G4DormandPrince745& G4DormandPrince745::operator=(const G4DormandPrince745& right)
>> 609 {
>> 610 // this->G4DormandPrince745(right.GetEquationOfMotion(),right.GetNumberOfVariables(), false);
>> 611
>> 612 int noVars = right.GetNumberOfVariables();
>> 613 for(int i =0; i< noVars; i++)
>> 614 {
>> 615 this->ak2[i] = right.ak2[i];
>> 616 this->ak3[i] = right.ak3[i];
>> 617 this->ak4[i] = right.ak4[i];
>> 618 this->ak5[i] = right.ak5[i];
>> 619 this->ak6[i] = right.ak6[i];
>> 620 this->ak7[i] = right.ak7[i];
>> 621 this->ak8[i] = right.ak8[i];
>> 622 this->ak9[i] = right.ak9[i];
>> 623
>> 624 this->fInitialDyDx[i] = right.fInitialDyDx[i];
>> 625 this->fLastInitialVector[i] = right.fLastInitialVector[i];
>> 626 this->fMidVector[i] = right.fMidVector[i];
>> 627 this->fMidError[i] = right.fMidError[i];
>> 628 }
>> 629
>> 630 this->fLastStepLength = right.fLastStepLength;
>> 631
>> 632 return *this;
446 } 633 }
447 634