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