80{
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
108 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
109
110 G4double Energy = std::sqrt(momentum_mag_square + mass*mass);
111 G4double inverse_velocity = Energy*inv_momentum_magnitude/c_light;
112
113 G4double cof1 = ElectroMagCof*inv_momentum_magnitude;
115 G4double cof3 = inv_momentum_magnitude*mass;
116
117 dydx[0] = y[3]*inv_momentum_magnitude;
118 dydx[1] = y[4]*inv_momentum_magnitude;
119 dydx[2] = y[5]*inv_momentum_magnitude;
120
121 dydx[3] = 0.;
122 dydx[4] = 0.;
123 dydx[5] = 0.;
124
125 G4double field[18] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
126
127 field[0] = Field[0];
128 field[1] = Field[1];
129 field[2] = Field[2];
130
131
132
133 if (fBfield)
134 {
135 if (charge != 0.)
136 {
137 dydx[3] += cof1*(y[4]*field[2] - y[5]*field[1]);
138 dydx[4] += cof1*(y[5]*field[0] - y[3]*field[2]);
139 dydx[5] += cof1*(y[3]*field[1] - y[4]*field[0]);
140 }
141 }
142
143
144
145 if (!fBfield)
146 {
147 field[3] = Field[0];
148 field[4] = Field[1];
149 field[5] = Field[2];
150 }
151 else
152 {
153 field[3] = Field[3];
154 field[4] = Field[4];
155 field[5] = Field[5];
156 }
157
158 if (fEfield)
159 {
160 if (charge != 0.)
161 {
162 dydx[3] += cof1*cof2*field[3];
163 dydx[4] += cof1*cof2*field[4];
164 dydx[5] += cof1*cof2*field[5];
165 }
166 }
167
168
169
170 if (!fBfield && !fEfield)
171 {
172 field[6] = Field[0];
173 field[7] = Field[1];
174 field[8] = Field[2];
175 }
176 else
177 {
178 field[6] = Field[6];
179 field[7] = Field[7];
180 field[8] = Field[8];
181 }
182
183 if (fGfield)
184 {
185 if (mass > 0.)
186 {
187 dydx[3] += field[6]*cof2*cof3/c_light;
188 dydx[4] += field[7]*cof2*cof3/c_light;
189 dydx[5] += field[8]*cof2*cof3/c_light;
190 }
191 }
192
193
194
195 if (!fBfield && !fEfield && !fGfield)
196 {
197 field[9] = Field[0];
198 field[10] = Field[1];
199 field[11] = Field[2];
200 field[12] = Field[3];
201 field[13] = Field[4];
202 field[14] = Field[5];
203 field[15] = Field[6];
204 field[16] = Field[7];
205 field[17] = Field[8];
206 }
207 else
208 {
209 field[9] = Field[9];
210 field[10] = Field[10];
211 field[11] = Field[11];
212 field[12] = Field[12];
213 field[13] = Field[13];
214 field[14] = Field[14];
215 field[15] = Field[15];
216 field[16] = Field[16];
217 field[17] = Field[17];
218 }
219
220 if (fgradB)
221 {
222 if (magMoment != 0.)
223 {
224 dydx[3] += magMoment*(y[9]*field[ 9]+y[10]*field[10]+y[11]*field[11])
225 *inv_momentum_magnitude*Energy;
226 dydx[4] += magMoment*(y[9]*field[12]+y[10]*field[13]+y[11]*field[14])
227 *inv_momentum_magnitude*Energy;
228 dydx[5] += magMoment*(y[9]*field[15]+y[10]*field[16]+y[11]*field[17])
229 *inv_momentum_magnitude*Energy;
230 }
231 }
232
233 dydx[6] = 0.;
234
235
236
237 dydx[7] = inverse_velocity;
238
239 if (fNvar == 12)
240 {
241 dydx[ 8] = 0.;
242
243 dydx[ 9] = 0.;
244 dydx[10] = 0.;
245 dydx[11] = 0.;
246 }
247
248 if (fSpin)
249 {
251 if (fBfield)
252 {
254 BField = F;
255 }
256
258 if (fEfield)
259 {
261 EField = F;
262 }
263
264 EField /= c_light;
265
267 u *= inv_momentum_magnitude;
268
269 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
270 G4double ucb = (anomaly+1./gamma)/beta;
271 G4double uce = anomaly + 1./(gamma+1.);
272
274
276 if (charge == 0.) pcharge = 1.;
277 else pcharge = charge;
278
280 if (
Spin.mag2() != 0.)
281 {
282 if (fBfield)
283 {
284 dSpin =
285 pcharge*omegac*( ucb*(
Spin.cross(BField))-udb*(
Spin.cross(u)) );
286 }
287 if (fEfield)
288 {
289 dSpin -= pcharge*omegac*( uce*(u*(
Spin*EField) - EField*(Spin*u)) );
290
291
292
293 }
294 }
295
296 dydx[ 9] = dSpin.x();
297 dydx[10] = dSpin.y();
298 dydx[11] = dSpin.z();
299 }
300
301 return;
302}