!**********************************************************************!
!
! File: newton.f90 (05-Jun-2008)
!
! Newton-Verfahren fuer
!
!   f(x) = cos(x)-x = 0
!
! x0.....Startwert
! xacc...Genauigkeit der Naeherungsloesung
!
!**********************************************************************!

module newton_m

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

  implicit none

  private

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

  public::df,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
    
!**********************************************************************!

  function df(x) result(y)

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

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

    y=-sin(x)-1.0

    return

  end function df

end module newton_m

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

program newton

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

  use newton_m

  implicit none

  integer,parameter::itmax=20

  integer::it
  real(kind=double)::dx,x,x0,xacc,xn

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

  x=x0

  iteration: do it=1,itmax

    dx=-f(x)/df(x)
    xn=x+dx

    write(unit=*,fmt="(a,es14.7,a,es14.7)") &
         " x=",x,"   |xn-x|=",abs(dx)

    if(abs(dx)<xacc) then
      write(unit=*,fmt="(a,es14.7)") " s=",x
      exit iteration
    end if

    x=xn

  end do iteration

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

end program newton

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