134{
135
136
137
138 G4double K1[7], K2[7], K3[7], K4[7];
139 G4double tTemp[8]={0.0}, yderiv[6]={0.0};
140
141
142
143
145 const G4double c1=0.5, c2=0.125, c3=1./6.;
146
147
148
149
150 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
151 inverse_mom = 1./mom;
152 for(auto i=0; i<3; ++i)
153 {
154 K1[i] = Step * dydx[i+3]*inverse_mom;
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 ;
157 }
158
160
161 for(auto i=0; i<3; ++i)
162 {
163 K2[i] = Step * yderiv[i+3]*inverse_mom;
164 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;
165 }
166
167
168
170
171 for(auto i=0; i<3; ++i)
172 {
173 K3[i] = Step * yderiv[i+3]*inverse_mom;
174 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ;
175 tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
176 }
177
178
179
181
182 for(auto i=0; i<3; ++i)
183 {
184 K4[i] = Step * yderiv[i+3]*inverse_mom;
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[i] + 2*K3[i] +K4[i])*c3 ;
187 }
188 tOut[6] = tIn[6];
189 tOut[7] = tIn[7];
190}
virtual void EvaluateRhsGivenB(const G4double y[], const G4double B[3], G4double dydx[]) const =0
void EvaluateRhsReturnB(const G4double y[], G4double dydx[], G4double Field[]) const
G4EquationOfMotion * GetEquationOfMotion()