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 // G4BogackiShampine45 implementation << 26 // Bogacki-Shampine's RK 5(4) non-FSAL implementation by Somnath Banerjee >> 27 // Supervision / code review: John Apostolakis 27 // 28 // 28 // Bogacki-Shampine's RK 5(4) non-FSAL interpo << 29 // Sponsored by Google in Google Summer of Code 2015. 29 // Definition of the stepper() method that eva << 30 // >> 31 // First version: 25 May 2015 >> 32 // >> 33 // History >> 34 // ----------------------------- >> 35 // Created by Somnath Banerjee on May-August 2015 >> 36 // Improvements by John Apostolakis, May 2016 >> 37 /////////////////////////////////////////////////////////////////////////////// >> 38 // >> 39 // This is the source file of G4BogackiShampine45 class containing the >> 40 // definition of the stepper() method that evaluates one step in 30 // field propagation. 41 // field propagation. 31 // 42 // 32 // The Butcher table of the Bogacki-Shampine-8 43 // The Butcher table of the Bogacki-Shampine-8-4-5 method is: 33 // 44 // 34 // 0 | 45 // 0 | 35 // 1/6 | 1/6 46 // 1/6 | 1/6 36 // 2/9 | 2/27 4/27 << 47 // 2/9 | 2/27 4/27 37 // 3/7 | 183/1372 -162/343 1053/1 48 // 3/7 | 183/1372 -162/343 1053/1372 38 // 2/3 | 68/297 -4/11 42/1 49 // 2/3 | 68/297 -4/11 42/143 1960/3861 39 // 3/4 | 597/22528 81/352 63099/5 50 // 3/4 | 597/22528 81/352 63099/585728 58653/366080 4617/20480 40 // 1 | 174197/959244 -30942/79937 8152137/1 51 // 1 | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219 41 // 1 | 587/8064 0 4440339/1 52 // 1 | 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 42 //-------------------------------------------- 53 //------------------------------------------------------------------------------------------------------------------- 43 // 587/8064 0 4440339/1 54 // 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 0 44 // 2479/34992 0 123/4 55 // 2479/34992 0 123/416 612941/3411720 43/1440 2272/6561 79937/1113912 3293/556956 45 // 56 // 46 // Coefficients have been obtained from: << 57 // Do NOT re-indent the lines above - their meaning becomes lost 47 // http://www.netlib.org/ode/rksuite/ << 58 // ******************************** 48 // << 59 // Coefficients have been obtained from rksuite.f : http://www.netlib.org/ode/rksuite/ 49 // Note on meaning of label "non-FSAL version" << 60 50 // This method calculates the deriviative dy << 61 // Note on meaning of label "non-FSAL version": 51 // integration interval at each step, as par << 62 // This method calculates the deriviative dy/dx at the endpoint of the integration interval at each step. 52 // endpoint and its error. So this value is << 63 // as part of its evaluation of the endpoint and its error. 53 // for re-use in case of a successful step. << 64 // So this value is available to be returned, for re-use in case of a successful step. 54 // (This is done in a 'later' version using << 65 // ( This is done in a 'later' version using a refined interface. ) 55 // << 56 // Created: Somnath Banerjee, Google Summer of << 57 // Revision: John Apostolakis, CERN, May 2016 << 58 // ------------------------------------------- << 59 66 60 #include <cassert> 67 #include <cassert> 61 68 62 #include "G4BogackiShampine45.hh" 69 #include "G4BogackiShampine45.hh" 63 #include "G4LineSection.hh" 70 #include "G4LineSection.hh" 64 71 65 G4bool G4BogackiShampine45::fPreparedConstants << 72 G4bool G4BogackiShampine45::fPreparedConstants= false; 66 G4double G4BogackiShampine45::bi[12][7]; 73 G4double G4BogackiShampine45::bi[12][7]; 67 74 68 // Constructor << 75 //Constructor 69 // << 70 G4BogackiShampine45::G4BogackiShampine45(G4Equ 76 G4BogackiShampine45::G4BogackiShampine45(G4EquationOfMotion *EqRhs, 71 G4int 77 G4int noIntegrationVariables, 72 G4boo 78 G4bool primary) 73 : G4MagIntegratorStepper(EqRhs, noIntegrati << 79 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables), >> 80 fLastStepLength(-1.0), >> 81 fAuxStepper(nullptr), >> 82 fPreparedInterpolation(false) 74 { 83 { >> 84 75 const G4int numberOfVariables = noIntegrat 85 const G4int numberOfVariables = noIntegrationVariables; 76 86 77 // New Chunk of memory being created for u << 87 //New Chunk of memory being created for use by the stepper 78 88 79 // aki - for storing intermediate RHS << 89 //aki - for storing intermediate RHS 80 ak2 = new G4double[numberOfVariables]; 90 ak2 = new G4double[numberOfVariables]; 81 ak3 = new G4double[numberOfVariables]; 91 ak3 = new G4double[numberOfVariables]; 82 ak4 = new G4double[numberOfVariables]; 92 ak4 = new G4double[numberOfVariables]; 83 ak5 = new G4double[numberOfVariables]; 93 ak5 = new G4double[numberOfVariables]; 84 ak6 = new G4double[numberOfVariables]; 94 ak6 = new G4double[numberOfVariables]; 85 ak7 = new G4double[numberOfVariables]; 95 ak7 = new G4double[numberOfVariables]; 86 ak8 = new G4double[numberOfVariables]; 96 ak8 = new G4double[numberOfVariables]; 87 ak9 = new G4double[numberOfVariables]; 97 ak9 = new G4double[numberOfVariables]; 88 ak10 = new G4double[numberOfVariables]; 98 ak10 = new G4double[numberOfVariables]; 89 ak11 = new G4double[numberOfVariables]; 99 ak11 = new G4double[numberOfVariables]; 90 100 91 for (auto & i : p) << 101 for (int i = 0; i < 6; i++) { 92 { << 102 p[i]= new G4double[numberOfVariables]; 93 i= new G4double[numberOfVariables]; << 94 } 103 } 95 104 96 assert ( GetNumberOfStateVariables() >= 8 105 assert ( GetNumberOfStateVariables() >= 8 ); 97 const G4int numStateVars = std::max(noInte 106 const G4int numStateVars = std::max(noIntegrationVariables, 98 GetNum 107 GetNumberOfStateVariables() ); 99 108 100 // Must ensure space extra 'state' variabl 109 // Must ensure space extra 'state' variables exists - i.e. yIn[7] 101 yTemp = new G4double[numStateVars]; 110 yTemp = new G4double[numStateVars]; 102 yIn = new G4double[numStateVars] ; 111 yIn = new G4double[numStateVars] ; 103 112 104 fLastInitialVector = new G4double[numState 113 fLastInitialVector = new G4double[numStateVars] ; 105 fLastFinalVector = new G4double[numStateVa 114 fLastFinalVector = new G4double[numStateVars] ; 106 fLastDyDx = new G4double[numberOfVariables 115 fLastDyDx = new G4double[numberOfVariables]; // Only derivatives 107 116 108 fMidVector = new G4double[numberOfVariable << 117 fMidVector = new G4double[numberOfVariables]; // new G4double[numStateVars]; 109 fMidError = new G4double[numberOfVariable << 118 fMidError = new G4double[numberOfVariables]; // new G4double[numStateVars]; 110 << 119 111 if( ! fPreparedConstants ) 120 if( ! fPreparedConstants ) 112 { << 113 PrepareConstants(); 121 PrepareConstants(); 114 } << 115 122 116 if( primary ) 123 if( primary ) 117 { 124 { 118 fAuxStepper = new G4BogackiShampine45(E 125 fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false); 119 } 126 } 120 } 127 } 121 128 122 // Destructor << 129 123 // << 130 //Destructor 124 G4BogackiShampine45::~G4BogackiShampine45() << 131 G4BogackiShampine45::~G4BogackiShampine45(){ 125 { << 132 //clear all previously allocated memory for stepper and DistChord 126 // Clear all previously allocated memory f << 133 delete[] ak2; 127 // << 134 delete[] ak3; 128 delete [] ak2; << 135 delete[] ak4; 129 delete [] ak3; << 136 delete[] ak5; 130 delete [] ak4; << 137 delete[] ak6; 131 delete [] ak5; << 138 delete[] ak7; 132 delete [] ak6; << 139 delete[] ak8; 133 delete [] ak7; << 140 delete[] ak9; 134 delete [] ak8; << 141 delete[] ak10; 135 delete [] ak9; << 142 delete[] ak11; 136 delete [] ak10; << 143 137 delete [] ak11; << 144 for (int i = 0; i < 6; i++) { 138 << 145 delete[] p[i]; 139 for (auto & i : p) << 146 } 140 { << 147 141 delete [] i; << 148 delete[] yTemp; 142 } << 149 delete[] yIn; 143 << 150 144 delete [] yTemp; << 151 delete[] fLastInitialVector; 145 delete [] yIn; << 152 delete[] fLastFinalVector; 146 << 153 delete[] fLastDyDx; 147 delete [] fLastInitialVector; << 154 delete[] fMidVector; 148 delete [] fLastFinalVector; << 155 delete[] fMidError; 149 delete [] fLastDyDx; << 150 delete [] fMidVector; << 151 delete [] fMidError; << 152 156 153 delete fAuxStepper; 157 delete fAuxStepper; 154 } 158 } 155 159 156 void G4BogackiShampine45::GetLastDydx( G4doubl << 160 // G4double* G4BogackiShampine45::getLastDydx(){ >> 161 // return ak8; >> 162 // } >> 163 >> 164 void >> 165 G4BogackiShampine45::GetLastDydx( G4double dyDxLast[] ) 157 { 166 { 158 const G4int numberOfVariables = GetNumberOf << 167 const G4int numberOfVariables= this->GetNumberOfVariables(); 159 168 160 for(G4int i=0; i < numberOfVariables; ++i ) << 169 for(G4int i=0; i < numberOfVariables; i++ ){ 161 { << 162 dyDxLast[i] = ak9[i]; 170 dyDxLast[i] = ak9[i]; 163 } 171 } 164 } 172 } 165 173 166 // Stepper << 174 //Stepper : 167 // << 175 168 // Passing in the value of yInput[],the first 176 // Passing in the value of yInput[],the first time dydx[] and Step length 169 // Giving back yOut and yErr arrays for output 177 // Giving back yOut and yErr arrays for output and error respectively 170 // << 178 171 void G4BogackiShampine45::Stepper( const G4dou 179 void G4BogackiShampine45::Stepper( const G4double yInput[], 172 const G4dou 180 const G4double DyDx[], 173 G4dou 181 G4double Step, 174 G4dou 182 G4double yOut[], 175 G4dou 183 G4double yErr[] ) 176 { 184 { 177 G4int i; 185 G4int i; 178 186 179 // Constants from the Butcher tableu 187 // Constants from the Butcher tableu 180 // << 181 const G4double 188 const G4double 182 b21 = 1.0/6.0 , 189 b21 = 1.0/6.0 , 183 b31 = 2.0/27.0 , b32 = 4.0/27.0, 190 b31 = 2.0/27.0 , b32 = 4.0/27.0, 184 191 185 b41 = 183.0/1372.0 , b42 = -162.0/343.0 192 b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0, 186 193 187 b51 = 68.0/297.0, b52 = -4.0/11.0, 194 b51 = 68.0/297.0, b52 = -4.0/11.0, 188 b53 = 42.0/143.0, b54 = 1960.0/3861. 195 b53 = 42.0/143.0, b54 = 1960.0/3861.0, 189 196 190 b61 = 597.0/22528.0, b62 = 81.0/352. 197 b61 = 597.0/22528.0, b62 = 81.0/352.0, 191 b63 = 63099.0/585728.0, b64 = 58653.0/3 198 b63 = 63099.0/585728.0, b64 = 58653.0/366080.0, 192 b65 = 4617.0/20480.0, 199 b65 = 4617.0/20480.0, 193 200 194 b71 = 174197.0/959244.0, b72 = -3094 201 b71 = 174197.0/959244.0, b72 = -30942.0/79937.0, 195 b73 = 8152137.0/19744439.0, b74 = 66610 202 b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0, 196 b75 = -29421.0/29068.0, b76 = 48204 203 b75 = -29421.0/29068.0, b76 = 482048.0/414219.0, 197 204 198 b81 = 587.0/8064.0, b82 = 0.0, 205 b81 = 587.0/8064.0, b82 = 0.0, 199 b83 = 4440339.0/15491840.0, b84 = 24353 206 b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0, 200 b85 = 387.0/44800.0, b86 = 2152. 207 b85 = 387.0/44800.0, b86 = 2152.0/5985.0, 201 b87 = 7267.0/94080.0; 208 b87 = 7267.0/94080.0; 202 209 203 // c1 = 2479.0/34992.0, 210 // c1 = 2479.0/34992.0, 204 // c2 = 0.0, 211 // c2 = 0.0, 205 // c3 = 123.0/416.0, 212 // c3 = 123.0/416.0, 206 // c4 = 612941.0/3411720.0, 213 // c4 = 612941.0/3411720.0, 207 // c5 = 43.0/1440.0, 214 // c5 = 43.0/1440.0, 208 // c6 = 2272.0/6561.0, 215 // c6 = 2272.0/6561.0, 209 // c7 = 79937.0/1113912.0, 216 // c7 = 79937.0/1113912.0, 210 // c8 = 3293.0/556956.0, 217 // c8 = 3293.0/556956.0, 211 218 212 // For the embedded higher order method on 219 // For the embedded higher order method only the difference of values 213 // taken and is used directly later (inste 220 // taken and is used directly later (instead of defining the last row 214 // of Butcher table in separate constants 221 // of Butcher table in separate constants and taking the 215 // difference) 222 // difference) 216 // << 217 const G4double 223 const G4double 218 dc1 = b81 - 2479.0 / 34992.0 , 224 dc1 = b81 - 2479.0 / 34992.0 , 219 dc2 = 0.0, 225 dc2 = 0.0, 220 dc3 = b83 - 123.0 / 416.0 , 226 dc3 = b83 - 123.0 / 416.0 , 221 dc4 = b84 - 612941.0 / 3411720.0, 227 dc4 = b84 - 612941.0 / 3411720.0, 222 dc5 = b85 - 43.0 / 1440.0, 228 dc5 = b85 - 43.0 / 1440.0, 223 dc6 = b86 - 2272.0 / 6561.0, 229 dc6 = b86 - 2272.0 / 6561.0, 224 dc7 = b87 - 79937.0 / 1113912.0, 230 dc7 = b87 - 79937.0 / 1113912.0, 225 dc8 = - 3293.0 / 556956.0; 231 dc8 = - 3293.0 / 556956.0; 226 232 227 const G4int numberOfVariables = GetNumberO << 233 const G4int numberOfVariables= this->GetNumberOfVariables(); 228 234 229 // The number of variables to be integrate 235 // The number of variables to be integrated over 230 // << 236 yOut[7] = yTemp[7] = yIn[7]; 231 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; << 237 // Saving yInput because yInput and yOut can be aliases for same array 232 << 238 233 // Saving yInput because yInput and yOut c << 239 for(i=0;i<numberOfVariables;i++) 234 // << 235 for(i=0; i<numberOfVariables; ++i) << 236 { 240 { 237 yIn[i]=yInput[i]; 241 yIn[i]=yInput[i]; 238 } 242 } 239 243 240 // RightHandSide(yIn, dydx) ; 244 // RightHandSide(yIn, dydx) ; 241 // 1st Step - Not doing, getting passed 245 // 1st Step - Not doing, getting passed 242 // << 246 243 for(i=0; i<numberOfVariables; ++i) << 247 for(i=0;i<numberOfVariables;i++) 244 { 248 { 245 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 249 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 246 } 250 } 247 RightHandSide(yTemp, ak2) ; / 251 RightHandSide(yTemp, ak2) ; // 2nd Step 248 252 249 for(i=0; i<numberOfVariables; ++i) << 253 for(i=0;i<numberOfVariables;i++) 250 { 254 { 251 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] 255 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ; 252 } 256 } 253 RightHandSide(yTemp, ak3) ; / 257 RightHandSide(yTemp, ak3) ; // 3rd Step 254 258 255 for(i=0; i<numberOfVariables; ++i) << 259 for(i=0;i<numberOfVariables;i++) 256 { 260 { 257 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] 261 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ; 258 } 262 } 259 RightHandSide(yTemp, ak4) ; / 263 RightHandSide(yTemp, ak4) ; // 4th Step 260 264 261 for(i=0; i<numberOfVariables; ++i) << 265 for(i=0;i<numberOfVariables;i++) 262 { 266 { 263 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] 267 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] + 264 b54*ak4[i]) 268 b54*ak4[i]) ; 265 } 269 } 266 RightHandSide(yTemp, ak5) ; / 270 RightHandSide(yTemp, ak5) ; // 5th Step 267 271 268 for(i=0; i<numberOfVariables; ++i) << 272 for(i=0;i<numberOfVariables;i++) 269 { 273 { 270 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] 274 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] + 271 b64*ak4[i] + 275 b64*ak4[i] + b65*ak5[i]) ; 272 } 276 } 273 RightHandSide(yTemp, ak6) ; / 277 RightHandSide(yTemp, ak6) ; // 6th Step 274 278 275 for(i=0; i<numberOfVariables; ++i) << 279 for(i=0;i<numberOfVariables;i++) 276 { 280 { 277 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] 281 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] + 278 b74*ak4[i] + 282 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]); 279 } 283 } 280 RightHandSide(yTemp, ak7); / << 284 RightHandSide(yTemp, ak7); //7th Step 281 285 282 for(i=0; i<numberOfVariables; ++i) << 286 for(i=0;i<numberOfVariables;i++) 283 { 287 { 284 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + 288 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] + 285 b84*ak4[i] + 289 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] + 286 b87*ak7[i]); 290 b87*ak7[i]); 287 } 291 } 288 RightHandSide(yOut, ak8); / << 292 RightHandSide(yOut, ak8); //8th Step - Final one Using FSAL 289 293 290 for(i=0; i<numberOfVariables; ++i) << 294 for(i=0;i<numberOfVariables;i++) 291 { 295 { 292 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[ 296 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] + 293 dc5*ak5[i] + dc6*ak6[i 297 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ; 294 298 295 // Store Input and Final values, for p 299 // Store Input and Final values, for possible use in calculating chord 296 // << 297 fLastInitialVector[i] = yIn[i] ; 300 fLastInitialVector[i] = yIn[i] ; 298 fLastFinalVector[i] = yOut[i]; 301 fLastFinalVector[i] = yOut[i]; 299 fLastDyDx[i] = DyDx[i]; 302 fLastDyDx[i] = DyDx[i]; 300 } 303 } 301 304 302 fLastStepLength = Step; 305 fLastStepLength = Step; 303 fPreparedInterpolation= false; 306 fPreparedInterpolation= false; 304 307 305 return ; 308 return ; 306 } 309 } 307 310 308 // DistChord << 311 309 // << 312 //The following has not been tested 310 G4double G4BogackiShampine45::DistChord() cons << 313 >> 314 //The DistChord() function fot the class - must define it here. >> 315 G4double G4BogackiShampine45::DistChord() const 311 { 316 { 312 G4double distLine, distChord; 317 G4double distLine, distChord; 313 G4ThreeVector initialPoint, finalPoint, mi 318 G4ThreeVector initialPoint, finalPoint, midPoint; 314 319 315 // Store last initial and final points << 320 // Store last initial and final points (they will be overwritten in self-Stepper call!) 316 // (they will be overwritten in self-Stepp << 321 initialPoint = G4ThreeVector( fLastInitialVector[0], 317 // << 318 initialPoint = G4ThreeVector(fLastInitialV << 319 fLastInitialV 322 fLastInitialVector[1], fLastInitialVector[2]); 320 finalPoint = G4ThreeVector(fLastFinalVec << 323 finalPoint = G4ThreeVector( fLastFinalVector[0], 321 fLastFinalVec 324 fLastFinalVector[1], fLastFinalVector[2]); 322 325 323 #if 1 326 #if 1 324 // Old method -- Do half a step using Step 327 // Old method -- Do half a step using StepNoErr 325 // << 328 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 326 fAuxStepper->Stepper( fLastInitialVector, << 329 fMidVector, fMidError); 327 fMidVector, fMidErro << 328 #else 330 #else 329 // New method -- Using interpolation, << 331 // New method -- Using interpolation, requires only 3 extra stages (ie 3 extra field evaluations ) 330 // requires only 3 extra stages (ie 3 extr << 331 332 332 // Use Interpolation, instead of auxiliary 333 // Use Interpolation, instead of auxiliary stepper to evaluate midpoint 333 // << 334 if( ! fPreparedInterpolation ) { 334 if( ! fPreparedInterpolation ) << 335 G4BogackiShampine45 *cThis= const_cast<G4BogackiShampine45 *>(this); 335 { << 336 cThis-> SetupInterpolationHigh(); // ( fLastInitialVector, fLastDyDx, fLastStepLength ); 336 G4BogackiShampine45* cThis = const_cast << 337 cThis-> SetupInterpolationHigh(); << 338 } 337 } 339 // For calculating the output at the tau f << 338 //For calculating the output at the tau fraction of Step 340 // << 341 G4double tau = 0.5; 339 G4double tau = 0.5; 342 InterpolateHigh( tau, fMidVector ); << 340 // cThis->InterpolateHigh( /* fLastInitialVector, fLastDyDx, fLastStepLength, */, fMidVector, tau); >> 341 // Old arguments: ( /*yInput, dydx, step,*/ yOut, tau ); >> 342 this->InterpolateHigh( tau, fMidVector ); 343 #endif 343 #endif 344 344 345 midPoint = G4ThreeVector( fMidVector[0], f 345 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]); 346 346 347 // Use stored values of Initial and Endpoi 347 // Use stored values of Initial and Endpoint + new Midpoint to evaluate 348 // distance of Chord << 348 // distance of Chord 349 349 350 if (initialPoint != finalPoint) 350 if (initialPoint != finalPoint) 351 { 351 { 352 distLine = G4LineSection::Distline( m << 352 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint ); 353 distChord = distLine; 353 distChord = distLine; 354 } 354 } 355 else 355 else 356 { 356 { 357 distChord = (midPoint-initialPoint).ma 357 distChord = (midPoint-initialPoint).mag(); 358 } 358 } 359 return distChord; 359 return distChord; 360 } 360 } 361 361 362 void G4BogackiShampine45::SetupInterpolationHi 362 void G4BogackiShampine45::SetupInterpolationHigh() >> 363 // ( const G4double *yInput, const G4double *dydx, const G4double Step) 363 { 364 { 364 // Coefficients for the additional stages << 365 //Coefficients for the additional stages : 365 // << 366 const G4double 366 const G4double 367 a91 = 455.0/6144.0 , 367 a91 = 455.0/6144.0 , 368 a92 = 0.0 , 368 a92 = 0.0 , 369 a93 = 10256301.0/35409920.0 , 369 a93 = 10256301.0/35409920.0 , 370 a94 = 2307361.0/17971200.0 , 370 a94 = 2307361.0/17971200.0 , 371 a95 = -387.0/102400.0 , 371 a95 = -387.0/102400.0 , 372 a96 = 73.0/5130.0 , 372 a96 = 73.0/5130.0 , 373 a97 = -7267.0/215040.0 , 373 a97 = -7267.0/215040.0 , 374 a98 = 1.0/32.0 , 374 a98 = 1.0/32.0 , 375 375 376 a101 = -837888343715.0/13176988637184.0 , 376 a101 = -837888343715.0/13176988637184.0 , 377 a102 = 30409415.0/52955362.0 , 377 a102 = 30409415.0/52955362.0 , 378 a103 = -48321525963.0/759168069632.0 , 378 a103 = -48321525963.0/759168069632.0 , 379 a104 = 8530738453321.0/197654829557760.0 , 379 a104 = 8530738453321.0/197654829557760.0 , 380 a105 = 1361640523001.0/1626788720640.0 , 380 a105 = 1361640523001.0/1626788720640.0 , 381 a106 = -13143060689.0/38604458898.0 , 381 a106 = -13143060689.0/38604458898.0 , 382 a107 = 18700221969.0/379584034816.0 , 382 a107 = 18700221969.0/379584034816.0 , 383 a108 = -5831595.0/847285792.0 , 383 a108 = -5831595.0/847285792.0 , 384 a109 = -5183640.0/26477681.0 , 384 a109 = -5183640.0/26477681.0 , 385 385 386 a111 = 98719073263.0/1551965184000.0 , 386 a111 = 98719073263.0/1551965184000.0 , 387 a112 = 1307.0/123552.0 , 387 a112 = 1307.0/123552.0 , 388 a113 = 4632066559387.0/70181753241600.0 , 388 a113 = 4632066559387.0/70181753241600.0 , 389 a114 = 7828594302389.0/382182512025600.0 , 389 a114 = 7828594302389.0/382182512025600.0 , 390 a115 = 40763687.0/11070259200.0 , 390 a115 = 40763687.0/11070259200.0 , 391 a116 = 34872732407.0/224610586200.0 , 391 a116 = 34872732407.0/224610586200.0 , 392 a117 = -2561897.0/30105600.0 , 392 a117 = -2561897.0/30105600.0 , 393 a118 = 1.0/10.0 , 393 a118 = 1.0/10.0 , 394 a119 = -1.0/10.0 , 394 a119 = -1.0/10.0 , 395 a1110 = -1403317093.0/11371610250.0 ; 395 a1110 = -1403317093.0/11371610250.0 ; 396 396 397 const G4int numberOfVariables= this->GetNu 397 const G4int numberOfVariables= this->GetNumberOfVariables(); 398 const G4double* dydx= fLastDyDx; << 398 >> 399 // const G4double *yIn= fLastInitialVector; >> 400 const G4double *dydx= fLastDyDx; 399 const G4double Step = fLastStepLength; 401 const G4double Step = fLastStepLength; 400 402 401 yTemp[7] = yIn[7]; 403 yTemp[7] = yIn[7]; 402 404 403 // Evaluate the extra stages << 405 //Evaluate the extra stages : 404 // << 406 for(int i=0; i<numberOfVariables; i++){ 405 for(G4int i=0; i<numberOfVariables; ++i) << 406 { << 407 yTemp[i] = yIn[i] + Step*(a91*dydx[i] 407 yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] + 408 a94*ak4[i] + 408 a94*ak4[i] + a95*ak5[i] + a96*ak6[i] + 409 a97*ak7[i] + 409 a97*ak7[i] + a98*ak8[i] ); 410 } 410 } 411 411 412 RightHandSide(yTemp, ak9); // 9th stage << 412 RightHandSide(yTemp, ak9); //9th stage 413 413 414 for(G4int i=0; i<numberOfVariables; ++i) << 414 for(int i=0; i<numberOfVariables; i++){ 415 { << 416 yTemp[i] = yIn[i] + Step*(a101*dydx[i] 415 yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] + 417 a104*ak4[i] 416 a104*ak4[i] + a105*ak5[i] + a106*ak6[i] + 418 a107*ak7[i] 417 a107*ak7[i] + a108*ak8[i] + a109*ak9[i] ); 419 } 418 } 420 419 421 RightHandSide(yTemp, ak10); // 10th stage << 420 RightHandSide(yTemp, ak10); //10th stage 422 421 423 for(G4int i=0; i<numberOfVariables; ++i) << 422 for(int i=0; i<numberOfVariables; i++){ 424 { << 425 yTemp[i] = yIn[i] + Step*(a111*dydx[i] 423 yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] + 426 a114*ak4[i] 424 a114*ak4[i] + a115*ak5[i] + a116*ak6[i] + 427 a117*ak7[i] 425 a117*ak7[i] + a118*ak8[i] + a119*ak9[i] + 428 a1110*ak10[i 426 a1110*ak10[i] ); 429 } 427 } 430 RightHandSide(yTemp, ak11); // 11th stage << 428 RightHandSide(yTemp, ak11); //11th stage 431 429 432 // In future we can restrict the number o 430 // In future we can restrict the number of variables interpolated 433 // << 431 int nwant = numberOfVariables; 434 G4int nwant = numberOfVariables; << 435 432 436 // Form the coefficients of the interpola << 433 // Form the coefficients of the interpolating polynomial in its shifted 437 // and scaled form. The terms are groupe << 434 // and scaled form. The terms are grouped to minimize the errors 438 // of the transformation, to cope with il << 435 // of the transformation, to cope with ill-conditioning. ( From RKSUITE ) 439 // << 436 // 440 for (G4int l = 0; l < nwant; ++l) << 437 for (int l = 0; l < nwant; l++) { 441 { << 442 // Coefficient of tau^6 438 // Coefficient of tau^6 443 p[5][l] = bi[5][6]*ak5[l] + 439 p[5][l] = bi[5][6]*ak5[l] + 444 ((bi[10][6]*ak10[l] + bi[8][ 440 ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) + 445 (bi[7][6]*ak7[l] + bi[6][6] 441 (bi[7][6]*ak7[l] + bi[6][6]*ak6[l])) + 446 ((bi[4][6]*ak4[l] + bi[9][6] 442 ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) + 447 (bi[3][6]*ak3[l] + bi[11][6 443 (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) + 448 bi[1][6]*dydx[l]); 444 bi[1][6]*dydx[l]); 449 // Coefficient of tau^5 445 // Coefficient of tau^5 450 p[4][l] = (bi[10][5]*ak10[l] + bi[9][5 446 p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l]) + 451 ((bi[7][5]*ak7[l] + bi[6][5]* 447 ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) + 452 bi[5][5]*ak5[l]) + ((bi[4 448 bi[5][5]*ak5[l]) + ((bi[4][5]*ak4[l] + 453 bi[8][5]*ak8[l]) + (bi[3][5 449 bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] + 454 bi[11][5]*ak11[l]) + bi[1][ 450 bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]); 455 // Coefficient of tau^4 451 // Coefficient of tau^4 456 p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4] 452 p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) + 457 (bi[7][4]*ak7[l] + bi[6][4] 453 (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) + 458 bi[5][4]*ak5[l]) + ((bi[10 454 bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] + 459 bi[9][4]*ak9[l]) + (bi[3] 455 bi[9][4]*ak9[l]) + (bi[3][4]*ak3[l] + 460 bi[11][4]*ak11[l]) + bi[1] 456 bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]); 461 // Coefficient of tau^3 457 // Coefficient of tau^3 462 p[2][l] = bi[5][3]*ak5[l] + bi[6][3]* 458 p[2][l] = bi[5][3]*ak5[l] + bi[6][3]*ak6[l] + 463 ((bi[3][3]*ak3[l] + bi[9][3]* 459 ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) + 464 (bi[10][3]*ak10[l]+ bi[8][3]* 460 (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) + 465 ((bi[4][3]*ak4[l] + bi[11][3] 461 ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]); 466 // Coefficient of tau^2 462 // Coefficient of tau^2 467 p[1][l] = bi[5][2]*ak5[l] + ((bi[6][2 463 p[1][l] = bi[5][2]*ak5[l] + ((bi[6][2]*ak6[l] + 468 bi[8][2]*ak8[l]) + bi[1][2 464 bi[8][2]*ak8[l]) + bi[1][2]*dydx[l]) + 469 ((bi[3][2]*ak3[l] + bi[9][2 465 ((bi[3][2]*ak3[l] + bi[9][2]*ak9[l]) + 470 bi[10][2]*ak10[l])+ ((bi[4][2 466 bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] + 471 bi[11][2]*ak2[l]) + bi[7][2 467 bi[11][2]*ak2[l]) + bi[7][2]*ak7[l]); 472 } 468 } 473 << 469 // 474 // Scale all the coefficients by the st << 470 // Scale all the coefficients by the step size. 475 // << 471 // 476 for (auto & i : p) << 472 for (int i = 0; i < 6; i++) { 477 { << 473 for (int l = 0; l < nwant; l++) { 478 for (G4int l = 0; l < nwant; ++l) << 474 p[i][l] *= Step; 479 { << 480 i[l] *= Step; << 481 } 475 } 482 } 476 } 483 477 484 fPreparedInterpolation = true; << 478 fPreparedInterpolation= true; 485 } 479 } 486 480 487 void G4BogackiShampine45::PrepareConstants() 481 void G4BogackiShampine45::PrepareConstants() 488 { 482 { 489 for(auto i=1; i<= 11; ++i) << 483 for(int i=1; i<= 11; i++) 490 { << 491 bi[i][1] = 0.0 ; 484 bi[i][1] = 0.0 ; 492 } << 493 485 494 for(auto i=1; i<=6; ++i) << 486 for(int i=1; i<=6; i++) 495 { << 496 bi[2][i] = 0.0 ; 487 bi[2][i] = 0.0 ; 497 } << 498 488 499 bi[1][6] = -12134338393.0 / 1050809760.0 , 489 bi[1][6] = -12134338393.0 / 1050809760.0 , 500 bi[1][5] = -1620741229.0 / 50038560.0 , 490 bi[1][5] = -1620741229.0 / 50038560.0 , 501 bi[1][4] = -2048058893.0 / 59875200.0 , 491 bi[1][4] = -2048058893.0 / 59875200.0 , 502 bi[1][3] = -87098480009.0 / 5254048800.0 , 492 bi[1][3] = -87098480009.0 / 5254048800.0 , 503 bi[1][2] = -11513270273.0 / 3502699200.0 , 493 bi[1][2] = -11513270273.0 / 3502699200.0 , 504 // 494 // 505 bi[3][6] = -33197340367.0 / 1218433216.0 495 bi[3][6] = -33197340367.0 / 1218433216.0 , 506 bi[3][5] = -539868024987.0 / 6092166080.0 496 bi[3][5] = -539868024987.0 / 6092166080.0 , 507 bi[3][4] = -39991188681.0 / 374902528.0 , 497 bi[3][4] = -39991188681.0 / 374902528.0 , 508 bi[3][3] = -69509738227.0 / 1218433216.0 498 bi[3][3] = -69509738227.0 / 1218433216.0 , 509 bi[3][2] = -29327744613.0 / 2436866432.0 499 bi[3][2] = -29327744613.0 / 2436866432.0 , 510 // 500 // 511 bi[4][6] = -284800997201.0 / 1990533916 501 bi[4][6] = -284800997201.0 / 19905339168.0 , 512 bi[4][5] = -7896875450471.0 / 16587782640 502 bi[4][5] = -7896875450471.0 / 165877826400.0 , 513 bi[4][4] = -333945812879.0 / 567103680 503 bi[4][4] = -333945812879.0 / 5671036800.0 , 514 bi[4][3] = -16209923456237.0 / 49763347920 504 bi[4][3] = -16209923456237.0 / 497633479200.0 , 515 bi[4][2] = -2382590741699.0 / 33175565280 505 bi[4][2] = -2382590741699.0 / 331755652800.0 , 516 // 506 // 517 bi[5][6] = -540919.0 / 741312.0 , 507 bi[5][6] = -540919.0 / 741312.0 , 518 bi[5][5] = -103626067.0 / 43243200.0 , 508 bi[5][5] = -103626067.0 / 43243200.0 , 519 bi[5][4] = -633779.0 / 211200.0 , 509 bi[5][4] = -633779.0 / 211200.0 , 520 bi[5][3] = -32406787.0 / 18532800.0 , 510 bi[5][3] = -32406787.0 / 18532800.0 , 521 bi[5][2] = -36591193.0 / 86486400.0 , 511 bi[5][2] = -36591193.0 / 86486400.0 , 522 // 512 // 523 bi[6][6] = 7157998304.0 / 374350977.0 , 513 bi[6][6] = 7157998304.0 / 374350977.0 , 524 bi[6][5] = 30405842464.0 / 623918295.0 , 514 bi[6][5] = 30405842464.0 / 623918295.0 , 525 bi[6][4] = 183022264.0 / 5332635.0 , 515 bi[6][4] = 183022264.0 / 5332635.0 , 526 bi[6][3] = -3357024032.0 / 1871754885.0 , 516 bi[6][3] = -3357024032.0 / 1871754885.0 , 527 bi[6][2] = -611586736.0 / 89131185.0 , 517 bi[6][2] = -611586736.0 / 89131185.0 , 528 // 518 // 529 bi[7][6] = -138073.0 / 9408.0 , 519 bi[7][6] = -138073.0 / 9408.0 , 530 bi[7][5] = -719433.0 / 15680.0 , 520 bi[7][5] = -719433.0 / 15680.0 , 531 bi[7][4] = -1620541.0 / 31360.0 , 521 bi[7][4] = -1620541.0 / 31360.0 , 532 bi[7][3] = -385151.0 / 15680.0 , 522 bi[7][3] = -385151.0 / 15680.0 , 533 bi[7][2] = -65403.0 / 15680.0 , 523 bi[7][2] = -65403.0 / 15680.0 , 534 // 524 // 535 bi[8][6] = 1245.0 / 64.0 , 525 bi[8][6] = 1245.0 / 64.0 , 536 bi[8][5] = 3991.0 / 64.0 , 526 bi[8][5] = 3991.0 / 64.0 , 537 bi[8][4] = 4715.0 / 64.0 , 527 bi[8][4] = 4715.0 / 64.0 , 538 bi[8][3] = 2501.0 / 64.0 , 528 bi[8][3] = 2501.0 / 64.0 , 539 bi[8][2] = 149.0 / 16.0 , 529 bi[8][2] = 149.0 / 16.0 , 540 bi[8][1] = 1.0 , 530 bi[8][1] = 1.0 , 541 // 531 // 542 bi[9][6] = 55.0 / 3.0 , 532 bi[9][6] = 55.0 / 3.0 , 543 bi[9][5] = 71.0 , 533 bi[9][5] = 71.0 , 544 bi[9][4] = 103.0 , 534 bi[9][4] = 103.0 , 545 bi[9][3] = 199.0 / 3.0 , 535 bi[9][3] = 199.0 / 3.0 , 546 bi[9][2] = 16.0 , 536 bi[9][2] = 16.0 , 547 // 537 // 548 bi[10][6] = -1774004627.0 / 75810735.0 538 bi[10][6] = -1774004627.0 / 75810735.0 , 549 bi[10][5] = -1774004627.0 / 25270245.0 539 bi[10][5] = -1774004627.0 / 25270245.0 , 550 bi[10][4] = -26477681.0 / 359975.0 540 bi[10][4] = -26477681.0 / 359975.0 , 551 bi[10][3] = -11411880511.0 / 379053675.0 541 bi[10][3] = -11411880511.0 / 379053675.0 , 552 bi[10][2] = -423642896.0 / 126351225.0 542 bi[10][2] = -423642896.0 / 126351225.0 , 553 // 543 // 554 bi[11][6] = 35.0 , 544 bi[11][6] = 35.0 , 555 bi[11][5] = 105.0 , 545 bi[11][5] = 105.0 , 556 bi[11][4] = 117.0 , 546 bi[11][4] = 117.0 , 557 bi[11][3] = 59.0 , 547 bi[11][3] = 59.0 , 558 bi[11][2] = 12.0 ; 548 bi[11][2] = 12.0 ; 559 << 549 fPreparedConstants= true; 560 fPreparedConstants = true; << 561 } 550 } 562 551 563 void G4BogackiShampine45::InterpolateHigh(G4do << 552 void G4BogackiShampine45::InterpolateHigh(G4double tau, G4double *yOut) const >> 553 // ( const G4double *yInput, const G4double *dydx, const G4double Step, G4double *yOut, G4double tau) 564 { 554 { 565 G4int numberOfVariables = GetNumberOfVaria << 555 G4int numberOfVariables = this->GetNumberOfVariables(); >> 556 assert( fPreparedConstants); 566 557 567 G4Exception("G4BogackiShampine45::Interpol 558 G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001", 568 FatalException, "Method is not << 559 FatalException, "Method is not yet validated."); 569 560 570 // const G4double *yIn= fLastInitialVecto 561 // const G4double *yIn= fLastInitialVector; 571 // const G4double *dydx= fLastDyDx; 562 // const G4double *dydx= fLastDyDx; 572 const G4double Step = fLastStepLength; 563 const G4double Step = fLastStepLength; 573 564 >> 565 // for(G4int i = 0; i< numberOfVariables; i++) yIn[i] = yInput[i]; 574 #if 1 566 #if 1 575 G4int nwant = numberOfVariables; 567 G4int nwant = numberOfVariables; 576 const G4int norder= 6; 568 const G4int norder= 6; 577 G4int l, k; 569 G4int l, k; 578 570 579 for (l = 0; l < nwant; ++l) << 571 for (l = 0; l < nwant; l++) { 580 { << 581 yOut[l] = p[norder-1][l] * tau; 572 yOut[l] = p[norder-1][l] * tau; 582 } 573 } 583 for (k = norder - 2; k >= 1; --k) << 574 for (k = norder - 2; k >= 1; k--) { 584 { << 575 for (l = 0; l < nwant; l++) { 585 for (l = 0; l < nwant; ++l) << 586 { << 587 yOut[l] = ( yOut[l] + p[k][l] ) * tau 576 yOut[l] = ( yOut[l] + p[k][l] ) * tau; 588 } 577 } 589 } 578 } 590 for (l = 0; l < nwant; ++l) << 579 for (l = 0; l < nwant; l++) { 591 { << 592 yOut[l] = ( yOut[l] + Step * ak8[l] ) * 580 yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l]; 593 } 581 } 594 // The derivative at the end-point is next 582 // The derivative at the end-point is nextDydx[i] = ak8[i]; 595 #else 583 #else 596 // The scheme tries to do the same as the << 584 // The scheme tries to do the same as the DormandPrince745 routine, but fails 597 // but fails << 598 << 599 G4double b[12]; 585 G4double b[12]; 600 const G4double* dydx = fLastDyDx; << 586 const G4double *dydx= fLastDyDx; 601 587 602 G4double tau0 = tau; 588 G4double tau0 = tau; 603 589 604 for(G4int iStage=1; iStage<=11; ++iStage) << 590 for(int iStage=1; iStage<=11; iStage++){ // iStage = stage number 605 { << 606 b[iStage] = 0.0; 591 b[iStage] = 0.0; 607 tau = tau0; 592 tau = tau0; 608 for(G4int j=6; j>=1; --j) << 593 for(int j=6; j>=1; j--){ // j reversed 609 { << 610 b[iStage] += bi[iStage][j] * tau; 594 b[iStage] += bi[iStage][j] * tau; 611 tau *= tau0; 595 tau *= tau0; 612 } 596 } 613 } 597 } 614 598 615 for(G4int i=0; i<numberOfVariables; ++i) << 599 for(int i=0; i<numberOfVariables; i++){ 616 { << 600 yOut[i] = yIn[i] + Step*(b[1] * dydx[i] + b[2] * ak2[i] + b[3] * ak3[i] + 617 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] << 601 b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] + 618 b[4]*ak4[i] + << 602 b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i] + 619 b[7]*ak7[i] + << 603 b[10] * ak10[i] + b[11] * ak11[i] ); 620 b[10]*ak10[i] << 621 } 604 } 622 #endif 605 #endif 623 } 606 } 624 607