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

module mcint_m

  implicit none

  public::f

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

contains

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

  function f(x) result (y)

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

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

    y=1.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,b="")",advance="no")
  read(unit=*,fmt=*) a,b
  write(unit=*,fmt="(""   n="")",advance="no")
  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="("" s="",es14.7,"" ("",es14.7"")"")") &
       s,atan(b)-atan(a)

end program mcint

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