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