]>
Commit | Line | Data |
---|---|---|
a65a7e70 | 1 | /* |
2 | ||
3 | Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna | |
4 | amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru | |
5 | November. 2, 2006 | |
6 | ||
7 | */ | |
8 | //This equation solver class is taken from GEANT4 and modified!! | |
9 | ||
10 | #include <iostream> | |
11 | #include <TMath.h> | |
12 | ||
13 | template <class Function> | |
14 | Bool_t EquationSolver<Function>::Brent(Function& theFunction) { | |
15 | const Double_t precision = 3.0e-8; | |
16 | ||
17 | // Check the interval before start | |
18 | if(fA > fB || TMath::Abs(fA-fB) <= fTolerance) { | |
19 | std::cerr << "EquationSolver::Brent: The interval must be properly set." << std::endl; | |
20 | return false; | |
21 | } | |
22 | Double_t fa = theFunction(fA); | |
23 | Double_t fb = theFunction(fB); | |
24 | if(fa*fb > 0.0) { | |
25 | std::cerr << "EquationSolver::Brent: The interval must include a root." << std::endl; | |
26 | return false; | |
27 | } | |
28 | ||
29 | Double_t c = fB; | |
30 | Double_t fc = fb; | |
31 | Double_t d = 0.0; | |
32 | //Double_t fd = 0.0; | |
33 | Double_t e = 0.0; | |
34 | //Double_t fe = 0.0; | |
35 | ||
36 | for(Int_t i=0; i < fMaxIter; i++) { | |
37 | // Rename a,b,c and adjust bounding interval d | |
38 | if(fb*fc > 0.0) { | |
39 | c = fA; | |
40 | fc = fa; | |
41 | d = fB - fA; | |
42 | e = d; | |
43 | } | |
44 | if(TMath::Abs(fc) < TMath::Abs(fb)) { | |
45 | fA = fB; | |
46 | fB = c; | |
47 | c = fA; | |
48 | fa = fb; | |
49 | fb = fc; | |
50 | fc = fa; | |
51 | } | |
52 | Double_t Tol1 = 2.0*precision*TMath::Abs(fB) + 0.5*fTolerance; | |
53 | Double_t xm = 0.5*(c-fB); | |
54 | if(TMath::Abs(xm) <= Tol1 || fb == 0.0) { | |
55 | fRoot = fB; | |
56 | return true; | |
57 | } | |
58 | // Inverse quadratic interpolation | |
59 | if(TMath::Abs(e) >= Tol1 && TMath::Abs(fa) > TMath::Abs(fb)) { | |
60 | Double_t s = fb/fa; | |
61 | Double_t p = 0.0; | |
62 | Double_t q = 0.0; | |
63 | if(fA == c) { | |
64 | p = 2.0*xm*s; | |
65 | q = 1.0 - s; | |
66 | } | |
67 | else { | |
68 | q = fa/fc; | |
69 | Double_t r = fb/fc; | |
70 | p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0)); | |
71 | q = (q-1.0)*(r-1.0)*(s-1.0); | |
72 | } | |
73 | // Check bounds | |
74 | if(p > 0.0) q = -q; | |
75 | p = TMath::Abs(p); | |
76 | Double_t min1 = 3.0*xm*q-TMath::Abs(Tol1*q); | |
77 | Double_t min2 = TMath::Abs(e*q); | |
78 | if (2.0*p < TMath::Min(min1,min2)) { | |
79 | // Interpolation | |
80 | e = d; | |
81 | d = p/q; | |
82 | } | |
83 | else { | |
84 | // Bisection | |
85 | d = xm; | |
86 | e = d; | |
87 | } | |
88 | } | |
89 | else { | |
90 | // Bounds decreasing too slowly, use bisection | |
91 | d = xm; | |
92 | e = d; | |
93 | } | |
94 | // Move last guess to a | |
95 | fA = fB; | |
96 | fa = fb; | |
97 | if (TMath::Abs(d) > Tol1) fB += d; | |
98 | else { | |
99 | if (xm >= 0.0) fB += TMath::Abs(Tol1); | |
100 | else fB -= TMath::Abs(Tol1); | |
101 | } | |
102 | fb = theFunction(fB); | |
103 | } | |
104 | std::cerr << "EquationSolver::Brent: Number of iterations exceeded." << std::endl; | |
105 | return false; | |
106 | } | |
107 | ||
108 | template <class Function> | |
109 | void EquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) { | |
110 | if(TMath::Abs(Limit1-Limit2) <= fTolerance) { | |
111 | std::cerr << "EquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl; | |
112 | return; | |
113 | } | |
114 | if(Limit1 < Limit2) { | |
115 | fA = Limit1; | |
116 | fB = Limit2; | |
117 | } | |
118 | else { | |
119 | fA = Limit2; | |
120 | fB = Limit1; | |
121 | } | |
122 | return; | |
123 | } |