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