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 // G4CashKarpRKF45 implementation << 26 // >> 27 // $Id: G4CashKarpRKF45.cc 97598 2016-06-06 07:19:46Z gcosmo $ 27 // 28 // 28 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 meth 29 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth 29 // order method (giving fifth-order accuracy) 30 // order method (giving fifth-order accuracy) for the solution of an ODE. 30 // Two different fourth order estimates are ca 31 // Two different fourth order estimates are calculated; their difference 31 // gives an error estimate. [ref. Numerical Re 32 // gives an error estimate. [ref. Numerical Recipes in C, 2nd Edition] 32 // It is used to integrate the equations of th 33 // It is used to integrate the equations of the motion of a particle 33 // in a magnetic field. 34 // in a magnetic field. 34 // 35 // 35 // [ref. Numerical Recipes in C, 2nd Edition] << 36 // [ref. Numerical Recipes in C, 2nd Edition] 36 // 37 // 37 // Authors: J.Apostolakis, V.Grichine - 30.01. << 38 // ------------------------------------------- 38 // ------------------------------------------------------------------- 39 39 40 #include "G4CashKarpRKF45.hh" 40 #include "G4CashKarpRKF45.hh" 41 #include "G4LineSection.hh" 41 #include "G4LineSection.hh" 42 42 43 ////////////////////////////////////////////// 43 ///////////////////////////////////////////////////////////////////// 44 // 44 // 45 // Constructor 45 // Constructor 46 // << 46 47 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfM 47 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfMotion *EqRhs, 48 G4int noInteg << 48 G4int noIntegrationVariables, 49 G4bool primar << 49 G4bool primary) 50 : G4MagIntegratorStepper(EqRhs, noIntegratio << 50 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables), >> 51 fLastStepLength(0.), fAuxStepper(0) 51 { 52 { 52 const G4int numberOfVariables = 53 const G4int numberOfVariables = 53 std::max( noIntegrationVariables, 54 std::max( noIntegrationVariables, 54 ( ( (noIntegrationVariables-1)/ 55 ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) ); 55 // For better alignment with cache-line 56 // For better alignment with cache-line 56 57 57 ak2 = new G4double[numberOfVariables] ; 58 ak2 = new G4double[numberOfVariables] ; 58 ak3 = new G4double[numberOfVariables] ; 59 ak3 = new G4double[numberOfVariables] ; 59 ak4 = new G4double[numberOfVariables] ; 60 ak4 = new G4double[numberOfVariables] ; 60 ak5 = new G4double[numberOfVariables] ; 61 ak5 = new G4double[numberOfVariables] ; 61 ak6 = new G4double[numberOfVariables] ; 62 ak6 = new G4double[numberOfVariables] ; 62 // ak7 = 0; 63 // ak7 = 0; 63 64 64 // Must ensure space extra 'state' variables 65 // Must ensure space extra 'state' variables exists - i.e. yIn[7] 65 const G4int numStateMax = std::max(GetNumbe 66 const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8); 66 const G4int numStateVars = std::max(noIntegr 67 const G4int numStateVars = std::max(noIntegrationVariables, 67 numState 68 numStateMax ); 68 // GetNumbe 69 // GetNumberOfStateVariables() ); 69 70 70 yTemp = new G4double[numStateVars] ; 71 yTemp = new G4double[numStateVars] ; 71 yIn = new G4double[numStateVars] ; 72 yIn = new G4double[numStateVars] ; 72 73 73 fLastInitialVector = new G4double[numStateVa 74 fLastInitialVector = new G4double[numStateVars] ; 74 fLastFinalVector = new G4double[numStateVars 75 fLastFinalVector = new G4double[numStateVars] ; 75 fLastDyDx = new G4double[numberOfVariables]; 76 fLastDyDx = new G4double[numberOfVariables]; 76 77 77 fMidVector = new G4double[numStateVars]; 78 fMidVector = new G4double[numStateVars]; 78 fMidError = new G4double[numStateVars]; 79 fMidError = new G4double[numStateVars]; 79 if( primary ) 80 if( primary ) 80 { 81 { 81 fAuxStepper = new G4CashKarpRKF45(EqRhs, n 82 fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary); 82 } 83 } 83 } 84 } 84 85 85 ////////////////////////////////////////////// 86 ///////////////////////////////////////////////////////////////////// 86 // 87 // 87 // Destructor 88 // Destructor 88 // << 89 89 G4CashKarpRKF45::~G4CashKarpRKF45() 90 G4CashKarpRKF45::~G4CashKarpRKF45() 90 { 91 { 91 delete [] ak2; << 92 delete[] ak2; 92 delete [] ak3; << 93 delete[] ak3; 93 delete [] ak4; << 94 delete[] ak4; 94 delete [] ak5; << 95 delete[] ak5; 95 delete [] ak6; << 96 delete[] ak6; 96 // delete [] ak7; << 97 // delete[] ak7; 97 delete [] yTemp; << 98 delete[] yTemp; 98 delete [] yIn; << 99 delete[] yIn; 99 << 100 100 delete [] fLastInitialVector; << 101 delete[] fLastInitialVector; 101 delete [] fLastFinalVector; << 102 delete[] fLastFinalVector; 102 delete [] fLastDyDx; << 103 delete[] fLastDyDx; 103 delete [] fMidVector; << 104 delete[] fMidVector; 104 delete [] fMidError; << 105 delete[] fMidError; 105 106 106 delete fAuxStepper; 107 delete fAuxStepper; 107 } 108 } 108 109 109 ////////////////////////////////////////////// 110 ////////////////////////////////////////////////////////////////////// 110 // 111 // 111 // Given values for n = 6 variables yIn[0,..., 112 // Given values for n = 6 variables yIn[0,...,n-1] 112 // known at x, use the fifth-order Cash-Karp 113 // known at x, use the fifth-order Cash-Karp Runge- 113 // Kutta-Fehlberg-4-5 method to advance the so 114 // Kutta-Fehlberg-4-5 method to advance the solution over an interval 114 // Step and return the incremented variables a 115 // Step and return the incremented variables as yOut[0,...,n-1]. Also 115 // return an estimate of the local truncation 116 // return an estimate of the local truncation error yErr[] using the 116 // embedded 4th-order method. The user supplie 117 // embedded 4th-order method. The user supplies routine 117 // RightHandSide(y,dydx), which returns deriva 118 // RightHandSide(y,dydx), which returns derivatives dydx for y . 118 // << 119 119 void 120 void 120 G4CashKarpRKF45::Stepper(const G4double yInput 121 G4CashKarpRKF45::Stepper(const G4double yInput[], 121 const G4double dydx[] 122 const G4double dydx[], 122 G4double Step, 123 G4double Step, 123 G4double yOut[] 124 G4double yOut[], 124 G4double yErr[] 125 G4double yErr[]) 125 { 126 { 126 // const G4int nvar = 6 ; 127 // const G4int nvar = 6 ; 127 // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 128 // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875; 128 G4int i; << 129 G4int i; 129 << 130 const G4double b21 = 0.2 , << 131 b31 = 3.0/40.0 , b32 = 9.0/4 << 132 b41 = 0.3 , b42 = -0.9 , b43 << 133 << 134 b51 = -11.0/54.0 , b52 = 2.5 << 135 b54 = 35.0/27.0 , << 136 << 137 b61 = 1631.0/55296.0 , b62 = << 138 b63 = 575.0/13824.0 , b64 = << 139 b65 = 253.0/4096.0 , << 140 << 141 c1 = 37.0/378.0 , c3 = 250.0 << 142 c6 = 512.0/1771.0 , dc5 = -2 << 143 130 144 const G4double dc1 = c1 - 2825.0/27648.0 , << 131 const G4double b21 = 0.2 , 145 dc4 = c4 - 13525.0/55296.0 , << 132 b31 = 3.0/40.0 , b32 = 9.0/40.0 , >> 133 b41 = 0.3 , b42 = -0.9 , b43 = 1.2 , >> 134 >> 135 b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 , >> 136 b54 = 35.0/27.0 , >> 137 >> 138 b61 = 1631.0/55296.0 , b62 = 175.0/512.0 , >> 139 b63 = 575.0/13824.0 , b64 = 44275.0/110592.0 , >> 140 b65 = 253.0/4096.0 , >> 141 >> 142 c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 , >> 143 c6 = 512.0/1771.0 , >> 144 dc5 = -277.0/14336.0 ; >> 145 >> 146 const G4double dc1 = c1 - 2825.0/27648.0 , dc3 = c3 - 18575.0/48384.0 , >> 147 dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ; >> 148 >> 149 // Initialise time to t0, needed when it is not updated by the integration. >> 150 // [ Note: Only for time dependent fields (usually electric) >> 151 // is it neccessary to integrate the time.] >> 152 yOut[7] = yTemp[7] = yIn[7]; >> 153 >> 154 const G4int numberOfVariables= this->GetNumberOfVariables(); >> 155 // The number of variables to be integrated over >> 156 >> 157 // Saving yInput because yInput and yOut can be aliases for same array >> 158 >> 159 for(i=0;i<numberOfVariables;i++) >> 160 { >> 161 yIn[i]=yInput[i]; >> 162 } >> 163 // RightHandSide(yIn, dydx) ; // 1st Step >> 164 >> 165 for(i=0;i<numberOfVariables;i++) >> 166 { >> 167 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; >> 168 } >> 169 RightHandSide(yTemp, ak2) ; // 2nd Step 146 170 147 // Initialise time to t0, needed when it is << 171 for(i=0;i<numberOfVariables;i++) 148 // [ Note: Only for time dependent fi << 172 { 149 // is it neccessary to inte << 150 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; << 151 << 152 const G4int numberOfVariables= this->GetNumb << 153 // The number of variables to be integrate << 154 << 155 // Saving yInput because yInput and yOut ca << 156 << 157 for(i=0; i<numberOfVariables; ++i) << 158 { << 159 yIn[i]=yInput[i]; << 160 } << 161 // RightHandSide(yIn, dydx) ; // << 162 << 163 for(i=0; i<numberOfVariables; ++i) << 164 { << 165 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; << 166 } << 167 RightHandSide(yTemp, ak2) ; // << 168 << 169 for(i=0; i<numberOfVariables; ++i) << 170 { << 171 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b3 173 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ; 172 } << 174 } 173 RightHandSide(yTemp, ak3) ; // << 175 RightHandSide(yTemp, ak3) ; // 3rd Step 174 176 175 for(i=0; i<numberOfVariables; ++i) << 177 for(i=0;i<numberOfVariables;i++) 176 { << 178 { 177 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b4 179 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ; 178 } << 180 } 179 RightHandSide(yTemp, ak4) ; // << 181 RightHandSide(yTemp, ak4) ; // 4th Step 180 182 181 for(i=0; i<numberOfVariables; ++i) << 183 for(i=0;i<numberOfVariables;i++) 182 { << 184 { 183 yTemp[i] = yIn[i] + Step*(b51*dydx[i] << 185 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + 184 + b52*ak2[i] + b53*ak3[i << 186 b54*ak4[i]) ; 185 } << 187 } 186 RightHandSide(yTemp, ak5) ; // << 188 RightHandSide(yTemp, ak5) ; // 5th Step >> 189 >> 190 for(i=0;i<numberOfVariables;i++) >> 191 { >> 192 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + >> 193 b64*ak4[i] + b65*ak5[i]) ; >> 194 } >> 195 RightHandSide(yTemp, ak6) ; // 6th Step 187 196 188 for(i=0; i<numberOfVariables; ++i) << 197 for(i=0;i<numberOfVariables;i++) 189 { << 198 { 190 yTemp[i] = yIn[i] + Step*(b61*dydx[i] << 191 + b62*ak2[i] + b63*ak3[i << 192 } << 193 RightHandSide(yTemp, ak6) ; // << 194 << 195 for(i=0; i<numberOfVariables; ++i) << 196 { << 197 // Accumulate increments with proper weigh 199 // Accumulate increments with proper weights 198 // << 200 199 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*a 201 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ; 200 202 201 // Estimate error as difference between 4t << 203 // Estimate error as difference between 4th and 202 // << 204 // 5th order methods 203 yErr[i] = Step*(dc1*dydx[i] << 205 204 + dc3*ak3[i] + dc4*ak4[i] + dc5*ak << 206 yErr[i] = Step*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] + >> 207 dc5*ak5[i] + dc6*ak6[i]) ; 205 208 206 // Store Input and Final values, for possi 209 // Store Input and Final values, for possible use in calculating chord 207 // << 208 fLastInitialVector[i] = yIn[i] ; 210 fLastInitialVector[i] = yIn[i] ; 209 fLastFinalVector[i] = yOut[i]; 211 fLastFinalVector[i] = yOut[i]; 210 fLastDyDx[i] = dydx[i]; 212 fLastDyDx[i] = dydx[i]; 211 } << 213 } 212 // NormaliseTangentVector( yOut ); // Not wa << 214 // NormaliseTangentVector( yOut ); // Not wanted 213 215 214 fLastStepLength = Step; << 216 fLastStepLength =Step; 215 217 216 return; << 218 return ; 217 } 219 } 218 220 219 ////////////////////////////////////////////// 221 /////////////////////////////////////////////////////////////////////////////// 220 // << 222 221 void 223 void 222 G4CashKarpRKF45::StepWithEst( const G4double*, 224 G4CashKarpRKF45::StepWithEst( const G4double*, 223 const G4double*, 225 const G4double*, 224 G4double, 226 G4double, 225 G4double*, 227 G4double*, 226 G4double&, 228 G4double&, 227 G4double&, 229 G4double&, 228 const G4double*, 230 const G4double*, 229 G4double* 231 G4double* ) 230 { 232 { 231 G4Exception("G4CashKarpRKF45::StepWithEst()" 233 G4Exception("G4CashKarpRKF45::StepWithEst()", "GeomField0001", 232 FatalException, "Method no longe 234 FatalException, "Method no longer used."); 233 return ; 235 return ; 234 } 236 } 235 237 236 ////////////////////////////////////////////// 238 ///////////////////////////////////////////////////////////////// 237 // << 239 238 G4double G4CashKarpRKF45::DistChord() const 240 G4double G4CashKarpRKF45::DistChord() const 239 { 241 { 240 G4double distLine, distChord; 242 G4double distLine, distChord; 241 G4ThreeVector initialPoint, finalPoint, midP 243 G4ThreeVector initialPoint, finalPoint, midPoint; 242 244 243 // Store last initial and final points << 245 // Store last initial and final points (they will be overwritten in self-Stepper call!) 244 // (they will be overwritten in self-Stepper << 245 // << 246 initialPoint = G4ThreeVector( fLastInitialVe 246 initialPoint = G4ThreeVector( fLastInitialVector[0], 247 fLastInitialVe 247 fLastInitialVector[1], fLastInitialVector[2]); 248 finalPoint = G4ThreeVector( fLastFinalVect 248 finalPoint = G4ThreeVector( fLastFinalVector[0], 249 fLastFinalVect 249 fLastFinalVector[1], fLastFinalVector[2]); 250 250 251 // Do half a step using StepNoErr 251 // Do half a step using StepNoErr 252 // << 252 253 fAuxStepper->Stepper( fLastInitialVector, fL << 253 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 254 0.5 * fLastStepLength, << 254 fMidVector, fMidError ); 255 255 256 midPoint = G4ThreeVector( fMidVector[0], fMi 256 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]); 257 257 258 // Use stored values of Initial and Endpoint 258 // Use stored values of Initial and Endpoint + new Midpoint to evaluate 259 // distance of Chord << 259 // distance of Chord 260 // << 260 >> 261 261 if (initialPoint != finalPoint) 262 if (initialPoint != finalPoint) 262 { 263 { 263 distLine = G4LineSection::Distline( midP 264 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint ); 264 distChord = distLine; 265 distChord = distLine; 265 } 266 } 266 else 267 else 267 { 268 { 268 distChord = (midPoint-initialPoint).mag() 269 distChord = (midPoint-initialPoint).mag(); 269 } 270 } 270 return distChord; 271 return distChord; 271 } 272 } >> 273 >> 274 272 275