module const

  implicit none
  
  real*8,parameter :: pi2 = 2.0d0*acos(-1.0d0)

end module


subroutine boxmuller(u1,u2)
  
  use const
  implicit none
  
  ! Output
  real*8 u1,u2
  
  ! Intern
  real*8 z1,z2,fac,arg
  
  call random_number(z1)
  call random_number(z2)
  
  fac = sqrt(-2.0d0*log(z1)) 
  arg = pi2*z2
  u1  = fac*cos(arg)
  u2  = fac*sin(arg)
  
end subroutine



program histogram

implicit none

! Anzahl der Zufallszahlen
  integer, Parameter :: n = 1000000
   character*80 :: outfile = 'hist.dat'

  integer i,ios,k
  real*8 dum,mini,maxi,h
  real*8 :: z(n),u1,u2

! nbin ist Anzahl der Bins, axis ist die x-Achse des Histogramms
  integer,parameter :: nbin = 100
  real*8,dimension(nbin) :: axis,hist

  do i=1,n,2
!! Box Muller method
    call boxmuller(u1,u2)
    z(i) = u1
    z(i+1) = u2
  end do  



  print*, 'Find minimum and maximum ...'
  mini = z(1)
  maxi = z(1)
  do i=2,n
    if (z(i)<mini) mini = z(i)
    if (z(i)>maxi) maxi = z(i)
  end do
  
  dum = 0.001d0*(maxi-mini)
  mini = mini - dum
  maxi = maxi + dum
  
  print*, 'Minimum:', mini
  print*, 'Maximum:', maxi
  print*, ''
  
! Bestimmung der Breite eines bins
  h = (maxi - mini) / real(nbin,8)
  
! Berechnung der Histogrammachse    
  do i=1,nbin
    axis(i) = mini + (real(i-1,8) + 0.5d0)*h 
  end do
  
  
!Berechnung des Histogramms  
  do i=1,n
    k = int( (z(i)-mini)/h ) + 1
    hist(k) = hist(k) + 1
  end do
  
!Normierung  
  hist = hist/(h*real(n,8))
  print*, 'Done.'
  print*, ''
  
  print*, 'Writing histogram to file' , outfile, '...'
  open(12,file=outfile,status='replace')
  do i=1,nbin
    write(12,*) axis(i), hist(i)
  end do
  print*, 'Done.'
  
end program