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