!**********************************************************************!
!
! File: bisect.f90 (05-Jun-2008)
!
! Bisektionsverfahren fuer
!
!   f(x) = cos(x)-x = 0
!
! a,b....Start- bzw. Bisektionsintervall
! xacc...Genauigkeit der Naeherungsloesung
!
!**********************************************************************!

module bisect_m

!**********************************************************************!

  implicit none

  private

  integer,parameter,public::double=selected_real_kind(15)

  public::f

contains

!**********************************************************************!

  function f(x) result (y)

!**********************************************************************!

    real(kind=double),intent(in)::x
    real(kind=double)::y

    y=cos(x)-x

    return

  end function f

end module bisect_m

!**********************************************************************!

program bisect

!**********************************************************************!

  use bisect_m

  implicit none

  integer,parameter::itmax=40

  integer::it
  real(kind=double)::a,b,fa,fb,fm,xacc,xm

  write(unit=*,fmt="(a)",advance="no") "  a,b="
  read(unit=*,fmt=*) a,b
  write(unit=*,fmt="(a)",advance="no") " xacc="
  read(unit=*,fmt=*) xacc

  fa=f(a)
  fb=f(b)

  if(fa*fb>0.0) then
    write(unit=*,fmt="(a)") &
         " bisect: f(a) und f(b) haben gleiches Vorzeichen"
    stop
  end if

  bisection: do it=1,itmax

    xm=a+0.5*(b-a)
    fm=f(xm)

    write(unit=*,fmt="(a,f12.8,a,es14.7,a,es14.7)") &
         " xm=", xm,"   |b-a|=",abs(b-a),"   f(xm)=",fm

    if(fa*fm<0.0) then   ! linkes Teilintervall
      b=xm
    else                 ! rechtes Teilintervall
      a=xm
      fa=fm
    end if

    if(abs(b-a)<xacc) then
      write(unit=*,fmt="(a,f12.8)") "  s=",a+0.5*(b-a)
      exit
    end if

  end do bisection

  if(it>itmax) then
    write(unit=*,fmt="(a)") " bisect: itmax ueberschritten"
  end if

end program bisect

!**********************************************************************!