c===========================================================
c     bisect: Uses bisection to find approximate root
c     of f(x) on interval [xmin .. xmax].  Finds
c     root to (relative) tolerance 'xtol'.  Return code 
c     'rc' set to 0 on success, non-zero on failure and 
c     routine succeeds (by definition) as long as initial
c     interval *does* bracket at least one root.  Routine
c     performs tracing of algorithm (on stderr) if input
c     argument 'trace' is .true.
c===========================================================
      real*8 function bisect(f,xmin,xmax,xtol,trace,rc)

         implicit     none

         real*8       drelabs

         real*8       f
         external     f

         real*8       xmin,          xmax,         xtol 
         logical      trace
         integer      rc
c-----------------------------------------------------------
c        Other variables needed for search.
c-----------------------------------------------------------
         integer      mxiter
         parameter  ( mxiter = 50 )

         real*8       xlo,           dx,           sfxmid,        
     &                sgn
         integer      iter

c-----------------------------------------------------------
c        Check that input interval is specified correctly
c        and that it manifestly brackets at least one root: 
c        (i.e. the fcn changes sign).
c-----------------------------------------------------------
         if( xmax .le. xmin  .or.  
     &       f(xmin) * f(xmax) .gt. 0.0d0 ) then
            write(0,*) 'bisect:  Input interval is not '//
     &                 'bracketing'
            rc = 1
c-----------------------------------------------------------
c           Returned value is meaningless in this case,
c           but have to return *some* value.
c-----------------------------------------------------------
            bisect = xmin
            return
         end if

c-----------------------------------------------------------
c        Compute 'sgn' such that sgn * f(xmin) < 0, and 
c        intialize bracketing interval
c-----------------------------------------------------------
         sgn  = 1.0d0
         if( f(xmin) .le. 0.0d0 ) then
            sgn =  1.0d0
         else
            sgn = -1.0d0 
         end if
         xlo  = xmin
         dx   = xmax - xmin

c-----------------------------------------------------------
c        Bisection loop: continue until root found to     
c        specfied tolerance or until maximum number of 
c        iterations taken
c-----------------------------------------------------------
         do iter = 1 , mxiter
            bisect = xlo + 0.5d0 * dx
            if( trace ) then
               write(0,*) xlo, xlo + dx, f(bisect)
            end if
            if( sgn * f(bisect) .lt. 0.0d0 ) then
               xlo = bisect
            end if
            if( drelabs(dx,bisect,1.0d-6) .le. xtol ) go to 900
            dx = 0.5d0 * dx
         end do

900      continue
         rc = 0 
         if( trace ) write(0,*)
         return

      end

c-----------------------------------------------------------
c     drelabs:  Function useful for 'relativizing' quantity
c     being monitored for detection of convergence.  
c-----------------------------------------------------------
      real*8 function drelabs(dx,x,xfloor)

         implicit      none

         real*8        dx,     x,    xfloor

         if( abs(x) .lt. abs(xfloor) ) then
            drelabs = abs(dx)
         else
            drelabs = abs(dx/x)
         end if

         return

      end