!------------------------------------------------------------------------
! This module contains the bisection algorithm to find roots of functions.
! Other methods such as Newton or secant can be added. 
! The bisection method can be used as follows:
! usage: bisect
!	(
!	f: function which zeroes are approximated
!	left, right: interval in which a zero will be found
!	error: maximum error in the obtained result
!       zero:  the value at which the function is zero within the required tolerance
!       ifail:   error handling
!                =0  successful search
!                =1  root is not bracketed f(left)*f(right)>0
!                =2  invalid tolerance     tolerance<0
!                =3  invalid interval      left>right
!	)
! Note: the borders of the interval are copied into local variables to ensure
!	they aren't changed. In this way, they can still be used in the calling
!	function.
!       The number of iterations is calculated by the well known formula:
!
!			n=(b-a)/log 2
!

module rootfinding
implicit none
contains
	subroutine bisect(f,left,right,tol,zero,ifail)
		real*8::zero
		real*8,intent(in)::left,right,tol
		real*8,external::f
                real*8,intent(out)::zero
		real*8::ileft,iright  ! local copies of left and right
		integer::i,iterations  ! i is a counter for loops
					! iterations is the number of intertions
		integer::ifail          !error indicator
		real*8::try		! middle of the interval

! the following errors cause immediate return since
! they would cause a floating exception in the log
                if(f(left)*f(right)>0) then
                 ifail=1
                 return
                endif

                if(tol < 0)            then
                 ifail=2
                 return
                endif

                if(left > right)       then
                 ifail=3
                 return
                endif
		ileft=left

		iright=right
		iterations=ceiling(log(tol/(iright-ileft))/log(0.5))
		do i=1,iterations
			try=(ileft+iright)/2
			if(f(ileft)*f(try)<0)then
				iright=try
			else
				ileft=try
			end if
		end do
		zero=(ileft+iright)/2
	end subroutine bisect

end module rootfinding