View Javadoc

1   package org.lcsim.recon.tracking.trfutil;
2   /**
3    * A linear root finder.
4    * @version 1.0
5    * @author  Norman A. Graf
6    */
7   public abstract class RootFindLinear implements RootFinder
8   {
9       // maximum # iterations
10      private int _max_iter;
11      
12      final public static int OK = 0;
13      final public static int INVALID_LIMITS=1;
14      final public static int INVALID_FUNCTION_CALL=2;
15      final public static int OUT_OF_RANGE=3;
16      final public static int TOO_MANY_ITERATIONS=4;
17      
18      /** Create a new instance of RootFindLinear */
19      public RootFindLinear()
20      {
21          _max_iter = 100;
22      }
23      
24      // Find a root.
25      public StatusDouble solve(double x1, double x2)
26      {
27          // Value to return for failure.
28          double xfail = 0.0;
29          
30          // Save limits.
31          // We require solution to be in the range (xmin,xmax).
32          double xmin = x1;
33          double xmax = x2;
34          
35          if ( xmin >= xmax ) return new StatusDouble(INVALID_LIMITS,xfail);
36          
37          // Evaluate the starting difference.
38          double adif = Math.abs(x2-x1);
39          
40          // Set the precision.
41          // Solution is said to be close if it changes by less than this value.
42          double xclose = 1.e-14 * Math.abs(x2-x1);
43          double yclose = 1.e-14;
44          
45          StatusDouble sd1 = evaluate(x1);
46          if ( sd1.status() != 0 ) return new StatusDouble(INVALID_FUNCTION_CALL,xfail);
47          double f1 = sd1.value();
48          StatusDouble sd2 = evaluate(x2);
49          if ( sd2.status() != 0 ) return new StatusDouble(INVALID_FUNCTION_CALL,xfail);
50          double f2 = sd2.value();
51          
52          // Loop.
53          for ( int count=0; count<_max_iter; ++count )
54          {
55              // interpolate for new solution
56              Assert.assertTrue( f1 != f2 );
57              if ( f1 == f2 ) return new StatusDouble(OUT_OF_RANGE,xfail);
58              double x0 = (f2*x1 - f1*x2) / (f2 - f1);
59              if ( x0 < xmin ) return new StatusDouble(OUT_OF_RANGE,xfail);
60              if ( x0 > xmax ) return new StatusDouble(OUT_OF_RANGE,xfail);
61              // Evaluate the function at the new solution.
62              StatusDouble sd0 = evaluate(x0);
63              if ( sd0.status() != 0 ) return new StatusDouble(INVALID_FUNCTION_CALL,xfail);
64              double f0 = sd0.value();
65              // Replace the most distant guess with x0.
66              double adif1 = Math.abs(x1-x0);
67              double adif2 = Math.abs(x2-x0);
68              if ( adif1 < adif2 )
69              {
70                  adif = adif1;
71                  x2 = x0;
72                  f2 = f0;
73              }
74              else
75              {
76                  adif = adif2;
77                  x1 = x0;
78                  f1 = f0;
79              }
80              // Return if solution is found.
81              if ( f0 == 0.0 ) return new StatusDouble(0,x0);
82              // If we are close in y, exit.
83              if ( Math.abs(f2-f1) < yclose ) return new StatusDouble(0,x0);
84              // If we are close in x, exit.
85              if ( adif < xclose ) return new StatusDouble(0,x0);
86          }
87          
88          // Too many iterations.
89          return new StatusDouble(TOO_MANY_ITERATIONS,xfail);
90          
91      }
92      
93  }