/// \ingroup newmat
///@{

/// \file solution.h
/// One dimensional solve routine.

#include "myexcept.h"

#ifdef use_namespace
namespace RBD_COMMON {
#endif


// Solve the equation f(x)=y for x where f is a monotone continuous
// function of x
// Essentially Brent s method

// You need to derive a class from R1_R1 and override "operator()"
// with the function you want to solve.
// Use an object from this class in OneDimSolve

class R1_R1
{
   // the prototype for a Real function of a Real variable
   // you need to derive your function from this one and put in your
   // function for operator() at least. You probably also want to set up a
   // constructor to put in additional parameter values (e.g. that will not
   // vary during a solve)

protected:
   Real x;                             // Current x value
   bool xSet;                          // true if a value assigned to x

public:
   Real minX, maxX;                    // range of value x
   bool minXinf, maxXinf;              // true if these are infinite
   R1_R1() : xSet(false), minXinf(true), maxXinf(true) {}
   virtual Real operator()() = 0;      // function value at current x
                                       // set current x
   virtual void Set(Real X);           // set x, check OK
   Real operator()(Real X) { Set(X); return operator()(); }
                                       // set x, return value
   virtual bool IsValid(Real X);
   operator Real();                    // implicit conversion
   virtual ~R1_R1() {}                 // keep gnu happy
};

class SolutionException : public BaseException
{
public:
   static unsigned long Select;
   SolutionException(const char* a_what = 0);
};

class OneDimSolve
{
   R1_R1& function;                     // reference to the function
   Real accX;                           // accuracy in X direction
   Real accY;                           // accuracy in Y direction
   int lim;                             // maximum number of iterations

public:
   OneDimSolve(R1_R1& f, Real AccY = 0.0001, Real AccX = 0.0)
      : function(f), accX(AccX), accY(AccY) {}
                       // f is an R1_R1 function
   Real Solve(Real Y, Real X, Real Dev, int Lim=100);
                       // Solve for x in Y=f(x)
                       // X is the initial trial value of x
                       // X+Dev is the second trial value
                       // program returns a value of x such that
                       // |Y-f(x)| <= accY or |f.inv(Y)-x| <= accX

private:
   Real x[3], y[3];                         // Trial values of X and Y
   int L,C,U,Last;                          // Locations of trial values
   int vpol, hpol;                          // polarities
   Real YY;                                 // target value
   int i;
   void LookAt(int);                        // get new value of function
   bool Finish;                             // true if LookAt finds conv.
   bool Captured;                           // true when target surrounded
   void VFlip();
   void HFlip();
   void Flip();
   void State(int I, int J, int K);
   void Linear(int, int, int);
   void Quadratic(int, int, int);
};


#ifdef use_namespace
}
#endif

// body file: solution.cpp


///@}