100{
101
102
103
104
105
106
107 G4int nstp, i, no_warnings = 0;
109
110#ifdef G4DEBUG_FIELD
111 static G4int dbg = 1;
112 static G4int nStpPr = 50;
115#endif
116
120 G4bool succeeded =
true, lastStepSucceeded;
121
123
124 G4int noFullIntegr = 0, noSmallIntegr = 0;
126 const G4int nvar = fNoVars;
127
129
130
131
132 if( hstep <= 0.0 )
133 {
134 if( hstep == 0.0 )
135 {
136 std::ostringstream message;
137 message << "Proposed step is zero; hstep = " << hstep << " !";
140 return succeeded;
141 }
142 else
143 {
144 std::ostringstream message;
145 message <<
"Invalid run condition." <<
G4endl
146 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
147 << "Requested step cannot be negative! Aborting event.";
150 return false;
151 }
152 }
153
155
157 x1= startCurveLength;
158 x2= x1 + hstep;
159
160 if ( (hinitial > 0.0) && (hinitial < hstep)
161 && (hinitial > perMillion * hstep) )
162 {
163 h = hinitial;
164 }
165 else
166 {
167 h = hstep;
168 }
169
170 x = x1;
171
172 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
173
175 nstp = 1;
176
177 do
178 {
180
181#ifdef G4DEBUG_FIELD
183 for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
184 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
185 yFldTrkStart.SetCurveLength(x);
186#endif
187
189 ++fNoTotalSteps;
190
191
192
193 if( h > fMinimumStep )
194 {
196
197 lastStepSucceeded = (hdid == h);
198#ifdef G4DEBUG_FIELD
199 if (dbg>2)
200 {
201
203 }
204#endif
205 }
206 else
207 {
210 G4double dchord_step, dyerr, dyerr_len;
211 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
212 yFldTrk.SetCurveLength( x );
213
214 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
215
216
217 yFldTrk.DumpToArray(y);
218
219#ifdef G4FLD_STATS
220 ++fNoSmallSteps;
221 if ( dyerr_len > fDyerr_max ) { fDyerr_max = dyerr_len; }
222 fDyerrPos_smTot += dyerr_len;
223 fSumH_sm += h;
224 if (nstp==1) { ++fNoInitialSmallSteps; }
225#endif
226#ifdef G4DEBUG_FIELD
227 if (dbg>1)
228 {
229 if(fNoSmallSteps<2) {
PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
230 G4cout <<
"Another sub-min step, no " << fNoSmallSteps
231 <<
" of " << fNoTotalSteps <<
" this time " << nstp <<
G4endl;
233 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
234 << " epsilon= " << eps << " hstep= " << hstep
235 <<
" h= " << h <<
" hmin= " << fMinimumStep <<
G4endl;
236 }
237#endif
238 if( h == 0.0 )
239 {
242 "Integration Step became Zero!");
243 }
244 dyerr = dyerr_len / h;
245 hdid = h;
246 x += hdid;
247
248
250
251
252 lastStepSucceeded = (dyerr<= eps);
253 }
254
255 if (lastStepSucceeded) { ++noFullIntegr; }
256 else { ++noSmallIntegr; }
257
259
260#ifdef G4DEBUG_FIELD
261 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
262 {
263 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
265 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
266 << "hnext=" << std::setw(12) << hnext << " "
267 << "hstep=" << std::setw(12) << hstep << " (requested) "
269 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
270 }
271#endif
272
273
274 G4double endPointDist= (EndPos-StartPos).mag();
275 if ( endPointDist >= hdid*(1.+perMillion) )
276 {
277 ++fNoBadSteps;
278
279
280
281 if ( endPointDist >= hdid*(1.+perThousand) )
282 {
283#ifdef G4DEBUG_FIELD
284 if (dbg)
285 {
287 G4cerr <<
" Total steps: bad " << fNoBadSteps
288 << " good " << noGoodSteps << " current h= " << hdid
290 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
291 }
292#endif
293 ++no_warnings;
294 }
295 }
296 else
297 {
298 ++noGoodSteps;
299 }
300
301
302
303 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
304 {
305
306 lastStep = true;
307 }
308 else
309 {
310
311 if(std::fabs(hnext) <=
Hmin())
312 {
313#ifdef G4DEBUG_FIELD
314
315 if( (x < x2 * (1-eps) ) &&
316 (std::fabs(hstep) >
Hmin()) )
317 {
318 if(dbg>0)
319 {
321 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
322 }
323 ++no_warnings;
324 }
325#endif
326
328 }
329 else
330 {
331 h = hnext;
332 }
333
334
335 if ( x+h > x2 )
336 {
337 h = x2 - x ;
338 }
339
340 if ( h == 0.0 )
341 {
342
343 lastStep = true;
344#ifdef G4DEBUG_FIELD
345 if (dbg>2)
346 {
347 int prec=
G4cout.precision(12);
348 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
350 << " Integration step 'h' became "
351 << h <<
" due to roundoff. " <<
G4endl
352 << " Calculated as difference of x2= "<< x2 << " and x=" << x
353 <<
" Forcing termination of advance." <<
G4endl;
355 }
356#endif
357 }
358 }
359 } while ( ((++nstp)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
360
361
362
363
364
365 succeeded = (x>=x2);
366
367 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
368
369
372
373 if(nstp > fMaxNoSteps)
374 {
375 ++no_warnings;
376 succeeded = false;
377#ifdef G4DEBUG_FIELD
378 if (dbg)
379 {
382 }
383#endif
384 }
385
386#ifdef G4DEBUG_FIELD
387 if( dbg && no_warnings )
388 {
389 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
391 }
392#endif
393
394 return succeeded;
395}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cerr
static void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, unsigned int subStepNo, unsigned int noIntegrationVariables)
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)
virtual G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr) override
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)
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent) override
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
void RightHandSide(const G4double y[], G4double dydx[]) const