c-------------------------------------------------------------------- 
c                                   Satellites in N-body simulations
c                                   Make dV-dR diagram
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      INTEGER   listfiles(4)
      DATA dR/10.d0/            ! Mpc/h: width of buffer for volume replication
      DATA listfiles/2,30,31,32/
c      Open(1,file='hlist_0.9996.20.dat')
c      Open(1,file='CatshortA.DAT')
      Open(1,file='CatFit.dat')
      Open(2,file='resultsAllClout.500.dat')
      Open(30,file='isoV180dx.dat')
      Open(31,file='isoV180dy.dat')
      Open(32,file='isoV180dz.dat')
      Open(4,file='sample2.dat')
      Open(10,file='part1.000.dat',form='unformatted')
c      Open(10,file='partL.9996.dat',form='unformatted')
c      Open(10,file='part10.dat')
 

      Do i=1,4
      write (listfiles(i),'(T10,"Hubble=",f8.2,T25,"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",
     &             /T10,"Isolation: Radii=",f8.1,"kpc/h ",f8.2,"Mpc/h",
     &             T55,"los Velocity= ",f8.3,"km/s",
     &            /T10,"Clusters: Vcirc >",f8.1,"km/s  Remove halos if",
     &             " distance <=",f8.2,"Mpc/h")')
     &               hubble,Box, aMasslimit,Vclimit,
     &               dR,Radlimit,RadLmax,Vlimit,Vcl,Dcl
      EndDo 
      Lines =22
      Do i=1,Lines
         read(1,'(A)')Line
         Do j=1,4
            if(i.le.4)write(listfiles(j),*)Line
         ENDDO 
      EndDo 

      Nh = 1
 10   read(1,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh),
     &                    vxc(Nh),vyc(Nh),vzc(Nh),aMass(Nh),
     &                    vvv,vc(Nh)
      if(vc(Nh).lt.Vclimit)goto 10  ! do not take halo if too small
      if(Nh.ge.Nmax)goto100
      Nh =Nh +1
      goto 10

 100  Nh =Nh -1
      write (2,'("  Number of halos read=",T30,i8)')Nh
      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)

      CALL ReadParticles(dR)
      CALL ExpandSample(dR)
      CALL ListClusters
      Vc1 =450.
      Vc2 =550.
      CALL Isolated3D(Vc1,Vc2,Nx,Ny,Nz)
c      CALL Satellites
c      CALL OutSatellites(Nx,Ny,Nz)
      CALL SatellitesP

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

      Lines =1
      Do i=1,Lines
c         read(10,'(A)')Line
         read(10)AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW
c     +                  EKIN,
c     +                  NROWC,NGRID,NRECL,Om0,Oml0,hub,
c     +                  Ocurv
         write(*,*)' a=',AEXPN,' dA=',ASTEP,' i=',ISTEP
c         if(i.le.4)write(2,*)' a=',AEXPN,' dA=',ASTEP,' i=',ISTEP
c         if(i.le.4)write(3,*)' a=',AEXPN,' dA=',ASTEP,' i=',ISTEP
      EndDo 

      Np = 1
c 10   read(10,*,end=100,err=100)xp(Np),yp(Np),zp(Np),
  10   read(10,end=100,err=100)xp(Np),yp(Np),zp(Np),
     &                    vxp(Np),vyp(Np),vzp(Np)
      if(Np.ge.Nmaxp)goto100
      Np =Np +1
      goto 10

 100  Np =Np -1
      NpBox = Np
      write (2,'("  Number of particles read=",T30,i8)')  Np
      write (*,*)  '<====== Number of particles read=',Np
      write (*,30)  (i,xp(i),yp(i),zp(i),i=1,2)
      write (*,30)  (i,xp(i),yp(i),zp(i),i=Np-3,Np)
 30   format(' particle=',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 = Np

      Do ip =1,Np    ! go over all particles
         x =xp(ip)
         y =yp(ip)
         z =zp(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.Nmaxp) STOP 'too many particles'
                     xp(Ne) =x+dx
                     yp(Ne) =y+dy
                     zp(Ne) =z+dz
                     vxp(Ne) =vxp(ip)
                     vyp(Ne) =vyp(ip)
                     vzp(Ne) =vzp(ip)
                  EndIf               ! end inside box
               EndIf                  ! end exclude central box

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

      RETURN
      END
c-------------------------------------------------------------------- 
c                        make sample larger by replicating
c                        some halos using periodical
c                        boundary conditions
c          dR is the width of the buffer in Mpc/h
c          Next = new number of halos
      SUBROUTINE ExpandSample(dR)
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      Xleft =-dR             ! new boundaries for the catalog
      Xright =Box+dR
      Yleft =-dR
      Yright =Box+dR
      Zleft =-dR
      Zright =Box+dR
      write (*,*)  Xleft ,Xright
      Next = 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
                     Next =Next +1
                     if(Next.gt.Nmax) STOP 'too many halos'
                     xc(Next) =x+dx
                     yc(Next) =y+dy
                     zc(Next) =z+dz
                     vxc(Next) =vxc(ip)
                     vyc(Next) =vyc(ip)
                     vzc(Next) =vzc(ip)
                     aMass(Next) =aMass(ip)
                     vc(Next)        =vc(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=',Next  
      write (2,'("  Nhalos in extended sample=",T30,i8)')Next  
      RETURN 
      END
c-------------------------------------------------------------------- 
c                  find 'clusters': halos with large Vcirc
c
      SUBROUTINE ListClusters
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

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

      RETURN 
      End

c-------------------------------------------------------------------- 
c                  select isolated halos
c                  range of circular velocities Vc1<Vc< Vc2
c                  do it in 3D, not in projection
c                  Lox,y,z = 1 if halo is isolated
c                  Nx,Ny,Nz = number of isolated halos in each direction
      SUBROUTINE Isolated3D(Vc1,Vc2,Nx,Ny,Nz)
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      Nx = 0
      Ny = 0
      Nz = 0
      R2     = (RadLmax)**2  ! square of max distance in Mpc
 
      Do ip =1,Nh            ! Isolation away from clusters
            vv =vc(ip)
         If(vv.gt.Vc1.and.vv.lt.Vc2)Then
            If(DistanceCluster(ip).gt.Dcl)Then
              Lox(ip) =1
              Loy(ip) =1
              Loz(ip) =1
              Nx = Nx +1
           EndIf 
        Else
          Lox(ip) =0
          Loy(ip) =0
          Loz(ip) =0
       EndIf 
      EndDo 
      write(*,'(" NHalos=",i6," in Vc range=",2f8.3,"km/s")')Nx,Vc1,Vc2
      Nx =0
      Do ip =1,Nh              ! Test against other halos
         If(Lox(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            Do jp =1,Next
               if(jp.ne.ip)Then
                                        ! distance
                  dd = (xc(jp)-x)**2+(yc(jp)-y)**2+(zc(jp)-z)**2
                 If(dd.lt.R2.and.vc(jp).gt.vv/Vratio)Lox(ip) =0
              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Lox(ip) =1
      EndDo      ! ip loop

      Do ip =1,Nh
         If(Lox(ip).eq.1)Nx =Nx+1
      EndDo
      write(*,'(" Nisolated=",3i7)') Nx
      write(2,'("  Nisolated halos =",T30,i4," Vc limits=",3f8.2)') 
     &             Nx,Vc1,Vc2

      RETURN 
      End
c-------------------------------------------------------------------- 
c                  select isolated halos
c                  range of circular velocities Vc1<Vc< Vc2
c                  do it in 3 projections: x y z
c                  Lox,y,z = 1 if halo is isolated
c                  Nx,Ny,Nz = number of isolated halos in each direction
      SUBROUTINE Isolated(Vc1,Vc2,Nx,Ny,Nz)
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      Nx = 0
      Ny = 0
      Nz = 0
      R2     = (RadLmax)**2  ! square of max distance in Mpc
 
      Do ip =1,Nh
            vv =vc(ip)
         If(vv.gt.Vc1.and.vv.lt.Vc2)Then
            If(DistanceCluster(ip).gt.Dcl)Then
              Lox(ip) =1
              Loy(ip) =1
              Loz(ip) =1
              Nx = Nx +1
           EndIf 
        Else
          Lox(ip) =0
          Loy(ip) =0
          Loz(ip) =0
       EndIf 
      EndDo 
      write(*,'(" NHalos=",i6," in Vc range=",2f8.3,"km/s")')Nx,Vc1,Vc2
      Nx =0
      Do ip =1,Nh
         If(Lox(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            vx =vxc(ip)
            vy =vyc(ip)
            vz =vzc(ip)
            vlx =hubble*x+vx
            vly =hubble*y+vy
            vlz =hubble*z+vz
c            if(ip.eq.69)Then
c               write(*,*)' It is it: vv=',vv
c               write(*,*)' x=',x,y,z
c            EndIf 
            Do jp =1,Next
               if(jp.ne.ip)Then
                                        ! x - direction
                  dv =abs(vlx-hubble*xc(jp)-vxc(jp))
                 If(dv.lt.Vlimit)Then
                     dx2 = (yc(jp)-y)**2+(zc(jp)-z)**2
                     
                     if(dx2.lt.R2)Then
                       if(vc(jp).gt.vv/Vratio)Lox(ip) =0
                    EndIf            
c            if(ip.eq.69.and.jp.eq.10415)Then
c               write(*,*)' : vc=',vc(jp),sqrt(dx2),sqrt(R2)
c               write(*,*)' Lox(ip)=',Lox(ip)
c            EndIf 
                 EndIf           
                                        ! y - direction
                  dv =abs(vly-hubble*yc(jp)-vyc(jp))
                 If(dv.lt.Vlimit)Then
                     dy2 = (xc(jp)-x)**2+(zc(jp)-z)**2
                     if(dy2.lt.R2)Then
                       if(vc(jp).gt.vv/Vratio)Loy(ip) =0
                    EndIf        
                 EndIf           
                                        ! z - direction
                  dv =abs(vlz-hubble*zc(jp)-vzc(jp))
                 If(dv.lt.Vlimit)Then
                     dz2 = (xc(jp)-x)**2+(yc(jp)-y)**2
                     if(dz2.lt.R2)Then
                       if(vc(jp).gt.vv/Vratio)Loz(ip) =0
                    EndIf        
                 EndIf           

              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Lox(ip) =1
      EndDo      ! ip loop

      Do ip =1,Nh
         If(Lox(ip).eq.1)Nx =Nx+1
         If(Loy(ip).eq.1)Ny =Ny+1
         If(Loz(ip).eq.1)Nz =Nz+1
      EndDo
      write(*,'(" Nisolated=",3i7)') Nx,Ny,Nz
      write(2,'("  Nisolated halos =",T30,3i4," Vc limits=",3f8.2)') 
     &             Nx,Ny,Nz,Vc1,Vc2
      write(30,'("  N Halos=",i6," in Vc range=",2f8.3," km/s")')
     &    Nx,Vc1,Vc2
      write(31,'("  N Halos=",i6," in Vc range=",2f8.3," km/s")')
     &      Ny,Vc1,Vc2
      write(32,'("  N Halos=",i6," in Vc range=",2f8.3," km/s")')
     &    Nz,Vc1,Vc2
 
c      Do ip =1,Next
c         write(4,20) xc(ip),yc(ip),zc(ip),vc(ip),Lox(ip),Loy(ip),Loz(ip)
c 20      format(4g12.4,3i3)
c      EndDo

      RETURN 
      End


c-------------------------------------------------------------------- 
c                  Find satellites of isolated halos
c                  do it in 3 directions: x y z
c                  Lox,y,z = 1 if halos is isolated
      SUBROUTINE Satellites
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      Nx = 0
      Ny = 0
      Nz = 0
      R2     = (Radlimit/1000.)**2  ! square of max distance in Mpc
      write(3,*)' r_p(kpc) dv(km/s) Vcprime Vcsat  iprime   isat'
      Do ip =1,Nh                         ! x - direction
         If(Lox(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            vx =vxc(ip)
            vlx =hubble*x+vx
            Do jp =1,Next
               if(jp.ne.ip)Then
                  dv =abs(vlx-hubble*xc(jp)-vxc(jp))
                  dx = abs(xc(jp)-x)
                 If(dv.lt.Vlimit.and.dx.lt.Box/2.)Then
                     dx2 = (yc(jp)-y)**2+(zc(jp)-z)**2
                     if(dx2.lt.R2.and.vc(jp).lt.vv/Vratio)Then
                        dd =sqrt(dx2 +(xc(jp)-x)**2)
c                        uv =(vx -vxc(jp))
                        uv =sqrt((vxc(ip)-vxc(jp))**2+
     &                  (vyc(ip) -vyc(jp))**2+(vzc(ip) -vzc(jp))**2)
 
c                    write(3,'(f8.2,f8.2,2f8.2,2f9.2,2i8,7f8.3)')
c     &                 1000.*sqrt(dx2)/(hubble/100.),dv,vv,vc(jp),
c     &                  1000.*dd,uv,ip,jp,
c     &                    x,y,z,xc(jp),yc(jp),zc(jp)
                     Nx = Nx +1
                    EndIf            
                 EndIf           
              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Lox(ip) =1
      EndDo      ! ip loop

      write (*,*)  ' Satellites in x-dir =',Nx
      write (2,*)  ' Satellites in x-dir =',Nx

      Do ip =1,Nh                         ! y - direction
         If(Loy(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            vy =vyc(ip)
            vly =hubble*y+vy
            Do jp =1,Next
               if(jp.ne.ip)Then
                  dv =abs(vly-hubble*yc(jp)-vyc(jp))
                  dy = abs(yc(jp)-y)
                 If(dv.lt.Vlimit.and.dy.lt.Box/2.)Then
                     dy2 = (xc(jp)-x)**2+(zc(jp)-z)**2
                     if(dy2.lt.R2.and.vc(jp).lt.vv/Vratio)Then
                        dd =sqrt(dy2 +(yc(jp)-y)**2)
                        uv =sqrt((vxc(ip)-vxc(jp))**2+
     &                  (vyc(ip) -vyc(jp))**2+(vzc(ip) -vzc(jp))**2)
                     write(3,'(f8.2,f8.2,4f9.2,2i8,7f8.3))')
     &                 1000.*sqrt(dy2)/(hubble/100.),dv,vv,vc(jp),
     &                   1000.*dd,uv,ip,jp,
     &                    x,y,z,xc(jp),yc(jp),zc(jp)
                     Ny = Ny +1
                    EndIf            
                 EndIf           
              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Loy(ip) =1
      EndDo      ! ip loop

      write (*,*)  ' Satellites in y-dir =',Ny
      write (2,*)  ' Satellites in y-dir =',Ny

      Do ip =1,Nh                         ! z - direction
         If(Loz(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            vz =vzc(ip)
            vlz =hubble*z+vz
            Do jp =1,Next
               if(jp.ne.ip)Then
                  dv =abs(vlz-hubble*zc(jp)-vzc(jp))
                  dz = abs(zc(jp)-z)
                 If(dv.lt.Vlimit.and.dz.lt.Box/2.)Then
                     dz2 = (xc(jp)-x)**2+(yc(jp)-y)**2
                     if(dz2.lt.R2.and.vc(jp).lt.vv/Vratio)Then
                        dd =sqrt(dz2 +(zc(jp)-z)**2)                        
c                      uv =(vz -vzc(jp))
                        uv =sqrt((vxc(ip)-vxc(jp))**2+
     &                  (vyc(ip) -vyc(jp))**2+(vzc(ip) -vzc(jp))**2)

c                     write(3,'(f8.2,f8.2,4f9.2,2i8,7f8.3)')
c     &                 1000.*sqrt(dz2)/(hubble/100.),dv,vv,vc(jp),
c     &                   1000.*dd,uv,ip,jp,
c     &                    x,y,z,xc(jp),yc(jp),zc(jp)
                     Nz = Nz +1
                    EndIf            
                 EndIf           
              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Loz(ip) =1
      EndDo      ! ip loop

      write (*,*)  ' Satellites in z-dir =',Nz
      write (2,*)  ' Satellites in z-dir =',Nz

      RETURN 
      End
c-------------------------------------------------------------------- 
c                  Find particles of isolated halos and
c                  Place them in a file
c                  Lox,y,z = 1 if halos is isolated
      SUBROUTINE OutSatellites(Nx,Ny,Nz)
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      Radmx = 15. ! in Mpc
      Nselect   = 100  ! number of particles to select
      Nseed = 12077121  ! initial seed for random numbers
      dRmax    = 0.500  ! max radius in projection
      dVmax    = 500.   ! velocity limit for particles km/s
      R2     = (Radmx)**2  ! square of max distance in Mpc
      Niso   = 0
      Do i=0,2
         write(30+i,*)' Number of particles per halo  =',Nselect
         write(30+i,*)' Max 3D radius around halo     =',Radmx,'Mpc/h'
         write(30+i,*)' Max projected radius          =',dRmax,'Mpc/h'
         write(30+i,*)' Format of tables:'
         write(30+i,*)' HaloNumb Vcirc  Mass(Msun/h) Coords(Mpc/h)'
         write(30+i,'("   dv(km/s) r_p(kpc/h)",3x,"dXYZ(kpc/h)",
     &          30x,"dVxyz(km/s)")')
      EndDo 

      Do ip =1,Nh
         If(Lox(ip).eq.1.or.Loy(ip).eq.1.or.Loz(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            Niso = Niso +1
            Nnx   = 0
            Nny   = 0
            Nnz   = 0
            write (*,'(i8,f8.2,g11.4,3f8.3)') ip,vv,aM,x,y,z
            If(Lox(ip).eq.1)
     &        write (30,'(i8,f8.2,g11.4,3f8.3)') ip,vv,aM,x,y,z
            If(Loy(ip).eq.1)
     &        write (31,'(i8,f8.2,g11.4,3f8.3)') ip,vv,aM,x,y,z
            If(Loz(ip).eq.1)
     &        write (32,'(i8,f8.2,g11.4,3f8.3)') ip,vv,aM,x,y,z
            Do jp =1,Np    ! Find how many particles are in cylinder
               dd =(xp(jp)-x)**2+(yp(jp)-y)**2+(zp(jp)-z)**2
               if(dd.lt.R2)Then
                     dx = xp(jp)-x
                     dy = yp(jp)-y
                     dz = zp(jp)-z
                     dd = sqrt(dd)
                     dvx = vxp(jp)-vxc(ip) +hubble*dx
                     dvy = vyp(jp)-vyc(ip) +hubble*dy
                     dvz = vzp(jp)-vzc(ip) +hubble*dz
                     dxp = sqrt(dy**2+dz**2)
                     dyp = sqrt(dx**2+dz**2)
                     dzp = sqrt(dx**2+dy**2)
                  If(abs(dvx).lt.dVmax.and.dxp.lt.dRmax)Nnx = Nnx+1  
                  If(abs(dvy).lt.dVmax.and.dyp.lt.dRmax)Nny = Nny+1  
                  If(abs(dvz).lt.dVmax.and.dzp.lt.dRmax)Nnz = Nnz+1
               EndIf     
            EndDo 
            FractionX = 100./Nnx*1.05
            FractionY = 100./Nny*1.05
            FractionZ = 100./Nnz*1.05
            Nseed_a = Nseed
 200        Nnx   = 0
            Nny   = 0
            Nnz   = 0
            Nnxp   = 0
            Nnyp   = 0
            Nnzp   = 0
             Do jp=1,Np     ! select Nselect particles in each direction
              If(Nnx.lt.Nselect.or.Nny.lt.Nselect.or.Nnz.lt.Nselect)Then
                  dd =(xp(jp)-x)**2+(yp(jp)-y)**2+(zp(jp)-z)**2
                  if(dd.lt.R2)Then
                        dx = xp(jp)-x
                        dy = yp(jp)-y
                        dz = zp(jp)-z
                        dd = sqrt(dd)
                        dvx = vxp(jp)-vxc(ip) +hubble*dx
                        dvy = vyp(jp)-vyc(ip) +hubble*dy
                        dvz = vzp(jp)-vzc(ip) +hubble*dz
                        dxp = sqrt(dy**2+dz**2)
                        dyp = sqrt(dx**2+dz**2)
                        dzp = sqrt(dx**2+dy**2)
                     If(Lox(ip).eq.1)Then
                     If(abs(dvx).lt.dVmax.and.dxp.lt.dRmax
     &                 .and.RANDd(Nseed).lt.FractionX)Then
                     If(Nnx.lt.Nselect)Then
                      Nnx = Nnx+1  
                      If(dvx.gt.0.)Nnxp = Nnxp+1
                       write(30,'(f9.2,f9.2,2x,3f10.2,5x,9f9.2)')
     &                 dvx,1000.*dxp,
     &                 1.e3*dx,1.e3*dy,1.e3*dz,dvx,dvy,dvz
                     EndIf                        
                     EndIf 
                     EndIf    
                     If(Loy(ip).eq.1)Then
                     If(abs(dvy).lt.dVmax.and.dyp.lt.dRmax
     &                 .and.RANDd(Nseed).lt.FractionY)Then
                     If(Nny.lt.Nselect)Then
                      Nny = Nny+1  
                      If(dvy.gt.0.)Nnyp = Nnyp+1
                       write(31,'(f9.2,f9.2,2x,3f10.2,5x,3f9.2)')
     &                 dvy,1000.*dyp,
     &                 1.e3*dx,1.e3*dy,1.e3*dz,dvx,dvy,dvz
                     EndIf 
                     EndIf  
                     EndIf
                     If(Loz(ip).eq.1)Then    
                     If(abs(dvz).lt.dVmax.and.dzp.lt.dRmax
     &                 .and.RANDd(Nseed).lt.FractionZ)Then
                     If(Nnz.lt.Nselect)Then
                      Nnz = Nnz+1  
                      If(dvz.gt.0.)Nnzp = Nnzp+1
                       write(32,'(f9.2,f9.2,2x,3f10.2,5x,3f9.2)')
     &                 dvz,1000.*dzp,
     &                 1.e3*dx,1.e3*dy,1.e3*dz,dvx,dvy,dvz
                     EndIf 
                     EndIf  
                     EndIf    
                  EndIf    
              EndIf
            EndDo   ! jp loop
            iflag = 0
c       If(Nnx.ne.Nselect)write(*,*)'  wrong number of selected Nnx=',Nnx
c       If(Nny.ne.Nselect)write(*,*)'  wrong number of selected Nny=',Nny
c       If(Nnz.ne.Nselect)write(*,*)'  wrong number of selected Nnz=',Nnz
       If(Nnx.lt.Nselect.and.Nnx.ne.0)Then
          iflag =1
          FractionX =FractionX*1.05
          IF(FractionX.ge.1.)stop '-- Not enough particles X'
       EndIf
       If(Nny.lt.Nselect.and.Nny.ne.0)Then
          iflag =1
          FractionY =FractionY*1.05
          IF(FractionY.ge.1.)stop '-- Not enough particles Y'
       EndIf
       If(Nnz.lt.Nselect.and.Nnz.ne.0)Then
          iflag =1
          FractionZ =FractionZ*1.05
          IF(FractionZ.ge.1.)stop '-- Not enough particles Z'
       EndIf
       If(iflag.eq.1)Then
          write (*,'(15x,"Change fractions=",3f8.3," Nselected=",3i4)')
     &          FractionX,FractionY,FractionZ,Nnx,Nny,Nnz
          Nseed =Nseed_a
          Do i =1,Nnx
             BACKSPACE 30
          EndDo 
          Do i =1,Nny
             BACKSPACE 31
          EndDo
          Do i =1,Nnz
             BACKSPACE 32
          EndDo
 
          goto200
       EndIf 

         write (*,'(" Number of positives=",3i4," Total=",3i4)')
     &           Nnxp,Nnyp,Nnzp,Nnx,Nny,Nnz
         EndIf         ! Lox(ip) =1

      EndDo      ! ip loop
      write (*,*)  ' Number of isolated halos=',Niso
      write (2,'("  Number of isolated halos =",T30,i9)')Niso
c     write (*,*)  ' Particles used as satellites=',Nselect
c     write (2,'("  Particles used as satellites =",T30,i9)')Nx

      RETURN 
      End
c-------------------------------------------------------------------- 
c                  Find particles of isolated halos
c
c                  Lox,y,z = 1 if halos is isolated
      SUBROUTINE SatellitesP
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz
      DIMENSION sigV(Nshells),sigVr(Nshells),sigdV(Nshells),ns(Nshells)
      DIMENSION sigVp(0:1,Nshells),sigVrp(0:1,Nshells),
     &                  sigdVp(0:1,Nshells),nsp(0:1,Nshells),Nisop(0:1)
      DATA sigVp/Nshell2*0./,sigVrp/Nshell2*0./,sigdVp/Nshell2*0./
      DATA nsp/Nshell2*0/

      dRad = 50.    ! in kpc
      Radmx = 5. ! in Mpc
      Nx   = 0
      Nisop(0) = 0       ! number of isolated halos
      Nisop(1) = 0       ! number of isolated halos
      R2     = (Radmx)**2  ! square of max distance in Mpc
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( k,jj,ip,jp,dd,dx,dy,dz,dvx,dvy,dvz)
C$OMP+PRIVATE (uv,uu,ur,ind,x,y,z,vv,aM)
      Do k=0,1    
      Do jj =1,Nh,2
         ip =jj+k
         If(Lox(ip).eq.1.and.Loy(ip).eq.1.and.Lox(ip).eq.1)Then
            Nisop(k) = Nisop(k) +1
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            vv =vc(ip)
            aM =aMass(ip)
            Do jp =1,Np
                  dd =(xp(jp)-x)**2+(yp(jp)-y)**2+(zp(jp)-z)**2
                     if(dd.lt.R2)Then
                        dx = xp(jp)-x
                        dy = yp(jp)-y
                        dz = zp(jp)-z
                        dd =sqrt(dd)
                        dvx = vxp(jp)-vxc(ip) +hubble*dx
                        dvy = vyp(jp)-vyc(ip) +hubble*dy
                        dvz = vzp(jp)-vzc(ip) +hubble*dz
                        uv =dvx**2+dvy*2+dvz**2
                        uu =vxp(jp)**2+vyp(jp)**2+vzp(jp)**2
                        ur =(dvx*dx+dvy*dy+dvz*dz)/max(dd,1.e-2)
                        ind =min(INT(dd*1000./dRad)+1,Nshells)
                        nsp(k,ind)     =nsp(k,ind) +1
                        sigVp(k,ind) =sigVp(k,ind) +uu
                        sigdVp(k,ind) =sigdVp(k,ind) +uv
                        sigVrp(k,ind) =sigVrp(k,ind) +ur
 
c                    write(3,'(f8.2,f8.2,2f8.2,2f9.2,2i8,7f8.3)')
c     &                 1000.*sqrt(dx2)/(hubble/100.),dv,vv,vc(jp),
c     &                  1000.*dd,uv,ip,jp,
c     &                    x,y,z,xc(jp),yc(jp),zc(jp)
c                     Nx = Nx +1
                    EndIf            ! dd < R2
            EndDo   ! jp loop

         EndIf         ! Lox(ip) =1
      EndDo      ! ip loop
            EndDo   ! k loop
      Niso =Nisop(0)+Nisop(1)
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( j)
      do j=1,Nshells
         sigV(j) =sigVp(0,j) + sigVp(1,j)
         ns(j)     =nsp(0,j)    + nsp(1,j)
         sigV(j)  =sigVp(0,j)+ sigVp(1,j)
         sigdV(j)  =sigdVp(0,j)+ sigdVp(1,j)
         sigVr(j)  =sigVrp(0,j)+ sigVrp(1,j)
      EndDo 

      write (*,*)  ' Number of isolated halos=',Niso
      write (2,'("  Number of isolated halos =",T30,i9)')Niso
c     write (*,*)  ' Particles used as satellites=',Nx
c     write (2,'("  Particles used as satellites =",T30,i9)')Nx
      write (*,*) ' n     R(kpc)  sigVpart     sigdV       Vr   Density'
      write (2,'(" n     R(kpc)  sigVpart     sigdV       Vr",
     &               "    OverDensity IntegratedoverDens  Nint")')
      nint =0
      Do i=1,Nshells
         Volume =4.1888*dRad**3*(i**3-max(0,i-1)**3)
         nint = nint +ns(i)        ! number of particles inside current radius 
         If(ns(i).gt.2) write(2,'(i6,6g11.4,i8,g12.4)')ns(i),
     &              (i-0.5)*dRad,sqrt(sigV(i)/ns(i)),
     &                   sqrt(sigdV(i)/ns(i)),sigVr(i)/ns(i),
     &              ns(i)/(Niso*float(NpBox)*Volume/(1000.*Box)**3),
     &              nint/(Niso*float(NpBox)*(4.1888*(dRad*i)**3)
     &                         /(1000.*Box)**3),
     &              nint

       EndDo
      RETURN 
      End
c-------------------------------------------------------------------- 
c                  find 'clusters': halos with large Vcirc
c
      Real Function DistanceCluster(i)
c-------------------------------------------------------------------- 
      include 'satellites.h' 
      Logical inx,iny,inz

      DistanceCluster =1.e10
      Do icl =1,Nclusters               ! look for clusters in extended halo list
         ip =ListCl(icl)                     ! pointer to cluster in full halo list 
         dd = (xc(i) -xc(ip))**2+
     &          (yc(i) -yc(ip))**2+(zc(i) -zc(ip))**2
         DistanceCluster =min(dd,DistanceCluster)
      EndDo 
      DistanceCluster =sqrt(DistanceCluster)

      RETURN 
      End

C------------------------------------------------
C				                                       random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END