program rootnr
c
c    written by: A Busy Code Developer, January 2003
c    modified 1:
c
c  PURPOSE:  This program finds the root of a given function that 
c  lies between the input bounds using the hybrid Newton-Raphson/
c  bisection method.
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       integer       niter
       real*8        xl      , xr      , root    , maxerr , err
       real*8        fofx    , dfdx
       real          cputot  , etime
c
       real          tdummy  (   2)
c
       external      binrph  , fofx    , dfdx
c
       data          maxerr
     1             / 5.0d-9 /
c
c-----------------------------------------------------------------------
c
c      Start up CPU time counter.
c
       cputot = etime ( tdummy )
c
c      Ask for initial limits/guesses.
c
       write(6, 2010)
       read (5,    *) xl, xr
c
c      Compute and report root, if any.  "niter" is returned as zero if 
c  no root is found.  Otherwise, it is the number of iterations taken to
c  converge on the root.
c
       call binrph ( xl, xr, root, fofx, dfdx, maxerr, err, niter )
c
       if (niter .gt. 0) then
         write(6, 2020) root, err, niter
       else
         write(6, 2030) xl, xr
       endif
c
c      End CPU time counter, and report cpu time used.
c
       cputot = etime ( tdummy ) - cputot
       write(6, 2040) cputot
c
       stop
c
2010   format('ROOTNR  : Enter xL and xR such that f(xL)*f(xR) < 0.'  
     1       ,'  Thus, xL < root < xR.')
2020   format('ROOTNR  : Root = ',1pg22.15,' +/- ',g9.2
     1       ,' after ',i3,' iterations.')
2030   format('ROOTNR  : No root found in domain (',1pg12.5
     1       ,', ',g12.5,').')
2040   format('ROOTNR  : Total cpu usage for this run is ',1pg12.5
     1       ,' seconds.')
c
       end
c
c=======================================================================
c
       subroutine binrph ( xl, xr, xn, f, fprime, maxerr, err, iter )
c
c    written by: A Busy Code Developer, January 2003
c    modified 1: 
c
c  PURPOSE:  This routine performs the hybrid bisection/Newton-Raphson
c  method to find the root of a function.
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       integer       iter    , maxiter
       real*8        xl      , xr      , xn      , xnp1    , fxl
     1             , fxr     , fxn     , dfxn    , err     , ferr
     2             , maxerr  , small
       real*8        f       , fprime
c
       data          maxiter , small   
     1             / 30    , 1.0d-99 /
c
c-----------------------------------------------------------------------
c
c      Initialise variables.
c
       iter = 0
       fxl  = f(xl)
       fxr  = f(xr)
c
c      Data validation:  Make sure root lies in between the left- and 
c  right-hand value.  If not, return iter=0 to the calling routine.
c
       if (fxl*fxr .gt. 0.0d0) then
         write (6, 2010) xl, xr
         return
       endif
c
c      Starting values for "xn", "fxn", and "dfxn".
c
       if (abs(fxl) .lt. abs(fxr)) then
         xn  = xl
         fxn = fxl
       else
         xn  = xr
         fxn = fxr
       endif
       dfxn = fprime(xn)
c
c-----------------------------------------------------------------------
c----- Top of hybrid loop. ---------------------------------------------
c-----------------------------------------------------------------------
c
10     continue
         iter = iter + 1
c
c      Evaluate the next Newton-Raphson guess.
c
         xnp1 = xn - fxn / ( dfxn + small )
c
c      If NR step is not justified, take a bisection step instead.
c
         if ( (xl .lt. xnp1) .and. (xnp1 .lt. xr) ) then
           write (6, 2020) iter, xnp1
         else
           xnp1 = 0.5d0 * ( xl + xr )
           write (6, 2030) iter, xnp1
         endif
c
c      Evaluate error estimate and set new guess.
c
         err = abs ( xnp1 - xn )
         xn  = xnp1
c
c      If the fractional error of "xn" is less than "maxerr", the root 
c  has been found---execution is returned to calling routine.
c
         ferr = err / ( abs (xn) + small )
         if ( ferr .le. maxerr ) go to 20
c
c      If the maximum number of iterations has been reached, abort.
c
         if ( iter .ge. maxiter ) then
           write (6, 2040) maxiter, xn, err
           stop
         endif
c
c      Prepare for next iteration.
c
         fxn  = f     (xn)
         dfxn = fprime(xn)
c
c      Determine which half of the interval (xl,xr) the root is in, and 
c  reset either xr or xl as appropriate.
c
         if ( fxl * fxn .le. 0.0d0 ) then
           xr  = xn
           fxr = fxn
         else
           xl  = xn
           fxl = fxn
         endif
         go to 10
20     continue
c
c-----------------------------------------------------------------------
c----- Bottom of hybrid loop. ------------------------------------------
c-----------------------------------------------------------------------
c
2010   format('BINRPH  : No unique root lies between the specified'
     1       ,' limits:',/
     2       ,'BINRPH  : xl = ', 1pg12.5,', and xr = ',g12.5)
2020   format('BINRPH  : Iteration: ',i3,' - NEWTON/RAPHSON: root ='
     1       ,1pg22.15)
2030   format('BINRPH  : Iteration: ',i3,' - BISECTION     : root ='
     1       ,1pg22.15)
2040   format('BINRPH  : Maximum number of iterations ',i9
     1       ,' exceeded.',/
     2       ,'BINRPH  : Best guess at root is ',1pg22.15
     3       ,' +/- ',g9.2,'.  Abort!')
c
       return
       end
c
c=======================================================================
c
       function fofx ( x )
c
c    written by: A Busy Code Developer, January 2003
c    modified 1: 
c
c  PURPOSE:  This routine returns the value of cos(x) - x.
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       real*8        x       , fofx
c
c-----------------------------------------------------------------------
c
       fofx = cos(x) - x
c
       return
       end
c
c=======================================================================
c
       function dfdx ( x )
c
c    written by: A Busy Code Developer, January 2003
c    modified 1: 
c
c  PURPOSE:  This routine returns the value of the derivative of fofx 
c  evaluated at x.
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       real*8        x       , dfdx
c
c-----------------------------------------------------------------------
c
       dfdx =-sin(x) - 1.0d0
c
       return
       end