!---------------------------------------------------
!
!               Get statistics for a halo catalog
!
!---------------------------------------------------
Program Stat
    Character :: Line*80 
    Integer, Parameter :: Nbin = 500
    Real*4,  Parameter :: dlog = 0.1

    Integer ::  ind(-Nbin:Nbin)  = 0
    Real*4  :: Mind(-Nbin:Nbin) = 0.


Open(1,file='Catshort.FullBox_z0.0_0512.0')
!Open(2,file='List.dat')

Do i=1,24
   read(1,'(a)') Line
   write(*,'(a)')Line
EndDo

Ncenter =0
aMax    =0.
aMin    =1.e+20 
iMax    =-Nbin-1
iMin    = Nbin+1

Do 
   Read(1,*,iostat=iStat) i,x,y,z,vx,vy,vz,aM,Rvir,Vrms,Vcirc
   If(iStat/= 0)Exit
      Ncenter = Ncenter +1
      iBin =  INT(log10(aM)/dlog+10000)-10000
      iBin =  max(min(iBin,Nbin),-Nbin)
      ind(iBin)  = ind(iBin)  + 1
      Mind(iBin) = Mind(iBin) + aM
      aMax  = max(aMax,aM)
      aMin  = min(aMin,aM)
      iMax  = max(iMax,iBin)
      iMin  = min(iMin,iBin)
End Do

print *,' Read Nhalos  =',Ncenter
print *,' Mass min/max =',aMin,aMax
print '(5x,a)','  N    >N       <Mass>      dN/dV      Mass      n(>Mass) '
Box = 500.
iFlag  = 0
Do i = iMin,iMax
   nn  = max(ind(i),1)
   aMass  = Mind(i)/nn
   aMbin  = 10.**(i*dlog)
   iUp  =sum(ind(i:iMax))
   write(*,'(2i8,1P,6g12.4)') ind(i),iUp,aMass,ind(i)/Box**3,aMbin,iUp/Box**3
EndDo


end Program Stat