/* --------------------------------------------------------------------- 
**  Figure 12.15:  A heat-transfer function. 
**  Figure 12.16:  Calling the bisection root-finder. 
**  Figure 12.18:  The bisect() function. 
**  Figure 12.19:  Iteration for the bisection method. 
** --------------------------------------------------------------------- 
*/
#include "tools.h"                                               

#define LIMIT 1.0e-10   /* How close to 0 a number must be to say it is 0.*/                                                   
#define EPS  .00001     /* Error tolerance. */           
#define PIPE_R 0.333    /* radius of pipe */

typedef enum { NO_ROOT, ROOT_OK } er_type;

/* ------------ Prototype for bisection root-finder --------------------- */      
er_type bisect( double r1, double r2, double* root ); 

/* ------------ function to be solved ----------------------------------- */    
double conductivity;                       /* Depends on insulation used. */
double
f( double r ){ 
    return (1.0 / conductivity) * log( r / PIPE_R ) + 1.0 / r - 2 / PIPE_R; 
}

/* ---------------------------------------------------------------------- */    
void swap( double* fp1, double* fp2 )
{   
    double swapper = *fp1;
    *fp1 = *fp2;
    *fp2 = swapper;
}   

/* ---------------------------------------------------------------------- */    
void main( void )                                                           
{   
    double r1, r2;          /* End points of interval */
    double root;            /* Value of root (to be computed) */           
    er_type success_code;   /* Result error code from bisect() */          

    banner();                                                         

    printf( "\n Enter thermal conductivity of insulation (BTU/h ft F): " );   
    scanf( "%lg", &conductivity );                                         

    do {         /* Enter search limits in correct order and valid range. */
        printf( " Enter interval limits (>= %.3f ft): ", PIPE_R );
        scanf ( "%lf%lf", &r1, &r2 );                                         
        if (r2 < r1) swap( &r1, &r2 ); 
        if (r1 < PIPE_R) puts( "\n Error: must be outside pipe radius." );
    } while (r1 < PIPE_R);
                   
    success_code = bisect( r1, r2, &root );              
    if (success_code == ROOT_OK)  {                 /* A root was found. */  
        printf( "\n The root is %g ft \n", root );  
        printf( " The thickness of the insulation is %.3g inches\n",
        (root  -PIPE_R) * 12 );
    }
    else printf( "\n No root found; try a different interval.\n" );
    bye();
}                                           

/* ------------ bisection method ---------------------------------------- */    
er_type                                                          
bisect(                                                            
    double xlow,               /* lower bound  */                 
    double xup,                /* upper bound  */                 
    double* rootp              /* address in which to store the root */
    )
{   
    double go_find ( double xlow, double fxlow, double xup ); 
    double fxlow, fxup;         /* function values at xlow and xup */                                             
    double result;              /* root of function */
    er_type answer;             /* success or failure code */

    answer = ROOT_OK;                                                   
    fxlow = f( xlow );   
    fxup = f( xup );
                                               
    if (fabs( fxlow ) < LIMIT) {        /* xlow is the actual root  */ 
        result = xlow;                                            
    }
    else if (fabs( fxup ) < LIMIT) {    /* xup is the actual root  */  
        result = xup;                                                
    }
    else if (fxup * fxlow < 0) {        /* true if they have opposite sign. */
        /* We know there is a root in the interval -- so find it. */   
        result = go_find( xlow, fxlow, xup );                           
    } 
    else answer = NO_ROOT;

    *rootp = result;
    return answer; 
}

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

    do {    
        xmid = .5 * (xup + xlow);                                         
        fxmid = f( xmid ); 
        // printf( "xlow=%g  fxlow=%g xmid=%g  fxmid=%g xup=%g\n", 
        //      xlow, fxlow, xmid, fxmid, xup );         /* For debugging */                                          
        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;                                                        
}