Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4RKG3_Stepper implementation << 27 // 23 // 28 // Created: J.Apostolakis, V.Grichine - 30.01. << 24 // $Id: G4RKG3_Stepper.cc,v 1.9 2003/10/31 14:35:55 gcosmo Exp $ >> 25 // GEANT4 tag $Name: geant4-07-01 $ >> 26 // 29 // ------------------------------------------- 27 // ------------------------------------------------------------------- 30 28 31 #include "G4RKG3_Stepper.hh" 29 #include "G4RKG3_Stepper.hh" 32 #include "G4LineSection.hh" 30 #include "G4LineSection.hh" 33 #include "G4Mag_EqRhs.hh" 31 #include "G4Mag_EqRhs.hh" 34 32 35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* Eq << 33 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs *EqRhs) 36 : G4MagIntegratorStepper(EqRhs,6) 34 : G4MagIntegratorStepper(EqRhs,6) 37 { 35 { >> 36 G4Exception("G4RKG3_Stepper::G4RKG3_Stepper()", "NotImplemented", >> 37 FatalException, "Stepper not yet available."); 38 } 38 } 39 39 40 G4RKG3_Stepper::~G4RKG3_Stepper() = default; << 40 G4RKG3_Stepper::~G4RKG3_Stepper() >> 41 { >> 42 } 41 43 42 void G4RKG3_Stepper::Stepper( const G4double y << 44 void G4RKG3_Stepper::Stepper( const G4double yInput[7], 43 const G4double d << 45 const G4double dydx[7], 44 G4double S << 46 G4double Step, 45 G4double y << 47 G4double yOut[7], 46 G4double y << 48 G4double yErr[]) 47 { 49 { 48 G4double B[3]; 50 G4double B[3]; >> 51 // G4double yderiv[6]; >> 52 // G4double alpha2, beta2; 49 G4int nvar = 6 ; 53 G4int nvar = 6 ; 50 G4double by15 = 1. / 15. ; // was 0.06666 << 54 // G4double beTemp2, beta2=0; 51 55 52 G4double yTemp[8], dydxTemp[6], yIn[8]; << 56 G4int i; >> 57 G4double by15 = 1. / 15. ; // was 0.066666666 ; >> 58 G4double yTemp[7], dydxTemp[6], yIn[7] ; >> 59 // Saving yInput because yInput and yOut can be aliases for same array >> 60 for(i=0;i<nvar;i++) yIn[i]=yInput[i]; 53 61 54 // Saving yInput because yInput and yOut ca << 55 // << 56 for(G4int i=0; i<nvar; ++i) << 57 { << 58 yIn[i]=yInput[i]; << 59 } << 60 yIn[6] = yInput[6]; << 61 yIn[7] = yInput[7]; << 62 G4double h = Step * 0.5; 62 G4double h = Step * 0.5; 63 hStep = Step; << 64 // Do two half steps << 65 63 66 StepNoErr(yIn, dydx,h, yTemp,B) ; << 64 // Do two half steps 67 << 65 68 // Store Bfld for DistChord Calculation << 66 69 // << 67 // To obtain B1 ... 70 for(auto i=0; i<3; ++i) << 68 // GetEquationOfMotion()->GetFieldValue(yIn,B); 71 { << 69 // G4RKG3_Stepper::StepWithEst(yIn, dydx, Step, yOut,alpha2, beta2, B1, B2 ); 72 BfldIn[i] = B[i]; << 73 } << 74 // RightHandSide(yTemp,dydxTemp) ; << 75 70 >> 71 StepNoErr(yIn, dydx,h, yTemp,B) ; >> 72 // RightHandSide(yTemp,dydxTemp) ; 76 GetEquationOfMotion()->EvaluateRhsGivenB(yT 73 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ; 77 StepNoErr(yTemp,dydxTemp,h,yOut,B); << 74 StepNoErr(yTemp,dydxTemp,h,yOut,B); // ,beTemp2) ; 78 << 75 // beta2 += beTemp2; >> 76 // beta2 *= 0.5; >> 77 79 // Store midpoint, chord calculation 78 // Store midpoint, chord calculation 80 79 81 fyMidPoint = G4ThreeVector(yTemp[0], yTemp << 80 fyMidPoint = G4ThreeVector( yTemp[0], yTemp[1], yTemp[2]); 82 81 83 // Do a full Step 82 // Do a full Step 84 // << 83 85 h *= 2 ; 84 h *= 2 ; 86 StepNoErr(yIn,dydx,h,yTemp,B); << 85 StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ; 87 for(G4int i=0; i<nvar; ++i) << 86 for(i=0;i<nvar;i++) 88 { 87 { 89 yErr[i] = yOut[i] - yTemp[i] ; 88 yErr[i] = yOut[i] - yTemp[i] ; 90 yOut[i] += yErr[i]*by15 ; // Pr 89 yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy 91 } 90 } 92 91 93 // Store values for DistChord method << 92 // for(i=0;i<ncomp;i++) 94 // << 93 // { 95 fyInitial = G4ThreeVector( yIn[0], yIn[1] << 94 // fyInitial[i] = yIn[i]; 96 fpInitial = G4ThreeVector( yIn[3], yIn[4] << 95 // fyFinal[i] = yOut[i]; >> 96 // } >> 97 >> 98 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]); 97 fyFinal = G4ThreeVector( yOut[0], yOut[1 99 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]); >> 100 // beta2 += beTemp2 ; >> 101 // beta2 *= 0.5 ; >> 102 // NormaliseTangentVector( yOut ); // Deleted 98 } 103 } 99 104 100 // ------------------------------------------- 105 // --------------------------------------------------------------------------- 101 106 102 // Integrator for RK from G3 with evaluation o 107 // Integrator for RK from G3 with evaluation of error in solution and delta 103 // geometry based on naive similarity with the 108 // geometry based on naive similarity with the case of uniform magnetic field. 104 // B1[3] is input and is the first magnetic f 109 // B1[3] is input and is the first magnetic field values 105 // B2[3] is output and is the final magnetic f 110 // B2[3] is output and is the final magnetic field values. 106 // << 111 107 void G4RKG3_Stepper::StepWithEst( const G4doub 112 void G4RKG3_Stepper::StepWithEst( const G4double*, 108 const G4doub 113 const G4double*, 109 G4doub 114 G4double, 110 G4doub 115 G4double*, 111 G4doub 116 G4double&, 112 G4doub 117 G4double&, 113 const G4doub 118 const G4double*, 114 G4doub 119 G4double* ) 115 120 116 { 121 { 117 G4Exception("G4RKG3_Stepper::StepWithEst()", << 122 G4Exception("G4RKG3_Stepper::StepWithEst()", "ObsoleteMethod", 118 FatalException, "Method no longe 123 FatalException, "Method no longer used."); 119 } 124 } 120 125 121 // ------------------------------------------- 126 // ----------------------------------------------------------------- 122 127 123 // Integrator RK Stepper from G3 with only two 128 // Integrator RK Stepper from G3 with only two field evaluation per Step. 124 // It is used in propagation initial Step by s 129 // It is used in propagation initial Step by small substeps after solution 125 // error and delta geometry considerations. B[ 130 // error and delta geometry considerations. B[3] is magnetic field which 126 // is passed from substep to substep. 131 // is passed from substep to substep. 127 // << 132 128 void G4RKG3_Stepper::StepNoErr(const G4double << 133 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7], 129 const G4double << 134 const G4double dydx[7], 130 G4double 135 G4double Step, 131 G4double << 136 G4double tOut[7], 132 G4double << 137 G4double B[3] ) // const 133 138 134 { << 139 { 135 << 140 // Copy and edit the routine above, to delete alpha2, beta2, ... 136 // Copy and edit the routine above, to dele << 141 G4double K1[7],K2[7],K3[7],K4[7] ; 137 // << 142 G4double tTemp[7], yderiv[6] ; 138 G4double K1[7], K2[7], K3[7], K4[7]; << 143 G4int i ; 139 G4double tTemp[8]={0.0}, yderiv[6]={0.0}; << 144 140 << 145 #ifdef END_CODE_G3STEPPER 141 // Need Momentum value to give correct valu << 146 G4Exception(" G4RKG3_Stepper::StepNoErr(): method to be no longer used."); 142 // equation. Integration on unit velocity, << 147 #else 143 << 148 // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ; 144 G4double mom, inverse_mom; << 149 145 const G4double c1=0.5, c2=0.125, c3=1./6.; << 150 for(i=0;i<3;i++) 146 << 147 // Correction for momentum not a velocity << 148 // Need the protection !!! must be not zero << 149 // << 150 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4] << 151 inverse_mom = 1./mom; << 152 for(auto i=0; i<3; ++i) << 153 { 151 { 154 K1[i] = Step * dydx[i+3]*inverse_mom; << 152 K1[i] = Step * dydx[i+3]; 155 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*in << 153 tTemp[i] = tIn[i] + Step*(0.5*tIn[i+3] + 0.125*K1[i]) ; 156 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ; << 154 tTemp[i+3] = tIn[i+3] + 0.5*K1[i] ; 157 } 155 } 158 << 159 GetEquationOfMotion()->EvaluateRhsReturnB(t 156 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; 160 << 157 // Calculates yderiv & returns B too! 161 for(auto i=0; i<3; ++i) << 158 >> 159 for(i=0;i<3;i++) 162 { 160 { 163 K2[i] = Step * yderiv[i+3]*inverse_mom; << 161 K2[i] = Step * yderiv[i+3]; 164 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ; << 162 tTemp[i+3] = tIn[i+3] + 0.5*K2[i] ; 165 } 163 } 166 << 164 167 // Given B, calculate yderiv ! << 165 // Given B, calculate yderiv ! 168 // << 169 GetEquationOfMotion()->EvaluateRhsGivenB(tT 166 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ; 170 << 167 171 for(auto i=0; i<3; ++i) << 168 for(i=0;i<3;i++) 172 { 169 { 173 K3[i] = Step * yderiv[i+3]*inverse_mom; << 170 K3[i] = Step * yderiv[i+3]; 174 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inver << 171 tTemp[i] = tIn[i] + Step*(tIn[i+3] + 0.5*K3[i]) ; 175 tTemp[i+3] = tIn[i+3] + K3[i]*mom ; << 172 tTemp[i+3] = tIn[i+3] + K3[i] ; 176 } 173 } 177 174 178 // Calculates y-deriv(atives) & returns B t << 175 // Calculates y-deriv(atives) & returns B too! 179 // << 180 GetEquationOfMotion()->EvaluateRhsReturnB(t 176 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; 181 177 182 for(auto i=0; i<3; ++i) // Output tr << 178 for(i=0;i<3;i++) // Output trajectory vector 183 { 179 { 184 K4[i] = Step * yderiv[i+3]*inverse_mom; << 180 K4[i] = Step * yderiv[i+3]; 185 tOut[i] = tIn[i] + Step*(tIn[i+3]*invers << 181 tOut[i] = tIn[i] + Step*(tIn[i+3] + (K1[i] + K2[i] + K3[i])/6.0) ; 186 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2 << 182 tOut[i+3] = tIn[i+3] + (K1[i] + 2*K2[i] + 2*K3[i] +K4[i])/6.0 ; 187 } 183 } 188 tOut[6] = tIn[6]; << 184 // NormaliseTangentVector( tOut ); 189 tOut[7] = tIn[7]; << 185 #endif 190 } 186 } 191 187 192 // ------------------------------------------- 188 // --------------------------------------------------------------------------- 193 << 194 G4double G4RKG3_Stepper::DistChord() const << 195 { << 196 // Soon: must check whether h/R > 2 pi !! << 197 // Method below is good only for < 2 pi << 198 189 199 G4double distChord,distLine; << 190 G4double G4RKG3_Stepper::DistChord() const 200 << 191 { 201 if (fyInitial != fyFinal) << 192 // Soon: must check whether h/R > 2 pi !! 202 { << 193 // Method below is good only for < 2 pi 203 distLine = G4LineSection::Distline(fyMid << 204 distChord = distLine; << 205 } << 206 else << 207 { << 208 distChord = (fyMidPoint-fyInitial).mag() << 209 } << 210 194 211 return distChord; << 195 return G4LineSection::Distline( fyMidPoint, fyInitial, fyFinal ); >> 196 // This is a class method that gives distance of Mid >> 197 // from the Chord between the Initial and Final points. 212 } 198 } 213 199