c-------------------------------------------------------------------- 
c                        Distribution functions
c                        in halo catalogs
c-------------------------------------------------------------------- 
      PARAMETER (Nmax = 300000) ! 
      PARAMETER (Nbinmax = 70) !
      PARAMETER (Nqual   = 13) ! 
      DIMENSION Nin(Nbinmax),NinC(Nbinmax),Vc(Nbinmax),
     &           Conc(Nbinmax), hals(Nqual)
      Character  File1*100,File2*100
      DATA      Nin,NinC,Vc,Conc/Nbinmax*0,Nbinmax*0,Nbinmax*0.,
     &                          Nbinmax*0./
      Character Header*100

      write (*,'(" Enter Name of halo Catalog=",$)')
      read(*,'(A)') File1
      File2 = 'distr'//File1
      Open(1,file=File1)
      Open(2,file=File2)

      Lines = 20
      dlog  = 0.05
      Do i=1,Lines
         read(1,'(a)') Header
         write (*,'(a)')  Header
         if(i.le.3)write (2,'(a)')  Header
      EndDo


      Nh =0
 100  Nh =Nh+1
 105  read(1,*,err=400,end=400) (hals(i),i=1,Nqual)
!      if(hals(13).lt.2.5)Goto105    ! skip halos wih too small C
c         write(*,'(" Halo=",i8," coords=",3f8.3," bins=",i3)')
c     &              Nh,(hals(i,Nh),i=1,3),Nhbin(Nh)
         ind =min(max(INT(log10(hals(10))/dlog),1),Nbinmax)
c         write(*,*) ' halo=',Nh,hals(10),ind
         Nin(ind) = Nin(ind) +1
         Vc(ind)  = Vc(ind)  +hals(10)
         If(hals(13).gt.3.0)Then   ! drop halos with too small C
            Conc(ind)= Conc(ind)+hals(13) 
            NinC(ind)= NinC(ind)+1
         EndIf 
      If(Nh.lt.Nmax)goto100
      Nh =Nh +1
 400  Nh =Nh -1

      write (*,*)  ' Number of halos read =',Nh
      write (2,*)  ' Number of halos =',Nh
      write(2,'(" Bin",6x,"Nhalos",5x,"V",8x,"Vhals",
     &                    5x,"dN/dV",8x,"<C>")')  
      Do i=1,Nbinmax
        V0  =10.**(i*dlog)      
        V1  =10.**((i+1.)*dlog) 
        If(V0.gt.30..and.V0.lt.1500.)Then
           dd  =Nin(i)/(V1-V0)
           If(Nin(i).gt.0)Then
              Vcav = Vc(i)/Nin(i)
           Else
              Vcav = (V1+V0)/2.
           EndIf
           write(2,'(2i6,10g11.4)') i,Nin(i),(V1+V0)/2.,
     &             Vcav,dd,Conc(i)/max(NinC(i),1)
           write(*,'(2i6,10g11.4)') i,Nin(i),(V1+V0)/2.,
     &             Vcav,dd,Conc(i)/max(NinC(i),1)
        EndIf
      EndDo 
      stop
      end