!**********************************************************************!
!
! File: histo.f90 (16-Jun-2001)
!                 (15-Jun-2005)
!
! Mittelwert, Varianz und Wahrscheinlichkeitsdichte der gleichverteilten
! Zufallsvariablen x auf dem Intervall [0,1]
!
! nbins.........Anzahl der Histogrammkanaele
! n.............Stichprobenumfang
! "histo.dat"...Ausgabefile fuer Histogramm
!
!**********************************************************************!

program histo

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

  implicit none

  integer,parameter::double=selected_real_kind(15)
  integer,parameter::iout=2
  real(kind=double),parameter::a=0.0,b=1.0

  integer::i,k,n,nbins
  integer,dimension(:),allocatable::ihist
  real(kind=double)::ax,axx,ran,x

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

  allocate(ihist(0:nbins-1))

  ax=0.0
  axx=0.0
  ihist=0

  do i=1,n
    call random_number(ran)
    x=a+(b-a)*ran
    ax=ax+x
    axx=axx+x**2
    k=min(int((x-a)/(b-a)*nbins),nbins-1)
    ihist(k)=ihist(k)+1
  end do

  write(unit=*,fmt="(""    <x>="",es14.7)") ax/n
  write(unit=*,fmt="("" Var(x)="",es14.7)") (axx-ax**2/n)/(n-1)

  open(unit=iout,file="histo.dat",status="replace",action="write", &
       form="formatted")

  do k=0,nbins-1
    x=a+(k+0.5)*(b-a)/nbins
    write(unit=iout,fmt="(2(tr1,es14.7))") &
         x,real(ihist(k))*nbins/(n*(b-a))
  end do

  close(unit=iout)

end program histo

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