131{
132
133
134
135
136
137
138 G4int nstp, i, no_warnings=0;
140
141#ifdef G4DEBUG_FIELD
143 static G4int nStpPr=50;
146#endif
147
151 G4bool succeeded =
true, lastStepSucceeded;
152
154
155 G4int noFullIntegr=0, noSmallIntegr = 0 ;
156 static G4int noGoodSteps =0 ;
157 const G4int nvar= fNoVars;
158
160
161
162
163 if( hstep <= 0.0 )
164 {
165 if(hstep==0.0)
166 {
167 std::ostringstream message;
168 message << "Proposed step is zero; hstep = " << hstep << " !";
171 return succeeded;
172 }
173 else
174 {
175 std::ostringstream message;
176 message <<
"Invalid run condition." <<
G4endl
177 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
178 << "Requested step cannot be negative! Aborting event.";
181 return false;
182 }
183 }
184
186
188 x1= startCurveLength;
189 x2= x1 + hstep;
190
191 if ( (hinitial > 0.0) && (hinitial < hstep)
192 && (hinitial > perMillion * hstep) )
193 {
194 h = hinitial;
195 }
196 else
197 {
198 h = hstep;
199 }
200
201 x = x1;
202
203 for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
204
206 nstp=1;
207
208 do
209 {
211
212#ifdef G4DEBUG_FIELD
214 for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
215 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
216 yFldTrkStart.SetCurveLength(x);
217#endif
218
219
220
221
223 fNoTotalSteps++;
224
225
226
227 if( h > fMinimumStep )
228 {
230
231 lastStepSucceeded= (hdid == h);
232#ifdef G4DEBUG_FIELD
233 if (dbg>2) {
234 PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp);
235 }
236#endif
237 }
238 else
239 {
242 G4double dchord_step, dyerr, dyerr_len;
243 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
244 yFldTrk.SetCurveLength( x );
245
246 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
247
248
249 yFldTrk.DumpToArray(y);
250
251#ifdef G4FLD_STATS
252 fNoSmallSteps++;
253 if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
254 fDyerrPos_smTot += dyerr_len;
255 fSumH_sm += h;
256 if (nstp<=1) { fNoInitialSmallSteps++; }
257#endif
258#ifdef G4DEBUG_FIELD
259 if (dbg>1)
260 {
261 if(fNoSmallSteps<2) {
PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
262 G4cout <<
"Another sub-min step, no " << fNoSmallSteps
263 <<
" of " << fNoTotalSteps <<
" this time " << nstp <<
G4endl;
265 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
266 << " epsilon= " << eps << " hstep= " << hstep
267 <<
" h= " << h <<
" hmin= " << fMinimumStep <<
G4endl;
268 }
269#endif
270 if( h == 0.0 )
271 {
274 "Integration Step became Zero!");
275 }
276 dyerr = dyerr_len / h;
277 hdid= h;
278 x += hdid;
279
280
282
283
284 lastStepSucceeded= (dyerr<= eps);
285 }
286
287 if (lastStepSucceeded) { noFullIntegr++; }
288 else { noSmallIntegr++; }
289
291
292#ifdef G4DEBUG_FIELD
293 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
294 {
295 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
297 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
298 << "hnext=" << std::setw(12) << hnext << " "
299 << "hstep=" << std::setw(12) << hstep << " (requested) "
301 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
302 }
303#endif
304
305
306 G4double endPointDist= (EndPos-StartPos).mag();
307 if ( endPointDist >= hdid*(1.+perMillion) )
308 {
309 fNoBadSteps++;
310
311
312
313 if ( endPointDist >= hdid*(1.+perThousand) )
314 {
315#ifdef G4DEBUG_FIELD
316 if (dbg)
317 {
319 G4cerr <<
" Total steps: bad " << fNoBadSteps
320 << " good " << noGoodSteps << " current h= " << hdid
322 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
323 }
324#endif
325 no_warnings++;
326 }
327 }
328 else
329 {
330 noGoodSteps ++;
331 }
332
333
334
335 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
336 {
337
338 lastStep = true;
339 }
340 else
341 {
342
343 if(std::fabs(hnext) <=
Hmin())
344 {
345#ifdef G4DEBUG_FIELD
346
347 if( (x < x2 * (1-eps) ) &&
348 (std::fabs(hstep) >
Hmin()) )
349 {
350 if(dbg>0)
351 {
353 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
354 }
355 no_warnings++;
356 }
357#endif
358
360 }
361 else
362 {
363 h = hnext;
364 }
365
366
367 if ( x+h > x2 )
368 {
369 h = x2 - x ;
370 }
371
372 if ( h == 0.0 )
373 {
374
375 lastStep = true;
376#ifdef G4DEBUG_FIELD
377 if (dbg>2)
378 {
379 int prec=
G4cout.precision(12);
380 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
382 << " Integration step 'h' became "
383 << h <<
" due to roundoff. " <<
G4endl
384 << " Calculated as difference of x2= "<< x2 << " and x=" << x
385 <<
" Forcing termination of advance." <<
G4endl;
387 }
388#endif
389 }
390 }
391 } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
392
393
394
395 succeeded= (x>=x2);
396
397 for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
398
399
402
403 if(nstp > fMaxNoSteps)
404 {
405 no_warnings++;
406 succeeded = false;
407#ifdef G4DEBUG_FIELD
408 if (dbg)
409 {
412 }
413#endif
414 }
415
416#ifdef G4DEBUG_FIELD
417 if( dbg && no_warnings )
418 {
419 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
421 }
422#endif
423
424 return succeeded;
425}
CLHEP::Hep3Vector G4ThreeVector
G4DLLIMPORT std::ostream G4cerr
G4double GetCurveLength() const
void SetCurveLength(G4double nCurve_s)
void DumpToArray(G4double valArr[ncompSVEC]) const
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)