63 : fDefaultDeltaChord( 0.25 *
mm ),
64 fDeltaChord( fDefaultDeltaChord ),
65 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
66 fMultipleRadius(15.0),
68 fRegularStepperOwned(nullptr),
70 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
88 : fDefaultDeltaChord( 0.25 *
mm ),
89 fDeltaChord( fDefaultDeltaChord ),
90 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
91 fMultipleRadius(15.0),
95 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
101 const char* NewFSALStepperName =
"G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
102 using RegularStepperType =
109 const char* RegularStepperName =
"G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded stepper";
113 G4bool forceFSALstepper=
false;
114 G4bool recallFSALflag = useFSALstepper;
115 useFSALstepper = forceFSALstepper || useFSALstepper;
118 G4cout <<
"G4ChordFinder 2nd Constructor called. " <<
G4endl;
120 G4cout <<
" useFSAL stepper= " << useFSALstepper
121 <<
" (request = " << recallFSALflag
122 <<
" force FSAL = " << forceFSALstepper <<
" )" <<
G4endl;
139 G4bool errorInStepperCreation =
false;
143 if( pItsStepper !=
nullptr )
149 else if ( !useFSALstepper )
152 auto regularStepper =
new RegularStepperType(pEquation);
159 if( regularStepper ==
nullptr )
161 message <<
"Stepper instantiation FAILED." <<
G4endl;
162 message <<
"G4ChordFinder: Attempted to instantiate "
163 << RegularStepperName <<
" type stepper " <<
G4endl;
166 errorInStepperCreation =
true;
173 regularStepper->GetNumberOfVariables() );
182 message <<
"Using G4IntegrationDriver with "
183 << RegularStepperName <<
" type stepper " <<
G4endl;
184 message <<
"Driver instantiation FAILED." <<
G4endl;
192 auto fsalStepper=
new NewFsalStepperType(pEquation);
196 if( fsalStepper ==
nullptr )
198 message <<
"Stepper instantiation FAILED." <<
G4endl;
199 message <<
"Attempted to instantiate "
200 << NewFSALStepperName <<
" type stepper " <<
G4endl;
203 errorInStepperCreation =
true;
210 fsalStepper->GetNumberOfVariables() );
215 message <<
"Using G4FSALIntegrationDriver with stepper type: "
216 << NewFSALStepperName <<
G4endl;
217 message <<
"Integration Driver instantiation FAILED." <<
G4endl;
233 if( errorInStepperCreation || (
fIntgrDriver ==
nullptr ))
235 std::ostringstream errmsg;
237 if( errorInStepperCreation )
239 errmsg <<
"ERROR> Failure to create Stepper object." << G4endl
240 <<
" --------------------------------" <<
G4endl;
244 errmsg <<
"ERROR> Failure to create Integration-Driver object." << G4endl
245 <<
" -------------------------------------------" <<
G4endl;
247 const std::string BoolName[2]= {
"False",
"True" };
248 errmsg <<
" Configuration: (constructor arguments) " << G4endl
249 <<
" provided Stepper = " << pItsStepper << G4endl
250 <<
" use FSAL stepper = " << BoolName[useFSALstepper]
251 <<
" (request = " << BoolName[recallFSALflag]
252 <<
" force FSAL = " << BoolName[forceFSALstepper] <<
" )" <<
G4endl;
253 errmsg << message.str();
254 errmsg <<
"Aborting.";
255 G4Exception(
"G4ChordFinder::G4ChordFinder() - constructor 2",
259 assert( ( pItsStepper !=
nullptr )
286 if( fractLast == -1.0 ) fractLast = 1.0;
287 if( fractNext == -1.0 ) fractNext = 0.98;
294 G4cout <<
" ChordFnd> Trying to set fractions: "
296 <<
" last " << fractLast
297 <<
" next " << fractNext
302 if( (fractLast > 0.0) && (fractLast <=1.0) )
309 message <<
"Invalid fraction Last = " << fractLast
310 <<
"; must be 0 < fractionLast <= 1 ";
311 G4Exception(
"G4ChordFinder::SetFractions_Last_Next()",
314 if( (fractNext > 0.0) && (fractNext <1.0) )
321 message <<
"Invalid fraction Next = " << fractNext
322 <<
"; must be 0 < fractionNext < 1 ";
323 G4Exception(
"G4ChordFinder::SetFractions_Last_Next()",
344 stepPossible=
FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
345 &nextStep, latestSafetyOrigin, latestSafetyRadius);
349 if ( dyErr < epsStep * stepPossible )
362 if ( ! good_advance )
388 G4double stepTrial, stepForAccuracy;
395 G4bool validEndPoint=
false;
396 G4double dChordStep, lastStepLength;
400 unsigned int noTrials=0;
401 const unsigned int maxTrials= 75;
414 dChordStep, dyErrPos);
420 lastStepLength = stepTrial;
423 stepForChord =
NewStep(stepTrial, dChordStep, newStepEst_Uncons );
425 if( ! validEndPoint )
429 stepTrial = stepForChord;
431 else if (stepForChord <= stepTrial)
443 while( (! validEndPoint) && (noTrials < maxTrials) );
446 if( noTrials >= maxTrials )
449 message <<
"Exceeded maximum number of trials= " << maxTrials <<
G4endl
450 <<
"Current sagita dist= " << dChordStep <<
G4endl
451 <<
"Step sizes (actual and proposed): " <<
G4endl
452 <<
"Last trial = " << lastStepLength <<
G4endl
453 <<
"Next trial = " << stepTrial <<
G4endl
454 <<
"Proposed for chord = " << stepForChord <<
G4endl
456 G4Exception(
"G4ChordFinder::FindNextChord()",
"GeomField0003",
460 if( newStepEst_Uncons > 0.0 )
462 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
467 if( pStepForAccuracy )
471 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
472 if( dyErr_relative > 1.0 )
479 stepForAccuracy = 0.0;
481 *pStepForAccuracy = stepForAccuracy;
493 G4double& stepEstimate_Unconstrained )
503 if (dChordStep > 0.0)
505 stepEstimate_Unconstrained =
506 stepTrialOld*std::sqrt(
fDeltaChord / dChordStep );
512 stepTrial = stepTrialOld * 2.;
515 if( stepTrial <= 0.001 * stepTrialOld)
519 stepTrial= stepTrialOld * 0.03;
525 stepTrial= stepTrialOld * 0.1;
529 stepTrial= stepTrialOld * 0.5;
533 else if (stepTrial > 1000.0 * stepTrialOld)
535 stepTrial= 1000.0 * stepTrialOld;
538 if( stepTrial == 0.0 )
547 stepTrial= stepTrialOld * 0.03;
553 stepTrial= stepTrialOld * 0.1;
557 stepTrial= stepTrialOld * 0.5;
597 if(!first){EndPoint= ApproxCurveV;}
610 ya=(PointG-Point_A).mag();
611 xb=(Point_A-CurrentF_Point).mag();
612 yb=-(PointG-CurrentF_Point).mag();
613 xc=(Point_A-Point_B).mag();
614 yc=-(CurrentE_Point-Point_B).mag();
619 ya=(Point_A-CurrentE_Point).mag();
620 xb=(Point_A-CurrentF_Point).mag();
621 yb=(PointG-CurrentF_Point).mag();
622 xc=(Point_A-Point_B).mag();
623 yc=-(Point_B-PointG).mag();
628 CurrentE_Point, eps_step);
634 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
650 test_step=(test_step-xb);
653 xb=(CurrentF_Point-Point_B).mag();
656 if(test_step<=0) { test_step=0.1*xb; }
657 if(test_step>=xb) { test_step=0.5*xb; }
658 if(test_step>=curve){ test_step=0.5*curve; }
660 if(curve*(1.+eps_step)<xb)
670 G4cout <<
"G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
671 << test_step <<
" EndPoint = " << EndPoint <<
G4endl;
677 CurveB_PointVelocity,
678 CurrentE_Point, eps_step );
680 G4cout <<
"G4ChordFinder::BrentApprox = " << EndPoint <<
G4endl;
681 G4cout <<
"G4ChordFinder::LinearApprox= " << TestTrack <<
G4endl;
699 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
715 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
718 G4cerr <<
" Warning in G4ChordFinder::ApproxCurvePoint: "
720 <<
" The two points are further apart than the curve length "
722 <<
" Dist = " << ABdist
723 <<
" curve length = " << curve_length
724 <<
" relativeDiff = " << (curve_length-ABdist)/ABdist
726 if( curve_length < ABdist * (1. - 10*eps_step) )
729 message <<
"Unphysical curve length." <<
G4endl
730 <<
"The size of the above difference exceeds allowed limits."
733 G4Exception(
"G4ChordFinder::ApproxCurvePointV()",
"GeomField0003",
746 AE_fraction = ChordAE_Vector.
mag() / ABdist;
752 G4cout <<
"Warning in G4ChordFinder::ApproxCurvePointV():"
753 <<
" A and B are the same point!" <<
G4endl
754 <<
" Chord AB length = " << ChordAE_Vector.
mag() <<
G4endl
759 if( (AE_fraction> 1.0 +
perMillion) || (AE_fraction< 0.) )
762 G4cerr <<
" G4ChordFinder::ApproxCurvePointV() - Warning:"
763 <<
" Anomalous condition:AE > AB or AE/AB <= 0 " <<
G4endl
764 <<
" AE_fraction = " << AE_fraction <<
G4endl
765 <<
" Chord AE length = " << ChordAE_Vector.
mag() <<
G4endl
767 G4cerr <<
" OK if this condition occurs after a recalculation of 'B'"
768 << G4endl <<
" Otherwise it is an error. " <<
G4endl ;
778 new_st_length= AE_fraction * curve_length;
780 if ( AE_fraction > 0.0 )
783 new_st_length, eps_step );
793 return Current_PointVelocity;
804 G4cout <<
"G4ChordFinder statistics report: " <<
G4endl;
827 G4cout <<
" ChF/fnc: notrial " << std::setw( 3) << noTrials
828 <<
" this_step= " << std::setw(10) << lastStepTrial;
829 if( std::fabs( (dChordStep /
fDeltaChord) - 1.0 ) < 0.001 )
837 G4cout <<
" dChordStep= " << std::setw(12) << dChordStep;
841 G4cout <<
" new_step= " << std::setw(10)
843 <<
" new_step_constr= " << std::setw(10)
844 << lastStepTrial <<
G4endl;
845 G4cout <<
" nextStepTrial = " << std::setw(10) << nextStepTrial <<
G4endl;
846 G4cout.precision(oldprec);
G4bool AcceptableMissDist(G4double dChordStep) const
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector &latestSafetyOrigin, G4double lasestSafetyRadius)
virtual G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErr, G4double epsStep, G4double *pNextStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4bool QuickAdvance(G4FieldTrack &track, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)=0
static constexpr double mm
G4double NewStep(G4double stepTrialOld, G4double dChordStep, G4double &stepEstimate_Unconstrained)
void message(RunManager *runmanager)
static constexpr double perMillion
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector ¤tEPoint, const G4ThreeVector ¤tFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
G4ThreeVector GetPosition() const
G4MagIntegratorStepper * fNewFSALStepperOwned
G4VIntegrationDriver * fIntgrDriver
G4double InvParabolic(const G4double xa, const G4double ya, const G4double xb, const G4double yb, const G4double xc, const G4double yc)
virtual void PrintStatistics()
G4double fLastStepEstimate_Unconstrained
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)=0
G4GLOB_DLL std::ostream G4cerr
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void TestChordPrint(G4int noTrials, G4int lastStepTrial, G4double dChordStep, G4double nextStepTrial)
G4MagIntegratorStepper * fRegularStepperOwned
G4GLOB_DLL std::ostream G4cout
G4EquationOfMotion * fEquation
G4int GetNumberOfVariables() const
G4double fFractionNextEstimate
virtual G4bool AccurateAdvance(G4FieldTrack &track, G4double hstep, G4double eps, G4double hinitial=0)=0
G4double GetCurveLength() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetFractions_Last_Next(G4double fractLast=0.90, G4double fractNext=0.95)
void AccumulateStatistics(G4int noTrials)
G4ChordFinder(G4VIntegrationDriver *pIntegrationDriver)