program histogram
      implicit none

! We use nref bins in the radial range of 0...rref
      real,    parameter :: rref=1.,pi=4.*atan(1.)
      integer, parameter :: nbin=100,nref=50
      integer, dimension(nbin) :: h
      real,    dimension(nbin) :: x
      real    :: r,ri,ro,v,s
      integer :: i,j,npart

      open(1,file='grid.txt',status='old')
    1 read(1,*,end=2) i,r
        j = max(ceiling(r/rref*nref),1)
        if(j.le.nbin) then
          h(j) = h(j)+1
          x(j) = x(j)+r
        endif
      goto 1
    2 close(1)
      npart = i

      ri = 0.
      do j=1,nbin
        if(h(j).ne.0) x(j) = x(j)/h(j) ! mean radius of particles in shell
        ro = j*rref/nref
        v  = 4.*pi/3.*(ro**3-ri**3) ! volume of shell
        s  = 1./(v*npart)
        write(*,'(2i6,5es12.4)') j,h(j),x(j),ri,ro,h(j)*s,sqrt(real(h(j)))*s
        ri = ro
      enddo

      end