!**********************************************************************!
!
! File: trapez.f90 (15-Mar-2001)
!                  (24-Mar-2006)
!                  (19-Mar-2008)
!
! Numerische Integration der Funktion
!
!   f(x) = 4/(1+x**2)
!
! mit der Trapezregel
!
! a,b...Intervallgrenzen
! n.....Anzahl der Stuetzpunkte
! s.....Wert des Integrals
!
!**********************************************************************!

module trapez_m

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

  implicit none

  private

  public::f

contains

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

  function f(x) result(y)

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

    real,intent(in)::x
    real::y

    y=4.0/(1.0+x**2)

    return

  end function f

end module trapez_m

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

program trapez

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

  use trapez_m

  implicit none

  integer::i,n
  real::a,b,dx,s,x

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

  dx=(b-a)/n

  s=0.5*(f(a)+f(b))
  do i=1,n-1
    x=a+i*dx
    s=s+f(x)
  end do
  s=dx*s

  write(unit=*,fmt="(a,es14.7,a,es14.7,a)") " s=",s, &
       " (",4.0*(atan(b)-atan(a)),")"

end program trapez

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