68 fInitialPoint = { y[0], y[1], y[2] };
75 const G4double S6 = hstep * one_sixth;
78 if (notEquals(momentum2, fMomentum2))
80 fMomentum = std::sqrt(momentum2);
81 fMomentum2 = momentum2;
82 fInverseMomentum = 1. / fMomentum;
83 fCoefficient = GetFCof() * fInverseMomentum;
88 fInverseMomentum * dydx[3],
89 fInverseMomentum * dydx[4],
90 fInverseMomentum * dydx[5]
95 y[0] + S5 * (dydx[0] + S4 * K1[0]),
96 y[1] + S5 * (dydx[1] + S4 * K1[1]),
97 y[2] + S5 * (dydx[2] + S4 * K1[2]),
101 GetFieldValue(p, field);
104 dydx[0] + S5 * K1[0],
105 dydx[1] + S5 * K1[1],
110 (A2[1] * field[2] - A2[2] * field[1]) * fCoefficient,
111 (A2[2] * field[0] - A2[0] * field[2]) * fCoefficient,
112 (A2[0] * field[1] - A2[1] * field[0]) * fCoefficient
115 fMidPoint = { p[0], p[1], p[2] };
119 dydx[0] + S5 * K2[0],
120 dydx[1] + S5 * K2[1],
125 (A3[1] * field[2] - A3[2] * field[1]) * fCoefficient,
126 (A3[2] * field[0] - A3[0] * field[2]) * fCoefficient,
127 (A3[0] * field[1] - A3[1] * field[0]) * fCoefficient
131 p[0] = y[0] + hstep * (dydx[0] + S5 * K3[0]);
132 p[1] = y[1] + hstep * (dydx[1] + S5 * K3[1]);
133 p[2] = y[2] + hstep * (dydx[2] + S5 * K3[2]);
135 GetFieldValue(p, field);
138 dydx[0] + hstep * K3[0],
139 dydx[1] + hstep * K3[1],
140 dydx[2] + hstep * K3[2]
144 (A4[1] * field[2] - A4[2] * field[1]) * fCoefficient,
145 (A4[2] * field[0] - A4[0] * field[2]) * fCoefficient,
146 (A4[0] * field[1] - A4[1] * field[0]) * fCoefficient
150 yOut[0] = y[0] + hstep * (dydx[0] + S6 * (K1[0] + K2[0] + K3[0]));
151 yOut[1] = y[1] + hstep * (dydx[1] + S6 * (K1[1] + K2[1] + K3[1]));
152 yOut[2] = y[2] + hstep * (dydx[2] + S6 * (K1[2] + K2[2] + K3[2]));
154 yOut[3] = dydx[0] + S6 * (K1[0] + K4[0] + 2. * (K2[0] + K3[0]));
155 yOut[4] = dydx[1] + S6 * (K1[1] + K4[1] + 2. * (K2[1] + K3[1]));
156 yOut[5] = dydx[2] + S6 * (K1[2] + K4[2] + 2. * (K2[2] + K3[2]));
161 fEndPoint = { yOut[0], yOut[1], yOut[2] };
164 yError[3] = hstep * std::fabs(K1[0] - K2[0] - K3[0] + K4[0]);
165 yError[4] = hstep * std::fabs(K1[1] - K2[1] - K3[1] + K4[1]);
166 yError[5] = hstep * std::fabs(K1[2] - K2[2] - K3[2] + K4[2]);
167 yError[0] = hstep * yError[3];
168 yError[1] = hstep * yError[4];
169 yError[2] = hstep * yError[5];
170 yError[3] *= fMomentum;
171 yError[4] *= fMomentum;
172 yError[5] *= fMomentum;