c-------------------------------------------------------------------- 
c                                   Satellites in N-body simulations
c                                   Make dV-dR diagram
c-------------------------------------------------------------------- 
      include 'subhalos.h' 
      INTEGER   listfiles(4)
      Character Name*80
      DATA dR/2.d0/            ! Mpc/h: width of buffer for volume replication
      DATA listfiles/1,30,31,32/
c      Open(1,file='hlist_0.9996.20.dat')
c      Open(1,file='CatshortA.DAT')
      write (*,*)  ' Enter 3 digits for snapshot '
      read (*,*) moment
      write(Name,'(a,i3.3,a)')'Subhalos.',moment,'.dat'
      Open(1,file=Name)
      write(Name,'(a,i3.3,a)')'Catshort.',moment,'.DAT'
      Open(10,file=Name)
 
      CALL ReadParticles(dR)
      write (*,*)  '<====== Number of halos read=',Nh
      write (*,30)  (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=1,1)
      write (*,30)  (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=Nh,Nh)
 30   format(' halo=',i8,' x=',3g11.3,' M=',g11.3,' Vc=',g11.3)

      Do i=1,1
      write (listfiles(i),'(T10,"Hubble=",f8.2,T28,"Box=",f8.2,"Mpc/h",
     &             /T10,"Selection: Mass limit=",g11.3,"Msun/h",
     &              T55,"Vc limit=",f8.2,"km/s",
     &              T77,"width of buffer=",f8.2,"Mpc/h")')
     &               hubble,Box, aMasslimit,Vclimit,
     &               dR
      EndDo 


      CALL ListPrimaries
      CALL Subhalos(1)


      stop
      end
c-------------------------------------------------------------------- 
c                      read particlesNx,Ny,Nz and
c                      replicate in buffer of width dR Mpc
      SUBROUTINE ReadParticles(dR)
c-------------------------------------------------------------------- 
      include 'subhalos.h' 
      Logical inx,iny,inz

      Lines =32
      Do i=1,Lines
         read(10,'(A)')Line
         If(i.le.3)write(1,*)Line
      EndDo 

      Nh = 1
 10   read(10,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh),
     &                          vxc(Nh),vyc(Nh),vzc(Nh),
     &                          aMass(Nh),rhalo(Nh),vrms(Nh),vc(Nh)
      if(Nh.ge.Nmax)goto100
      Nh =Nh +1
      goto 10

 100  Nh =Nh -1
      NpBox = Nh
      write (1,'("  Number of halos read=",T30,i8)')  Nh
      write (*,*)  '<====== Number of particles read=',Nh
      write (*,30)  (i,xc(i),yc(i),zc(i),i=1,2)
      write (*,30)  (i,xc(i),yc(i),zc(i),i=Nh-3,Nh)
 30   format(' halo=',i8,' x=',3g11.3)

c                              now relpicate periodically

      Xleft =-dR             ! new boundaries for the catalog
      Xright =Box+dR
      Yleft =-dR
      Yright =Box+dR
      Zleft =-dR
      Zright =Box+dR
      write (*,*)  Xleft ,Xright
      Ne = Nh

      Do ip =1,Nh    ! go over all halos
         x =xc(ip)
         y =yc(ip)
         z =zc(ip)

      Do k=-1,1       ! go over all 27 periodical images
         dz = k*Box
         inz =(z+dz.gt.Zleft.and.z+dz.lt.Zright)
         Do j=-1,1
           dy = j*Box
           iny =(y+dy.gt.Yleft.and.y+dy.lt.Yright)
           Do i=-1,1
               dx =i*Box
               If(i**2+j**2+k**2.ne.0)Then
                  inx =(x+dx.gt.Xleft.and.x+dx.lt.Xright)
                  If(inx.and.iny.and.inz)Then  ! add this halo to the list
                     Ne =Ne +1
                     if(Ne.gt.Nmax) STOP 'too many particles'
                     xc(Ne) =x+dx
                     yc(Ne) =y+dy
                     zc(Ne) =z+dz
                     vxc(Ne) =vxc(ip)
                     vyc(Ne) =vyc(ip)
                     vzc(Ne) =vzc(ip)
                     vc(Ne)  =vc(ip)
                     vrms(Ne)=vrms(ip)
                     aMass(Ne)=aMass(ip)
                  EndIf               ! end inside box
               EndIf                  ! end exclude central box

            EndDo       ! end i
          EndDo       ! end j
        EndDo       ! end k
      EndDo         ! end ip
      write (*,*)' Nhalos in extended sample=',Ne  
      write (1,'("  Nhals in extended sample=",T30,i8)')Ne  
      Next =Ne     ! Np is the number of particles in extened set

      close (10)
      RETURN
      END
c-------------------------------------------------------------------- 
c                  find 'clusters': halos with large Vcirc
c
      SUBROUTINE ListPrimaries
c-------------------------------------------------------------------- 
      include 'subhalos.h' 

      Nprimaries  = 0
      Do ip =1,Next               ! look for primaries in extended halo list
         If(vc(ip).gt.Vcl)Then
            Nprimaries = Nprimaries + 1
            ListCl(Nprimaries) = ip
       EndIf 
      EndDo 
      write (*,'(" N primaries=",i6," with Vc >=",2f8.3)')Nprimaries,Vcl

      RETURN 
      End

c-------------------------------------------------------------------- 
c
      SUBROUTINE Subhalos(iHalo)
c-------------------------------------------------------------------- 
      include 'subhalos.h' 
      Logical inx,iny,inz
      DIMENSION Nprof(0:Nshells),rbin(0:Nshells),
     &          Nvd(0:Nshells)

      dlogR = 0.1    ! in units of Mpc/h
      dVc = 2.5      ! in km/s
      If(iHalo.gt.Nprimaries)Then
        write(*,'(a,2i5)')' Too large index for halo:',iHalo,Nprimaries
        return
      EndIf 

      Nx = 0
      N2 = 0
      N30= 0
      N40= 0
      xx = xc(ListCl(iHalo))
      yy = yc(ListCl(iHalo))
      zz = zc(ListCl(iHalo))
      rvir = rhalo(ListCl(iHalo))
      vcprime =vc(ListCl(iHalo))
      write (1,'(T7,a,T30,a,T55,a,T72,a,T89,a,T101,a,T113,a)')
     &      'Number','Coord(Mpc/h)','Vcirc(kms)','Vrms3D','M(sun/h)',
     &      'Rvir(kpc/h)'
      write (1,'(a,i5,7g14.4)')' Halo=',ListCl(iHalo),xx,yy,zz,
     &            vcprime,vrms(ListCl(iHalo)),
     &            aMass(ListCl(iHalo)),rvir
      Do i=0,Nshells
         Nprof(i) =0
         Nvd(i)   =0
         rbin(i)  =0.
      EndDo 
      Do ip =1,Next              
        x =xc(ip)
        y =yc(ip)
        z =zc(ip)
        vv =vc(ip)
        aM =aMass(ip)
        If(vv.gt.Vclimit.and.am.gt.aMasslimit)Then
           dd   = sqrt((x-xx)**2+(y-yy)**2+(z-zz)**2)*1000.
           indx = MIN(INT(LOG10(dd)/dlogR),Nshells)
           If(dd.le.rvir)Nx  = Nx +1
           Nprof(indx) = Nprof(indx) +1
           rbin(indx)  = rbin(indx)  +dd
c           If(dd.le.rvir.and.vv.gt.40.)
c     &    write (*,'(3f8.4,3g12.4)')x,y,z,vv,aM,dd
           If(dd.lt.30.)Then
              N30 =N30+1
              write (*,'(3f8.4,3g12.4)')x,y,z,vv,aM,dd
           EndIf 
           If(dd.lt.40.)N40 =N40+1
           If(dd.le.2.*rvir)Then
              N2 =N2+1
              iv =MIN(INT(vv/dVc),Nshells)
              Nvd(iv) = Nvd(iv)+1
           EndIf 
        EndIf 
      EndDo      ! ip loop
      v0 = 0.
      write (1,'(a,i4,a,i4)')' Halos inside Rvir=',Nx,
     &                       '   Halos inside 2Rvir=',N2
      write (1,'(" R(kpc/h)   R/Rvir     dN/dV     N")')
      write (*,*)  ' Inside 30kpc n=',N30
      write (*,*)  ' Inside 40kpc n=',N40
      Do i=1,Nshells
         rout = 10.**(dlogR*i)
         rin  = 10.**(dlogR*(i-1))
         v1   = 4.1888*rout**3
         v0   = 4.1888*rin**3
         dv   = v1-v0
         Nbin = Nprof(i-1)
         dd   = Nbin/dv
         rad= rbin(i-1)/max(Nbin,1)
         if(Nbin.le.3)rad=(rout+rin)/2.
         if(rout.lt.1.e3.and.Nbin.gt.0.and.rin.gt.20.)
     &    write (1,'(2f8.3,g12.3,i5,3g12.3)')
     &        rad,rad/rvir,
     &        dd,Nprof(i-1),rin,rout
      EndDo 


      write (1,'(" V_c(km/s)      N(>Vc)  N(Vc)   dV=",f8.3)')dVc
      Nx  =0
      Do i=Nshells-1,0,-1
         Vout = dVc*(i+1)
         Vin  = dVc*i
         If(Vin.lt. 150.)Then
         Nx = Nx + Nvd(i)
         write (1,'(f8.3,2i5)')
     &        Vin,Nx,Nvd(i)
         EndIf 
      EndDo 



      RETURN 
      End

