100{
101
102
103
104
105
106
109
110#ifdef G4DEBUG_FIELD
111 G4int no_warnings = 0;
112 static G4int dbg = 1;
113 static G4int nStpPr = 50;
116#endif
117
122
124
125 const G4int nvar = fNoVars;
126
128
129
130
131 if( hstep <= 0.0 )
132 {
133 if( hstep == 0.0 )
134 {
135 std::ostringstream message;
136 message << "Proposed step is zero; hstep = " << hstep << " !";
139 return succeeded;
140 }
141 else
142 {
143 std::ostringstream message;
144 message <<
"Invalid run condition." <<
G4endl
145 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
146 << "Requested step cannot be negative! Aborting event.";
149 return false;
150 }
151 }
152
154
156 x1= startCurveLength;
157 x2= x1 + hstep;
158
159 if ( (hinitial > 0.0) && (hinitial < hstep)
160 && (hinitial > perMillion * hstep) )
161 {
162 h = hinitial;
163 }
164 else
165 {
166 h = hstep;
167 }
168
169 x = x1;
170
171 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
172
174 nstp = 1;
175
176 do
177 {
179
180#ifdef G4DEBUG_FIELD
182 for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
183 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
184 yFldTrkStart.SetCurveLength(x);
185#endif
186
188 ++fNoTotalSteps;
189
190
191
192 if( h > fMinimumStep )
193 {
195
196#ifdef G4DEBUG_FIELD
197 if (dbg>2)
198 {
199
201 }
202#endif
203 }
204 else
205 {
208 G4double dchord_step, dyerr, dyerr_len;
209 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
210 yFldTrk.SetCurveLength( x );
211
212 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
213
214
215 yFldTrk.DumpToArray(y);
216
217#ifdef G4FLD_STATS
218 ++fNoSmallSteps;
219 if ( dyerr_len > fDyerr_max ) { fDyerr_max = dyerr_len; }
220 fDyerrPos_smTot += dyerr_len;
221 fSumH_sm += h;
222 if (nstp==1) { ++fNoInitialSmallSteps; }
223#endif
224#ifdef G4DEBUG_FIELD
225 if (dbg>1)
226 {
227 if(fNoSmallSteps<2) {
PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
228 G4cout <<
"Another sub-min step, no " << fNoSmallSteps
229 <<
" of " << fNoTotalSteps <<
" this time " << nstp <<
G4endl;
231 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
232 << " epsilon= " << eps << " hstep= " << hstep
233 <<
" h= " << h <<
" hmin= " << fMinimumStep <<
G4endl;
234 }
235#endif
236 if( h == 0.0 )
237 {
240 "Integration Step became Zero!");
241 }
242 dyerr = dyerr_len / h;
243 hdid = h;
244 x += hdid;
245
246
248 }
249
251
252#ifdef G4DEBUG_FIELD
253 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
254 {
255 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
257 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
258 << "hnext=" << std::setw(12) << hnext << " "
259 << "hstep=" << std::setw(12) << hstep << " (requested) "
261 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
262 }
263#endif
264
265
266 G4double endPointDist= (EndPos-StartPos).mag();
267 if ( endPointDist >= hdid*(1.+perMillion) )
268 {
269 ++fNoBadSteps;
270
271
272
273 if ( endPointDist >= hdid*(1.+perThousand) )
274 {
275#ifdef G4DEBUG_FIELD
276 if (dbg)
277 {
279 G4cerr <<
" Total steps: bad " << fNoBadSteps
280 <<
" current h= " << hdid <<
G4endl;
281 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
282 }
283 ++no_warnings;
284#endif
285 }
286 }
287
288
289 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
290 {
291
292 lastStep = true;
293 }
294 else
295 {
296
297 if(std::fabs(hnext) <=
Hmin())
298 {
299#ifdef G4DEBUG_FIELD
300
301 if( (x < x2 * (1-eps) ) &&
302 (std::fabs(hstep) >
Hmin()) )
303 {
304 if(dbg>0)
305 {
307 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
308 }
309 ++no_warnings;
310 }
311#endif
312
314 }
315 else
316 {
317 h = hnext;
318 }
319
320
321 if ( x+h > x2 )
322 {
323 h = x2 - x ;
324 }
325
326 if ( h == 0.0 )
327 {
328
329 lastStep = true;
330#ifdef G4DEBUG_FIELD
331 if (dbg>2)
332 {
333 int prec=
G4cout.precision(12);
334 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
336 << " Integration step 'h' became "
337 << h <<
" due to roundoff. " <<
G4endl
338 << " Calculated as difference of x2= "<< x2 << " and x=" << x
339 <<
" Forcing termination of advance." <<
G4endl;
341 }
342#endif
343 }
344 }
345 } while ( ((++nstp)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
346
347
348
349
350
351 succeeded = (x>=x2);
352
353 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
354
355
358
359 if(nstp > fMaxNoSteps)
360 {
361 succeeded = false;
362#ifdef G4DEBUG_FIELD
363 ++no_warnings;
364 if (dbg)
365 {
368 }
369#endif
370 }
371
372#ifdef G4DEBUG_FIELD
373 if( dbg && no_warnings )
374 {
375 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
377 }
378#endif
379
380 return succeeded;
381}
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