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

module func_m

  implicit none

  public::f

contains

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

  function f(x) result(y)

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

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

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

    return

  end function f

end module func_m

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

program trapez

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

  use func_m

  implicit none

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

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

  dx=(b-a)/n

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

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

end program trapez

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