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 // G4RKG3_Stepper implementation 26 // G4RKG3_Stepper implementation 27 // 27 // 28 // Created: J.Apostolakis, V.Grichine - 30.01. 28 // Created: J.Apostolakis, V.Grichine - 30.01.1997 29 // ------------------------------------------- 29 // ------------------------------------------------------------------- 30 30 31 #include "G4RKG3_Stepper.hh" 31 #include "G4RKG3_Stepper.hh" 32 #include "G4LineSection.hh" 32 #include "G4LineSection.hh" 33 #include "G4Mag_EqRhs.hh" 33 #include "G4Mag_EqRhs.hh" 34 34 35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* Eq 35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* EqRhs) 36 : G4MagIntegratorStepper(EqRhs,6) 36 : G4MagIntegratorStepper(EqRhs,6) 37 { 37 { 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[], // [8] 43 const G4double d 45 const G4double dydx[], // [6] 44 G4double S 46 G4double Step, 45 G4double y 47 G4double yOut[], // [8] 46 G4double y 48 G4double yErr[] ) 47 { 49 { 48 G4double B[3]; 50 G4double B[3]; 49 G4int nvar = 6 ; 51 G4int nvar = 6 ; 50 G4double by15 = 1. / 15. ; // was 0.06666 52 G4double by15 = 1. / 15. ; // was 0.066666666 ; 51 53 52 G4double yTemp[8], dydxTemp[6], yIn[8]; 54 G4double yTemp[8], dydxTemp[6], yIn[8]; 53 55 54 // Saving yInput because yInput and yOut ca 56 // Saving yInput because yInput and yOut can be aliases for same array 55 // 57 // 56 for(G4int i=0; i<nvar; ++i) 58 for(G4int i=0; i<nvar; ++i) 57 { 59 { 58 yIn[i]=yInput[i]; 60 yIn[i]=yInput[i]; 59 } 61 } 60 yIn[6] = yInput[6]; 62 yIn[6] = yInput[6]; 61 yIn[7] = yInput[7]; 63 yIn[7] = yInput[7]; 62 G4double h = Step * 0.5; 64 G4double h = Step * 0.5; 63 hStep = Step; 65 hStep = Step; 64 // Do two half steps 66 // Do two half steps 65 67 66 StepNoErr(yIn, dydx,h, yTemp,B) ; 68 StepNoErr(yIn, dydx,h, yTemp,B) ; 67 69 68 // Store Bfld for DistChord Calculation 70 // Store Bfld for DistChord Calculation 69 // 71 // 70 for(auto i=0; i<3; ++i) 72 for(auto i=0; i<3; ++i) 71 { 73 { 72 BfldIn[i] = B[i]; 74 BfldIn[i] = B[i]; 73 } 75 } 74 // RightHandSide(yTemp,dydxTemp) ; 76 // RightHandSide(yTemp,dydxTemp) ; 75 77 76 GetEquationOfMotion()->EvaluateRhsGivenB(yT 78 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ; 77 StepNoErr(yTemp,dydxTemp,h,yOut,B); 79 StepNoErr(yTemp,dydxTemp,h,yOut,B); 78 80 79 // Store midpoint, chord calculation 81 // Store midpoint, chord calculation 80 82 81 fyMidPoint = G4ThreeVector(yTemp[0], yTemp 83 fyMidPoint = G4ThreeVector(yTemp[0], yTemp[1], yTemp[2]); 82 84 83 // Do a full Step 85 // Do a full Step 84 // 86 // 85 h *= 2 ; 87 h *= 2 ; 86 StepNoErr(yIn,dydx,h,yTemp,B); 88 StepNoErr(yIn,dydx,h,yTemp,B); 87 for(G4int i=0; i<nvar; ++i) 89 for(G4int i=0; i<nvar; ++i) 88 { 90 { 89 yErr[i] = yOut[i] - yTemp[i] ; 91 yErr[i] = yOut[i] - yTemp[i] ; 90 yOut[i] += yErr[i]*by15 ; // Pr 92 yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy 91 } 93 } 92 94 93 // Store values for DistChord method 95 // Store values for DistChord method 94 // 96 // 95 fyInitial = G4ThreeVector( yIn[0], yIn[1] 97 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]); 96 fpInitial = G4ThreeVector( yIn[3], yIn[4] 98 fpInitial = G4ThreeVector( yIn[3], yIn[4], yIn[5]); 97 fyFinal = G4ThreeVector( yOut[0], yOut[1 99 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]); 98 } 100 } 99 101 100 // ------------------------------------------- 102 // --------------------------------------------------------------------------- 101 103 102 // Integrator for RK from G3 with evaluation o 104 // Integrator for RK from G3 with evaluation of error in solution and delta 103 // geometry based on naive similarity with the 105 // geometry based on naive similarity with the case of uniform magnetic field. 104 // B1[3] is input and is the first magnetic f 106 // B1[3] is input and is the first magnetic field values 105 // B2[3] is output and is the final magnetic f 107 // B2[3] is output and is the final magnetic field values. 106 // 108 // 107 void G4RKG3_Stepper::StepWithEst( const G4doub 109 void G4RKG3_Stepper::StepWithEst( const G4double*, 108 const G4doub 110 const G4double*, 109 G4doub 111 G4double, 110 G4doub 112 G4double*, 111 G4doub 113 G4double&, 112 G4doub 114 G4double&, 113 const G4doub 115 const G4double*, 114 G4doub 116 G4double* ) 115 117 116 { 118 { 117 G4Exception("G4RKG3_Stepper::StepWithEst()", 119 G4Exception("G4RKG3_Stepper::StepWithEst()", "GeomField0001", 118 FatalException, "Method no longe 120 FatalException, "Method no longer used."); 119 } 121 } 120 122 121 // ------------------------------------------- 123 // ----------------------------------------------------------------- 122 124 123 // Integrator RK Stepper from G3 with only two 125 // Integrator RK Stepper from G3 with only two field evaluation per Step. 124 // It is used in propagation initial Step by s 126 // It is used in propagation initial Step by small substeps after solution 125 // error and delta geometry considerations. B[ 127 // error and delta geometry considerations. B[3] is magnetic field which 126 // is passed from substep to substep. 128 // is passed from substep to substep. 127 // 129 // 128 void G4RKG3_Stepper::StepNoErr(const G4double 130 void G4RKG3_Stepper::StepNoErr(const G4double tIn[8], 129 const G4double 131 const G4double dydx[6], 130 G4double 132 G4double Step, 131 G4double 133 G4double tOut[8], 132 G4double 134 G4double B[3] ) 133 135 134 { 136 { 135 137 136 // Copy and edit the routine above, to dele 138 // Copy and edit the routine above, to delete alpha2, beta2, ... 137 // 139 // 138 G4double K1[7], K2[7], K3[7], K4[7]; 140 G4double K1[7], K2[7], K3[7], K4[7]; 139 G4double tTemp[8]={0.0}, yderiv[6]={0.0}; 141 G4double tTemp[8]={0.0}, yderiv[6]={0.0}; 140 142 141 // Need Momentum value to give correct valu 143 // Need Momentum value to give correct values to the coefficients in 142 // equation. Integration on unit velocity, 144 // equation. Integration on unit velocity, but tIn[3,4,5] is momentum 143 145 144 G4double mom, inverse_mom; 146 G4double mom, inverse_mom; 145 const G4double c1=0.5, c2=0.125, c3=1./6.; 147 const G4double c1=0.5, c2=0.125, c3=1./6.; 146 148 147 // Correction for momentum not a velocity 149 // Correction for momentum not a velocity 148 // Need the protection !!! must be not zero 150 // Need the protection !!! must be not zero 149 // 151 // 150 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4] 152 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]); 151 inverse_mom = 1./mom; 153 inverse_mom = 1./mom; 152 for(auto i=0; i<3; ++i) 154 for(auto i=0; i<3; ++i) 153 { 155 { 154 K1[i] = Step * dydx[i+3]*inverse_mom; 156 K1[i] = Step * dydx[i+3]*inverse_mom; 155 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*in 157 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*inverse_mom + c2*K1[i]) ; 156 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ; 158 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ; 157 } 159 } 158 160 159 GetEquationOfMotion()->EvaluateRhsReturnB(t 161 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; 160 162 161 for(auto i=0; i<3; ++i) 163 for(auto i=0; i<3; ++i) 162 { 164 { 163 K2[i] = Step * yderiv[i+3]*inverse_mom; 165 K2[i] = Step * yderiv[i+3]*inverse_mom; 164 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ; 166 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ; 165 } 167 } 166 168 167 // Given B, calculate yderiv ! 169 // Given B, calculate yderiv ! 168 // 170 // 169 GetEquationOfMotion()->EvaluateRhsGivenB(tT 171 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ; 170 172 171 for(auto i=0; i<3; ++i) 173 for(auto i=0; i<3; ++i) 172 { 174 { 173 K3[i] = Step * yderiv[i+3]*inverse_mom; 175 K3[i] = Step * yderiv[i+3]*inverse_mom; 174 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inver 176 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ; 175 tTemp[i+3] = tIn[i+3] + K3[i]*mom ; 177 tTemp[i+3] = tIn[i+3] + K3[i]*mom ; 176 } 178 } 177 179 178 // Calculates y-deriv(atives) & returns B t 180 // Calculates y-deriv(atives) & returns B too! 179 // 181 // 180 GetEquationOfMotion()->EvaluateRhsReturnB(t 182 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; 181 183 182 for(auto i=0; i<3; ++i) // Output tr 184 for(auto i=0; i<3; ++i) // Output trajectory vector 183 { 185 { 184 K4[i] = Step * yderiv[i+3]*inverse_mom; 186 K4[i] = Step * yderiv[i+3]*inverse_mom; 185 tOut[i] = tIn[i] + Step*(tIn[i+3]*invers 187 tOut[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom+ (K1[i]+K2[i]+K3[i])*c3) ; 186 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2 188 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ; 187 } 189 } 188 tOut[6] = tIn[6]; 190 tOut[6] = tIn[6]; 189 tOut[7] = tIn[7]; 191 tOut[7] = tIn[7]; 190 } 192 } 191 193 192 // ------------------------------------------- 194 // --------------------------------------------------------------------------- 193 195 194 G4double G4RKG3_Stepper::DistChord() const 196 G4double G4RKG3_Stepper::DistChord() const 195 { 197 { 196 // Soon: must check whether h/R > 2 pi !! 198 // Soon: must check whether h/R > 2 pi !! 197 // Method below is good only for < 2 pi 199 // Method below is good only for < 2 pi 198 200 199 G4double distChord,distLine; 201 G4double distChord,distLine; 200 202 201 if (fyInitial != fyFinal) 203 if (fyInitial != fyFinal) 202 { 204 { 203 distLine = G4LineSection::Distline(fyMid 205 distLine = G4LineSection::Distline(fyMidPoint,fyInitial,fyFinal); 204 distChord = distLine; 206 distChord = distLine; 205 } 207 } 206 else 208 else 207 { 209 { 208 distChord = (fyMidPoint-fyInitial).mag() 210 distChord = (fyMidPoint-fyInitial).mag(); 209 } 211 } 210 212 211 return distChord; 213 return distChord; 212 } 214 } 213 215