#include <G4CashKarpRKF45.hh>
Inheritance diagram for G4CashKarpRKF45:
Public Member Functions | |
| G4CashKarpRKF45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true) | |
| ~G4CashKarpRKF45 () | |
| void | Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) |
| G4double | DistChord () const |
| G4int | IntegratorOrder () const |
Definition at line 50 of file G4CashKarpRKF45.hh.
| G4CashKarpRKF45::G4CashKarpRKF45 | ( | G4EquationOfMotion * | EqRhs, | |
| G4int | numberOfVariables = 6, |
|||
| G4bool | primary = true | |||
| ) |
Definition at line 47 of file G4CashKarpRKF45.cc.
00050 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables), 00051 fLastStepLength(0.), fAuxStepper(0) 00052 { 00053 const G4int numberOfVariables = noIntegrationVariables; 00054 00055 ak2 = new G4double[numberOfVariables] ; 00056 ak3 = new G4double[numberOfVariables] ; 00057 ak4 = new G4double[numberOfVariables] ; 00058 ak5 = new G4double[numberOfVariables] ; 00059 ak6 = new G4double[numberOfVariables] ; 00060 ak7 = 0; 00061 yTemp = new G4double[numberOfVariables] ; 00062 yIn = new G4double[numberOfVariables] ; 00063 00064 fLastInitialVector = new G4double[numberOfVariables] ; 00065 fLastFinalVector = new G4double[numberOfVariables] ; 00066 fLastDyDx = new G4double[numberOfVariables]; 00067 00068 fMidVector = new G4double[numberOfVariables]; 00069 fMidError = new G4double[numberOfVariables]; 00070 if( primary ) 00071 { 00072 fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary); 00073 } 00074 }
| G4CashKarpRKF45::~G4CashKarpRKF45 | ( | ) |
Definition at line 80 of file G4CashKarpRKF45.cc.
00081 { 00082 delete[] ak2; 00083 delete[] ak3; 00084 delete[] ak4; 00085 delete[] ak5; 00086 delete[] ak6; 00087 // delete[] ak7; 00088 delete[] yTemp; 00089 delete[] yIn; 00090 00091 delete[] fLastInitialVector; 00092 delete[] fLastFinalVector; 00093 delete[] fLastDyDx; 00094 delete[] fMidVector; 00095 delete[] fMidError; 00096 00097 delete fAuxStepper; 00098 }
| G4double G4CashKarpRKF45::DistChord | ( | ) | const [virtual] |
Implements G4MagIntegratorStepper.
Definition at line 230 of file G4CashKarpRKF45.cc.
References G4LineSection::Distline(), and Stepper().
00231 { 00232 G4double distLine, distChord; 00233 G4ThreeVector initialPoint, finalPoint, midPoint; 00234 00235 // Store last initial and final points (they will be overwritten in self-Stepper call!) 00236 initialPoint = G4ThreeVector( fLastInitialVector[0], 00237 fLastInitialVector[1], fLastInitialVector[2]); 00238 finalPoint = G4ThreeVector( fLastFinalVector[0], 00239 fLastFinalVector[1], fLastFinalVector[2]); 00240 00241 // Do half a step using StepNoErr 00242 00243 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 00244 fMidVector, fMidError ); 00245 00246 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]); 00247 00248 // Use stored values of Initial and Endpoint + new Midpoint to evaluate 00249 // distance of Chord 00250 00251 00252 if (initialPoint != finalPoint) 00253 { 00254 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint ); 00255 distChord = distLine; 00256 } 00257 else 00258 { 00259 distChord = (midPoint-initialPoint).mag(); 00260 } 00261 return distChord; 00262 }
| G4int G4CashKarpRKF45::IntegratorOrder | ( | ) | const [inline, virtual] |
| void G4CashKarpRKF45::Stepper | ( | const G4double | y[], | |
| const G4double | dydx[], | |||
| G4double | h, | |||
| G4double | yout[], | |||
| G4double | yerr[] | |||
| ) | [virtual] |
Implements G4MagIntegratorStepper.
Definition at line 111 of file G4CashKarpRKF45.cc.
References G4MagIntegratorStepper::GetNumberOfVariables(), and G4MagIntegratorStepper::RightHandSide().
Referenced by DistChord().
00116 { 00117 // const G4int nvar = 6 ; 00118 // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875; 00119 G4int i; 00120 00121 const G4double b21 = 0.2 , 00122 b31 = 3.0/40.0 , b32 = 9.0/40.0 , 00123 b41 = 0.3 , b42 = -0.9 , b43 = 1.2 , 00124 00125 b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 , 00126 b54 = 35.0/27.0 , 00127 00128 b61 = 1631.0/55296.0 , b62 = 175.0/512.0 , 00129 b63 = 575.0/13824.0 , b64 = 44275.0/110592.0 , 00130 b65 = 253.0/4096.0 , 00131 00132 c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 , 00133 c6 = 512.0/1771.0 , 00134 dc5 = -277.0/14336.0 ; 00135 00136 const G4double dc1 = c1 - 2825.0/27648.0 , dc3 = c3 - 18575.0/48384.0 , 00137 dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ; 00138 00139 // Initialise time to t0, needed when it is not updated by the integration. 00140 // [ Note: Only for time dependent fields (usually electric) 00141 // is it neccessary to integrate the time.] 00142 yOut[7] = yTemp[7] = yIn[7]; 00143 00144 const G4int numberOfVariables= this->GetNumberOfVariables(); 00145 // The number of variables to be integrated over 00146 00147 // Saving yInput because yInput and yOut can be aliases for same array 00148 00149 for(i=0;i<numberOfVariables;i++) 00150 { 00151 yIn[i]=yInput[i]; 00152 } 00153 // RightHandSide(yIn, dydx) ; // 1st Step 00154 00155 for(i=0;i<numberOfVariables;i++) 00156 { 00157 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; 00158 } 00159 RightHandSide(yTemp, ak2) ; // 2nd Step 00160 00161 for(i=0;i<numberOfVariables;i++) 00162 { 00163 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ; 00164 } 00165 RightHandSide(yTemp, ak3) ; // 3rd Step 00166 00167 for(i=0;i<numberOfVariables;i++) 00168 { 00169 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ; 00170 } 00171 RightHandSide(yTemp, ak4) ; // 4th Step 00172 00173 for(i=0;i<numberOfVariables;i++) 00174 { 00175 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + 00176 b54*ak4[i]) ; 00177 } 00178 RightHandSide(yTemp, ak5) ; // 5th Step 00179 00180 for(i=0;i<numberOfVariables;i++) 00181 { 00182 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + 00183 b64*ak4[i] + b65*ak5[i]) ; 00184 } 00185 RightHandSide(yTemp, ak6) ; // 6th Step 00186 00187 for(i=0;i<numberOfVariables;i++) 00188 { 00189 // Accumulate increments with proper weights 00190 00191 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ; 00192 00193 // Estimate error as difference between 4th and 00194 // 5th order methods 00195 00196 yErr[i] = Step*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] + 00197 dc5*ak5[i] + dc6*ak6[i]) ; 00198 00199 // Store Input and Final values, for possible use in calculating chord 00200 fLastInitialVector[i] = yIn[i] ; 00201 fLastFinalVector[i] = yOut[i]; 00202 fLastDyDx[i] = dydx[i]; 00203 } 00204 // NormaliseTangentVector( yOut ); // Not wanted 00205 00206 fLastStepLength =Step; 00207 00208 return ; 00209 }
1.4.7