Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4ConstRK4 implementation 27 // 28 // Created: J.Apostolakis, T.Nikitina - 18.09. 29 // ------------------------------------------- 30 31 #include "G4ConstRK4.hh" 32 #include "G4ThreeVector.hh" 33 #include "G4LineSection.hh" 34 35 ////////////////////////////////////////////// 36 // 37 // Constructor sets the number of *State* vari 38 // The number of variables integrated is alw 39 // 40 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4i 41 : G4MagErrorStepper(EqRhs, 6, numStateVariab 42 { 43 // const G4int numberOfVariables= 6; 44 if( numStateVariables < 8 ) 45 { 46 std::ostringstream message; 47 message << "The number of State variables 48 << "Instead it is - numStateVariab 49 G4Exception("G4ConstRK4::G4ConstRK4()", "G 50 FatalException, message, "Use 51 } 52 53 fEq = EqRhs; 54 yMiddle = new G4double[8]; 55 dydxMid = new G4double[8]; 56 yInitial = new G4double[8]; 57 yOneStep = new G4double[8]; 58 59 dydxm = new G4double[8]; 60 dydxt = new G4double[8]; 61 yt = new G4double[8]; 62 Field[0]=0.; Field[1]=0.; Field[2]=0.; 63 } 64 65 ////////////////////////////////////////////// 66 // 67 // Destructor 68 69 G4ConstRK4::~G4ConstRK4() 70 { 71 delete [] yMiddle; 72 delete [] dydxMid; 73 delete [] yInitial; 74 delete [] yOneStep; 75 delete [] dydxm; 76 delete [] dydxt; 77 delete [] yt; 78 } 79 80 ////////////////////////////////////////////// 81 // 82 // Given values for the variables y[0,..,n-1] 83 // dydx[0,...,n-1] known at x, use the classic 84 // method to advance the solution over an inte 85 // incremented variables as yout[0,...,n-1], w 86 // array from y. The user supplies the routine 87 // which returns derivatives dydx at x. The so 88 // NRC p. 712-713 . 89 // 90 void G4ConstRK4::DumbStepper( const G4double y 91 const G4double d 92 G4double h 93 G4double y 94 { 95 G4double hh = h*0.5 , h6 = h/6.0 ; 96 97 // 1st Step K1=h*dydx 98 yt[5] = yIn[5] + hh*dydx[5] ; 99 yt[4] = yIn[4] + hh*dydx[4] ; 100 yt[3] = yIn[3] + hh*dydx[3] ; 101 yt[2] = yIn[2] + hh*dydx[2] ; 102 yt[1] = yIn[1] + hh*dydx[1] ; 103 yt[0] = yIn[0] + hh*dydx[0] ; 104 RightHandSideConst(yt,dydxt) ; 105 106 // 2nd Step K2=h*dydxt 107 yt[5] = yIn[5] + hh*dydxt[5] ; 108 yt[4] = yIn[4] + hh*dydxt[4] ; 109 yt[3] = yIn[3] + hh*dydxt[3] ; 110 yt[2] = yIn[2] + hh*dydxt[2] ; 111 yt[1] = yIn[1] + hh*dydxt[1] ; 112 yt[0] = yIn[0] + hh*dydxt[0] ; 113 RightHandSideConst(yt,dydxm) ; 114 115 // 3rd Step K3=h*dydxm 116 // now dydxm=(K2+K3)/h 117 yt[5] = yIn[5] + h*dydxm[5] ; 118 dydxm[5] += dydxt[5] ; 119 yt[4] = yIn[4] + h*dydxm[4] ; 120 dydxm[4] += dydxt[4] ; 121 yt[3] = yIn[3] + h*dydxm[3] ; 122 dydxm[3] += dydxt[3] ; 123 yt[2] = yIn[2] + h*dydxm[2] ; 124 dydxm[2] += dydxt[2] ; 125 yt[1] = yIn[1] + h*dydxm[1] ; 126 dydxm[1] += dydxt[1] ; 127 yt[0] = yIn[0] + h*dydxm[0] ; 128 dydxm[0] += dydxt[0] ; 129 RightHandSideConst(yt,dydxt) ; 130 131 // 4th Step K4=h*dydxt 132 yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*d 133 yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*d 134 yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*d 135 yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*d 136 yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*d 137 yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*d 138 139 } // end of DumbStepper ..................... 140 141 ////////////////////////////////////////////// 142 // 143 // Stepper 144 145 void 146 G4ConstRK4::Stepper( const G4double yInput[], 147 const G4double dydx[], 148 G4double hstep, 149 G4double yOutput[], 150 G4double yError [] 151 { 152 const G4int nvar = 6; // number of variabl 153 const G4int maxvar = GetNumberOfStateVariab 154 155 // Correction for Richardson extrapolation 156 G4double correction = 1. / ( (1 << Integra 157 158 G4int i; 159 160 // Saving yInput because yInput and yOutput 161 for (i=0; i<maxvar; ++i) { yInitial[i]= 162 163 // Must copy the part of the state *not* in 164 for (i=nvar; i<maxvar; ++i) { yOutput[i]= 165 166 // yInitial[7]= yInput[7]; // The time is 167 yMiddle[7] = yInput[7]; // Copy the time 168 yOneStep[7] = yInput[7]; // As it contrib 169 // yOutput[7] = yInput[7]; // -> dumb step 170 yError[7] = 0.0; 171 172 G4double halfStep = hstep * 0.5; 173 174 // Do two half steps 175 // 176 GetConstField(yInitial,Field); 177 DumbStepper (yInitial, dydx, halfStep, 178 RightHandSideConst(yMiddle, dydxMid); 179 DumbStepper (yMiddle, dydxMid, halfStep, y 180 181 // Store midpoint, chord calculation 182 // 183 fMidPoint = G4ThreeVector( yMiddle[0], yMi 184 185 // Do a full Step 186 // 187 DumbStepper(yInitial, dydx, hstep, yOneStep 188 for(i=0; i<nvar; ++i) 189 { 190 yError [i] = yOutput[i] - yOneStep[i] ; 191 yOutput[i] += yError[i]*correction ; 192 // Provides accuracy increased by 1 or 193 // Richardson extrapolation 194 } 195 196 fInitialPoint = G4ThreeVector( yInitial[0], 197 fFinalPoint = G4ThreeVector( yOutput[0], 198 199 return; 200 } 201 202 ////////////////////////////////////////////// 203 // 204 // Estimate the maximum distance from the curv 205 // 206 // We estimate this using the distance of the 207 // The method below is good only for angle dev 208 // this restriction should not be a problem fo 209 // which generally cannot integrate accurately 210 // 211 G4double G4ConstRK4::DistChord() const 212 { 213 G4double distLine, distChord; 214 215 if (fInitialPoint != fFinalPoint) 216 { 217 distLine= G4LineSection::Distline( fMidPo 218 // This is a class method that gives di 219 // from the Chord between the Initial a 220 distChord = distLine; 221 } 222 else 223 { 224 distChord = (fMidPoint-fInitialPoint).mag 225 } 226 return distChord; 227 } 228