!**********************************************************************!
!
! File: mcint.f90 (16-Jun-2001)
!                 (23-Jun-2006)
!                 (06-Jun-2008)
!
! Monte Carlo-Integration der Funktion
!
!   f(x) = 4/(1+x**2)
!
! ueber dem Interval [a,b]
!
! a,b...Intervallgrenzen
! n.....Anzahl der Stuetzpunkte
!
!**********************************************************************!

module mcint_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=4.0/(1.0+x**2)

  end function f

end module mcint_m

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

program mcint

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

  use mcint_m

  implicit none

  integer::i,n
  real(kind=double)::a,b,ran,s,x

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

  s=0.0

  do i=1,n

    call random_number(ran)
    x=a+(b-a)*ran

    s=s+f(x)

  end do

  s=(b-a)*s/n

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

end program mcint

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