91 if( (std::fabs(R_1) < 1
e-10)||(Bmag<1
e-12) )
103 Bnorm = (1.0/Bmag)*Bfld;
107 B_x_P = Bnorm.
cross(initTangent);
111 B_d_P = Bnorm.
dot(initTangent);
113 vpar = B_d_P * Bnorm;
114 vperp= initTangent - vpar;
116 B_v_P = std::sqrt( 1 - B_d_P * B_d_P);
124 if( std::fabs(Theta) > approc_limit )
126 SinT = std::sin(Theta);
127 CosT = std::cos(Theta);
134 SinT = Theta - 1.0/6.0 * Theta3;
135 CosT = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4;
142 positionMove = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar;
143 endTangent = CosT * vperp + SinT * B_x_P + vpar;
147 yHelix[0] = yIn[0] + positionMove.
x();
148 yHelix[1] = yIn[1] + positionMove.
y();
149 yHelix[2] = yIn[2] + positionMove.
z();
150 yHelix[3] = velocityVal * endTangent.
x();
151 yHelix[4] = velocityVal * endTangent.
y();
152 yHelix[5] = velocityVal * endTangent.
z();
158 SinT2 = 2.0 * SinT * CosT;
159 CosT2 = 1.0 - 2.0 * SinT * SinT;
160 endTangent = (CosT2 * vperp + SinT2 * B_x_P + vpar);
161 positionMove = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar;
163 yHelix2[0] = yIn[0] + positionMove.
x();
164 yHelix2[1] = yIn[1] + positionMove.
y();
165 yHelix2[2] = yIn[2] + positionMove.
z();
166 yHelix2[3] = velocityVal * endTangent.
x();
167 yHelix2[4] = velocityVal * endTangent.
y();
168 yHelix2[5] = velocityVal * endTangent.
z();
194 const G4int nvar = 6;
204 for(
G4int i=0; i<nvar; ++i)
226 for(
G4int i=0; i<nvar; ++i)
228 yErr[i] = yOut[i] - yTemp[i] ;