/* ---------------------------------------------------------------------- 
** Use the bisection method to find a root of a function named f().
** Figure 23.9:  The bisect() function.    
** Figure 23.10: The go_find() function.   
** ---------------------------------------------------------------------- 
*/
#include "tools.h"
#include "root.h"       /* Interface between roots module and main. */
#include "ffam.h"       /* Prototype for function to solve.   */

extern bool verbose;    /* Imported from main. */ 

static double go_find( double xlow, double fxlow, double xup ); 

/* ----------------------------------------------------------------------- 
** Bisection method   
*/    
er_type                                                          
bisect( double xlow, double xup, double *rootp ) 
{   
    double fxlow, fxup;         /* Function values at xlow and xup. */                 
    er_type answer = ROOT_OK;   /* Default result is successful search. */                                     
    double swap;                /* To interchange interval end points. */ 

    if (xlow > xup) { swap = xlow; xlow = xup; xup = swap; }
                                               
    fxlow = f( xlow );  fxup = f( xup );  
    if (verbose) {
        printf( "    x1 = %g, f(x1) = %g\n", xlow, fxlow );                                      
        printf( "    x2 = %g, f(x2) = %g\n", xup, fxup );        
    }
                                    
    if (fabs( fxlow ) < LIMIT) {     /* xlow is the actual root  */ 
        *rootp = xlow;                                            
    }
    else if (fabs( fxup ) < LIMIT) { /* xup is the actual root  */  
            *rootp = xup;                                                
    }
    else if (fxup * fxlow < 0) { /* true if they have opposite sign. */
            /* We know a root is in the interval -- so find it. */   
            *rootp = go_find( xlow, fxlow, xup );                           
    } 
    else answer = NO_ROOT;
    return answer; 
}

/* ----------------------------------------------------------------------- 
** Iteration for the bisection method.   
*/
double 
go_find( double xlow, double fxlow, double xup ) 
{ 
    double xmid;     /* midpoint of interval */                       
    double fxmid;    /* value of function at midpoint of interval */ 

    do {             /* Iterate until interval is too small. */
        xmid = .5 * (xup + xlow);                                         
        fxmid = f(xmid);
        if (verbose) printf( "    xmid = %g, f(xmid) = %g\n", xmid, fxmid ); 
                                               
        if (fabs( fxmid ) < LIMIT) break;   /* quit; we are close enough! */ 

        else if(fxmid * fxlow < 0) {        /* root is between xlow and xmid */
            xup = xmid;                     /* use xmid as new upper bound */
        }             
        else  {                             /* root is between xmid and xup */  
            xlow = xmid;                    /* use xmid as new lower bound */ 
            fxlow = fxmid;                                                 
        }                    
    } while ((xup - xlow) > EPS);   
    return xmid;                                                        
}