00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id: G4ClassicalRK4.cc 69786 2013-05-15 09:38:51Z gcosmo $ 00028 // 00029 // ------------------------------------------------------------------- 00030 00031 #include "G4ClassicalRK4.hh" 00032 #include "G4ThreeVector.hh" 00033 00035 // 00036 // Constructor sets the number of variables (default = 6) 00037 00038 G4ClassicalRK4:: 00039 G4ClassicalRK4(G4EquationOfMotion* EqRhs, G4int numberOfVariables) 00040 : G4MagErrorStepper(EqRhs, numberOfVariables) 00041 { 00042 unsigned int noVariables= std::max(numberOfVariables,8); // For Time .. 7+1 00043 00044 dydxm = new G4double[noVariables]; 00045 dydxt = new G4double[noVariables]; 00046 yt = new G4double[noVariables]; 00047 } 00048 00050 // 00051 // Destructor 00052 00053 G4ClassicalRK4::~G4ClassicalRK4() 00054 { 00055 delete[] dydxm; 00056 delete[] dydxt; 00057 delete[] yt; 00058 } 00059 00061 // 00062 // Given values for the variables y[0,..,n-1] and their derivatives 00063 // dydx[0,...,n-1] known at x, use the classical 4th Runge-Kutta 00064 // method to advance the solution over an interval h and return the 00065 // incremented variables as yout[0,...,n-1], which not be a distinct 00066 // array from y. The user supplies the routine RightHandSide(x,y,dydx), 00067 // which returns derivatives dydx at x. The source is routine rk4 from 00068 // NRC p. 712-713 . 00069 00070 void 00071 G4ClassicalRK4::DumbStepper( const G4double yIn[], 00072 const G4double dydx[], 00073 G4double h, 00074 G4double yOut[]) 00075 { 00076 const G4int nvar = this->GetNumberOfVariables(); // fNumberOfVariables(); 00077 G4int i; 00078 G4double hh = h*0.5 , h6 = h/6.0 ; 00079 00080 // Initialise time to t0, needed when it is not updated by the integration. 00081 // [ Note: Only for time dependent fields (usually electric) 00082 // is it neccessary to integrate the time.] 00083 yt[7] = yIn[7]; 00084 yOut[7] = yIn[7]; 00085 00086 for(i=0;i<nvar;i++) 00087 { 00088 yt[i] = yIn[i] + hh*dydx[i] ; // 1st Step K1=h*dydx 00089 } 00090 RightHandSide(yt,dydxt) ; // 2nd Step K2=h*dydxt 00091 00092 for(i=0;i<nvar;i++) 00093 { 00094 yt[i] = yIn[i] + hh*dydxt[i] ; 00095 } 00096 RightHandSide(yt,dydxm) ; // 3rd Step K3=h*dydxm 00097 00098 for(i=0;i<nvar;i++) 00099 { 00100 yt[i] = yIn[i] + h*dydxm[i] ; 00101 dydxm[i] += dydxt[i] ; // now dydxm=(K2+K3)/h 00102 } 00103 RightHandSide(yt,dydxt) ; // 4th Step K4=h*dydxt 00104 00105 for(i=0;i<nvar;i++) // Final RK4 output 00106 { 00107 yOut[i] = yIn[i]+h6*(dydx[i]+dydxt[i]+2.0*dydxm[i]); //+K1/6+K4/6+(K2+K3)/3 00108 } 00109 if ( nvar == 12 ) { NormalisePolarizationVector ( yOut ); } 00110 00111 } // end of DumbStepper .................................................... 00112 00114 // 00115 // StepWithEst 00116 00117 void 00118 G4ClassicalRK4::StepWithEst( const G4double*, 00119 const G4double*, 00120 G4double, 00121 G4double*, 00122 G4double&, 00123 G4double&, 00124 const G4double*, 00125 G4double* ) 00126 { 00127 G4Exception("G4ClassicalRK4::StepWithEst()", "GeomField0001", 00128 FatalException, "Method no longer used."); 00129 00130 } // end of StepWithEst ...................................................... 00131