source: trunk/source/geometry/magneticfield/src/G4ChordFinder.cc@ 852

Last change on this file since 852 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 22.8 KB
RevLine 
[831]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
[850]27// $Id: G4ChordFinder.cc,v 1.49 2008/07/15 14:02:06 japost Exp $
28// GEANT4 tag $Name: HEAD $
[831]29//
30//
31// 25.02.97 John Apostolakis, design and implimentation
32// 05.03.97 V. Grichine , style modification
33// -------------------------------------------------------------------
34
35#include "G4ChordFinder.hh"
36#include "G4MagneticField.hh"
37#include "G4Mag_UsualEqRhs.hh"
38#include "G4ClassicalRK4.hh"
39
40#include <iomanip>
41
42// ..........................................................................
43
44G4ChordFinder::G4ChordFinder(G4MagInt_Driver* pIntegrationDriver)
45 : fDefaultDeltaChord( 0.25 * mm ),
46 fDeltaChord( fDefaultDeltaChord ),
47 fAllocatedStepper(false),
48 fEquation(0),
49 fDriversStepper(0),
50 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
51 fMultipleRadius(15.0),
52 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
53 fStatsVerbose(0)
54{
55 // Simple constructor which does not create equation, ..
56 // fDeltaChord= fDefaultDeltaChord;
57 fIntgrDriver= pIntegrationDriver;
58 fAllocatedStepper= false;
59 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
60
61 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
62 // check the values and set the other parameters
63}
64
65// ..........................................................................
66
67G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
68 G4double stepMinimum,
69 G4MagIntegratorStepper* pItsStepper )
70 : fDefaultDeltaChord( 0.25 * mm ),
71 fDeltaChord( fDefaultDeltaChord ),
72 fAllocatedStepper(false),
73 fEquation(0),
74 fDriversStepper(0),
75 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
76 fMultipleRadius(15.0),
77 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
78 fStatsVerbose(0)
79{
80 // Construct the Chord Finder
81 // by creating in inverse order the Driver, the Stepper and EqRhs ...
82 G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
83 fEquation = pEquation;
84 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
85 // G4FieldTrack ??
86
87 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
88 // check the values and set the other parameters
89
90 // --->> Charge Q = 0
91 // --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
92
93 if( pItsStepper == 0 )
94 {
95 pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
96 fAllocatedStepper= true;
97 }
98 else
99 {
100 fAllocatedStepper= false;
101 }
102 fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
103 pItsStepper->GetNumberOfVariables() );
104}
105
106// ......................................................................
107
108void
109G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
110{
111 // Use -1.0 as request for Default.
112 if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
113 if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
114
115 // fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
116 // fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
117
118 if( fStatsVerbose ) {
119 G4cout << " ChordFnd> Trying to set fractions: "
120 << " first " << fFirstFraction
121 << " last " << fractLast
122 << " next " << fractNext
123 << " and multiple " << fMultipleRadius
124 << G4endl;
125 }
126
127 if( (fractLast > 0.0) && (fractLast <=1.0) )
128 { fFractionLast= fractLast; }
129 else
130 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
131 << " fraction Last = " << fractLast
132 << " must be 0 < fractionLast <= 1 " << G4endl;
133 if( (fractNext > 0.0) && (fractNext <1.0) )
134 { fFractionNextEstimate = fractNext; }
135 else
136 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
137 << " fraction Next = " << fractNext
138 << " must be 0 < fractionNext < 1 " << G4endl;
139}
140
141// ......................................................................
142
143G4ChordFinder::~G4ChordFinder()
144{
145 delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
146 if( fAllocatedStepper)
147 {
148 delete fDriversStepper;
149 } // fIntgrDriver->pIntStepper;}
150 delete fIntgrDriver;
151
152 if( fStatsVerbose ) { PrintStatistics(); }
153}
154
155void
156G4ChordFinder::PrintStatistics()
157{
158 // Print Statistics
159 G4cout << "G4ChordFinder statistics report: " << G4endl;
160 G4cout
161 << " No trials: " << fTotalNoTrials_FNC
162 << " No Calls: " << fNoCalls_FNC
163 << " Max-trial: " << fmaxTrials_FNC
164 << G4endl;
165 G4cout
166 << " Parameters: "
167 << " fFirstFraction " << fFirstFraction
168 << " fFractionLast " << fFractionLast
169 << " fFractionNextEstimate " << fFractionNextEstimate
170 << G4endl;
171}
172
173// ......................................................................
174
175G4double
176G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
177 G4double stepMax,
178 G4double epsStep,
179 const G4ThreeVector latestSafetyOrigin,
180 G4double latestSafetyRadius
181 )
182{
183 G4double stepPossible;
184 G4double dyErr;
185 G4FieldTrack yEnd( yCurrent);
186 G4double startCurveLen= yCurrent.GetCurveLength();
187 G4double nextStep;
188 // *************
189 stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
190 , latestSafetyOrigin, latestSafetyRadius
191 );
192 // *************
[850]193 // G4cout<<"Exit Find Next Chord Err= "<<dyErr<<" eps= "<<epsStep<<"stepPos= "<<stepPossible<<G4endl;
[831]194 G4bool good_advance;
[850]195
[831]196 if ( dyErr < epsStep * stepPossible )
[850]197
198 { //G4cout<<"err comparison = "<<dyErr<<" eps= "<<epsStep<<G4endl;
[831]199 // Accept this accuracy.
200 yCurrent = yEnd;
201 good_advance = true;
202 }
203 else
[850]204 {
205 // G4cout<<"Entering Accurate Advance"<<G4endl;
[831]206 // Advance more accurately to "end of chord"
207 // ***************
208 good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
209 // ***************
210 if ( ! good_advance ){
211 // In this case the driver could not do the full distance
212 stepPossible= yCurrent.GetCurveLength()-startCurveLen;
213 }
214 }
215
216#ifdef G4DEBUG_FIELD
[850]217 //G4cout << "Exiting FindNextChord Limited with:" << G4endl
218 // << " yCurrent: " << yCurrent<< G4endl;
[831]219#endif
220
221 return stepPossible;
222}
223
224// #define TEST_CHORD_PRINT 1
225
226// ............................................................................
227
228G4double
[850]229G4ChordFinder::FindNextChord( const G4FieldTrack& yStart,
[831]230 G4double stepMax,
231 G4FieldTrack& yEnd, // Endpoint
232 G4double& dyErrPos, // Error of endpoint
233 G4double epsStep,
234 G4double* pStepForAccuracy,
235 const G4ThreeVector, // latestSafetyOrigin,
236 G4double // latestSafetyRadius
237 )
238 // Returns Length of Step taken
239{
[850]240 //G4cout<<"Inter Find Next Chord with step="<<stepMax<<G4endl;
[831]241 // G4int stepRKnumber=0;
242 G4FieldTrack yCurrent= yStart;
243 G4double stepTrial, stepForAccuracy;
244 G4double dydx[G4FieldTrack::ncompSVEC];
245
246 // 1.) Try to "leap" to end of interval
247 // 2.) Evaluate if resulting chord gives d_chord that is good enough.
248 // 2a.) If d_chord is not good enough, find one that is.
249
250 G4bool validEndPoint= false;
251 G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
252
253 fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
[850]254 //for (G4int i=0;i<G4FieldTrack::ncompSVEC;i++){
255 // dydx[i]=0.;
256 //}
[831]257 G4int noTrials=0;
258 const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
259
260 stepTrial = std::min( stepMax,
261 safetyFactor * fLastStepEstimate_Unconstrained );
262
263 G4double newStepEst_Uncons= 0.0;
264 do
265 {
266 G4double stepForChord;
267 yCurrent = yStart; // Always start from initial point
[850]268
[831]269 // ************
270 fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
271 dChordStep, dyErrPos);
272 // ************
273
[850]274 // G4cout<<"AfterQuickAdv step="<<stepTrial<<" dC="<<dChordStep<<" yErr="<<dyErrPos<<G4endl;
275
276 // We check whether the criterion is met here.
277 validEndPoint = AcceptableMissDist(dChordStep);
278 // if(validEndPoint){G4cout<<"validEndPoint"<<fDeltaChord<<G4endl;}
279 // else{G4cout<<"No__validEndPoint"<<G4endl;}
[831]280
[850]281
282 //&& (dyErrPos < eps) ;
283 // validEndPoint = AcceptableMissDist(dChordStep) && (dyErrPos < epsStep) ;
284
[831]285 lastStepLength = stepTrial;
286
287 // This method estimates to step size for a good chord.
288 stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
289
290 if( ! validEndPoint ) {
291 if( stepTrial<=0.0 )
292 stepTrial = stepForChord;
293 else if (stepForChord <= stepTrial)
294 // Reduce by a fraction, possibly up to 20%
295 stepTrial = std::min( stepForChord,
296 fFractionLast * stepTrial);
297 else
298 stepTrial *= 0.1;
299
300 // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
301 }
302 // #ifdef TEST_CHORD_PRINT
303 // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
304 // #endif
305
306 noTrials++;
307 }
308 while( ! validEndPoint ); // End of do-while RKD
309
310 if( newStepEst_Uncons > 0.0 ){
311 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
312 }
313
314 AccumulateStatistics( noTrials );
315
316 // stepOfLastGoodChord = stepTrial;
317
318 if( pStepForAccuracy ){
319 // Calculate the step size required for accuracy, if it is needed
320 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
321 if( dyErr_relative > 1.0 ) {
322 stepForAccuracy =
323 fIntgrDriver->ComputeNewStepSize( dyErr_relative,
324 lastStepLength );
325 }else{
326 stepForAccuracy = 0.0; // Convention to show step was ok
327 }
328 *pStepForAccuracy = stepForAccuracy;
329 }
330
331#ifdef TEST_CHORD_PRINT
332 static int dbg=0;
333 if( dbg )
334 G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
335 << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
336#endif
[850]337 //G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
338 // << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
[831]339 yEnd= yCurrent;
340 return stepTrial;
341}
342
343// ----------------------------------------------------------------------------
344#if 0
345// #ifdef G4VERBOSE
346if( dbg ) {
347 G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
348 G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
349}
350#endif
351// ----------------------------------------------------------------------------
352
353// ...........................................................................
354
355G4double G4ChordFinder::NewStep(G4double stepTrialOld,
356 G4double dChordStep, // Curr. dchord achieved
357 G4double& stepEstimate_Unconstrained )
358//
359// Is called to estimate the next step size, even for successful steps,
360// in order to predict an accurate 'chord-sensitive' first step
361// which is likely to assist in more performant 'stepping'.
362//
363
364{
365 G4double stepTrial;
366 static G4double lastStepTrial = 1., lastDchordStep= 1.;
367
368#if 1
369 // const G4double threshold = 1.21, multiplier = 0.9;
370 // 0.9 < 1 / std::sqrt(1.21)
371
372 if (dChordStep > 0.0)
373 {
374 stepEstimate_Unconstrained = stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
375 // stepTrial = 0.98 * stepEstimate_Unconstrained;
376 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
377 }
378 else
379 {
380 // Should not update the Unconstrained Step estimate: incorrect!
381 stepTrial = stepTrialOld * 2.;
382 }
383
384 // if ( dChordStep < threshold * fDeltaChord ){
385 // stepTrial= stepTrialOld * multiplier;
386 // }
387 if( stepTrial <= 0.001 * stepTrialOld)
388 {
389 if ( dChordStep > 1000.0 * fDeltaChord ){
390 stepTrial= stepTrialOld * 0.03;
391 }else{
392 if ( dChordStep > 100. * fDeltaChord ){
393 stepTrial= stepTrialOld * 0.1;
394 }else{
395 // Try halving the length until dChordStep OK
396 stepTrial= stepTrialOld * 0.5;
397 }
398 }
399 }else if (stepTrial > 1000.0 * stepTrialOld)
400 {
401 stepTrial= 1000.0 * stepTrialOld;
402 }
403
404 if( stepTrial == 0.0 ){
405 stepTrial= 0.000001;
406 }
407
408 lastStepTrial = stepTrialOld;
409 lastDchordStep= dChordStep;
410#else
411 if ( dChordStep > 1000. * fDeltaChord ){
412 stepTrial= stepTrialOld * 0.03;
413 }else{
414 if ( dChordStep > 100. * fDeltaChord ){
415 stepTrial= stepTrialOld * 0.1;
416 }else{
417 // Keep halving the length until dChordStep OK
418 stepTrial= stepTrialOld * 0.5;
419 }
420 }
421#endif
422
423 // A more sophisticated chord-finder could figure out a better
424 // stepTrial, from dChordStep and the required d_geometry
425 // eg
426 // Calculate R, r_helix (eg at orig point)
427 // if( stepTrial < 2 pi R )
428 // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
429 // else
430 // ??
431
432 return stepTrial;
433}
434
[850]435// ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
436// Use Brent Algorithm(or InvParabolic) when it possible.
437// Given a starting curve point A (CurveA_PointVelocity),
438// curve point B (CurveB_PointVelocity), a point E which is (generally)
439// not on the curve and a point F which is on the curve(first approximation)
440// From this information find new point S on the curve closer to point E.
441// While advancing towards S utilise eps_step
[831]442// as a measure of the relative accuracy of each Step.
443G4FieldTrack
[850]444G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
445 const G4FieldTrack& CurveB_PointVelocity,
446 const G4ThreeVector& CurrentE_Point,
447 const G4ThreeVector& CurrentF_Point,
448 const G4ThreeVector& PointG,
449 G4bool first, G4double eps_step)
450{
451
452
453 G4FieldTrack EndPoint( CurveA_PointVelocity);
454 G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
455 G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
456 G4double xa,xb,xc,ya,yb,yc;
457
458//InverseParabolic
459//AF Intersects (First Part of Curve)
460 if(first){
461 xa=0.;
462 ya=(PointG-Point_A).mag();
463 xb=(Point_A-CurrentF_Point).mag();
464 yb=-(PointG-CurrentF_Point).mag();
465 xc=(Point_A-Point_B).mag();
466 yc=-(CurrentE_Point-Point_B).mag();
467 }
468 else{
469 xa=0.;
470 ya=(Point_A-PointG).mag();
471 xb=(Point_B-Point_A).mag();
472 yb=-(PointG-Point_B).mag();
473 xc=-(Point_A-CurrentF_Point).mag();
474 yc=-(Point_A-CurrentE_Point).mag();
475
476 }
477 const G4double tolerance= 1.e-12;
478 if(ya<=tolerance||std::abs(yc)<=tolerance){
479 ; //What to do for the moment return the same point as in begin
480 //Then PropagatorInField will take care
481 }
482 else{
483
484 G4double test_step =InvParabolic(xa,ya,xb,yb,xc,yc);
485 G4double curve=std::abs(EndPoint.GetCurveLength()-CurveB_PointVelocity.GetCurveLength());
486 G4double dist= (EndPoint.GetPosition()-Point_B).mag();
487 if(test_step<=0) { test_step=0.1*xb;}
488 if(test_step>=xb){ test_step=0.5*xb;}
489
490
491 if(curve*(1.+eps_step)<dist){
492 test_step=0.5*dist;
493 }
494
495 G4bool goodAdvance;
496 goodAdvance=
497 fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
498 // ***************
499
500 #ifdef G4DEBUG_FIELD
501 G4cout<<"G4ChordFinder:: test-step ShF="<<test_step<<" EndPoint="<<EndPoint<<G4endl;
502 // Test Track
503 G4FieldTrack TestTrack( CurveA_PointVelocity);
504 TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
505 CurveB_PointVelocity,
506 CurrentE_Point,
507 eps_step );
508 G4cout.precision(14);
509 G4cout<<"G4ChordFinder:: BrentApprox="<<EndPoint<<G4endl;
510 G4cout<<"G4ChordFinder::LinearApprox="<<TestTrack<<G4endl;
511 #endif
512 }
513return EndPoint;
514}
515
516G4FieldTrack
[831]517G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
518 const G4FieldTrack& CurveB_PointVelocity,
519 const G4ThreeVector& CurrentE_Point,
520 G4double eps_step)
521{
522 // 1st implementation:
523 // if r=|AE|/|AB|, and s=true path lenght (AB)
524 // return the point that is r*s along the curve!
[850]525 /////////////////////////////
526 //
527 //2st implementation : Inverse Parabolic Extrapolation by D.C.Williams
528 //
529 // Uses InvParabolic (xa,ya,xb,yb,xc,yc)
[831]530
[850]531 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
[831]532
533 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
534 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
535
536 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
537 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
538
539 G4double ABdist= ChordAB_Vector.mag();
540 G4double curve_length; // A curve length of AB
541 G4double AE_fraction;
542
543 curve_length= CurveB_PointVelocity.GetCurveLength()
544 - CurveA_PointVelocity.GetCurveLength();
545
546 // const
547 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
548 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
549#ifdef G4DEBUG_FIELD
550 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
551 << G4endl
552 << " The two points are further apart than the curve length "
553 << G4endl
554 << " Dist = " << ABdist
555 << " curve length = " << curve_length
556 << " relativeDiff = " << (curve_length-ABdist)/ABdist
557 << G4endl;
558 if( curve_length < ABdist * (1. - 10*eps_step) ) {
559 G4cerr << " ERROR: the size of the above difference"
560 << " exceeds allowed limits. Aborting." << G4endl;
561 G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
562 FatalException, "Unphysical curve length.");
563 }
564#endif
565 // Take default corrective action:
566 // --> adjust the maximum curve length.
567 // NOTE: this case only happens for relatively straight paths.
568 curve_length = ABdist;
569 }
570
571 G4double new_st_length;
572
573 if ( ABdist > 0.0 ){
574 AE_fraction = ChordAE_Vector.mag() / ABdist;
575 }else{
576 AE_fraction = 0.5; // Guess .. ?;
577#ifdef G4DEBUG_FIELD
578 G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
579 << " A and B are the same point!" << G4endl
580 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
581 << G4endl;
582#endif
583 }
584
585 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ){
586#ifdef G4DEBUG_FIELD
587 G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
588 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
589 << " AE_fraction = " << AE_fraction << G4endl
590 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
591 << " Chord AB length = " << ABdist << G4endl << G4endl;
592 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
593 << G4endl << " Otherwise it is an error. " << G4endl ;
594#endif
595 // This course can now result if B has been re-evaluated,
596 // without E being recomputed (1 July 99)
597 // In this case this is not a "real error" - but it undesired
598 // and we cope with it by a default corrective action ...
599 AE_fraction = 0.5; // Default value
600 }
601
602 new_st_length= AE_fraction * curve_length;
603
604 G4bool good_advance;
605 if ( AE_fraction > 0.0 ) {
606 good_advance =
607 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
608 new_st_length,
609 eps_step ); // Relative accuracy
610 // In this case it does not matter if it cannot advance the full distance
611 }
612
613 // If there was a memory of the step_length actually require at the start
614 // of the integration Step, this could be re-used ...
[850]615 G4cout.precision(14);
616
617 // G4cout<<"G4ChordFinder::LinearApprox="<<Current_PointVelocity<<G4endl;
[831]618 return Current_PointVelocity;
619}
620
621void
622G4ChordFinder::TestChordPrint( G4int noTrials,
623 G4int lastStepTrial,
624 G4double dChordStep,
625 G4double nextStepTrial )
626{
627 G4int oldprec= G4cout.precision(5);
628 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
629 << " this_step= " << std::setw(10) << lastStepTrial;
630 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
631 G4cout.precision(8);
632 }else{ G4cout.precision(6); }
633 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
634 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
635 else { G4cout << " d-"; }
636 G4cout.precision(5);
637 G4cout << " new_step= " << std::setw(10)
638 << fLastStepEstimate_Unconstrained
639 << " new_step_constr= " << std::setw(10)
640 << lastStepTrial << G4endl;
641 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
642 G4cout.precision(oldprec);
643}
Note: See TracBrowser for help on using the repository browser.