頻度分布用サブルーチンのサンプル

program main

! just test for mkhist

implicit none
integer, parameter :: hist_nmax=100
integer, parameter :: data_n   = 10
integer            :: hist(hist_nmax), hist_n
real               :: min, max, dx,  dat(data_n)
integer            :: i, c, r, m


!  call system_clock( count=c, count_rate=r, count_max=m)
!  write(*,*) c, r, m
!  call random_seed(c)
  call random_number(dat)

  write(*,*) 'DATA'
  do i=1, data_n
    write(*,*) i, dat(i)
  end do

  min = 0.1; max = 0.9; dx = 0.1
  call mkhist( dat, data_n, hist, hist_nmax, hist_n, min, max, dx )

  write(*,*) 'histgram'
 i=1
  write(*,*) i, '---------------' , min+dx*(i-1), hist(i)
 do i=2, hist_n-1
  write(*,*) i, min+dx*(i-2), min+dx*(i-1), hist(i)
 end do
 i=hist_n
  write(*,*) i, min+dx*(i-2), '---------------', hist(i)


end

subroutine mkhist( x, n, hist, hist_nmax, hist_n, min, max, dx )
!  hist(  i )  = the number of data x(n) 
!  in the range : min+dx*(i-2) <= x < min+dx*(i-1)
!  hist(1) for the range x < min; hist(hist_n) for min+dx*(hist_n-2) <= x

implicit none

integer :: n, hist_nmax, hist_n, hist(hist_nmax)
real    :: x(n), min, max, dx

integer :: i, k

  hist_n = (max-min)/dx + 2

 hist = 0

 do i=1, n
  k = (x(i)-min)/dx + 2
  if( k < 1 ) then
   k = 1
  else if( k > hist_n) then
   k = hist_n
  end if

  hist(k) = hist(k) + 1
 end do
end