!----------------------------------------------------------------------------------------
!                    Find properties of large distinct halos
!                   
!     Input:  Halo Catalog (short list)
!                
!     Output:  Distribution functions
!     Parameters:
!                     
!                      
!
MODULE   setHalos
PUBLIC      

Real,   PARAMETER    ::                                        &
                         dLog              = 0.075,                    & ! bin size in log(r)
                         dMlog            =0.05
Integer, ALLOCATABLE    :: LDistinct(:),Lst(:),Label(:,:,:)
Real,     ALLOCATABLE, DIMENSION (:)     :: Xc,Yc,Zc,         &
                                           VXc,VYc,Vzc,      &
                                           Amc,Rmc,Vrms,Vcirc
Integer :: Ncenter,Nmain,Nsample,Nmx,Nbx,Nmy,Nby,Nmz,Nbz
Real      :: Box, Cell , dBuffer,AEXPN                                          
Contains

!--------------------------------------------------------------
!                                Make list of  distinct halos
!    Input:            Ncenter - numer of halos, {Xc ...} = halo data     
!
!    Output:         LDistinct(i) = 1 for distinct halo
!                          
      SUBROUTINE HaloDistinct
!--------------------------------------------------------------

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
        Do i=1,Ncenter
           LDistinct(i) =1
        EndDo

        !write(*,*) 'Box=',Box,dBuffer,Ncenter
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (ic,jc,x,y,z,R0,Vcc,dd, &
!$OMP                     N0x,N1x,N0y,N1y,N0z,N1z, i,j,k) 
        Do ic=1,Nmain
           If(mod(ic,10000)==0)write(*,*)' Distinct: ',ic
           x = Xc(ic); y = Yc(ic); z = Zc(ic)
           R0    =Rmc(ic)
           Vcc =Vcirc(ic)
           N0x = max(INT((x-dBuffer)/Cell+100)-100, Nmx)  ! limits for linker list
           N1x = min(INT((x+dBuffer)/Cell+100)-100,  Nbx)  
           N0y = max(INT((y-dBuffer)/Cell+100)-100, Nmy) 
           N1y = min(INT((y+dBuffer)/Cell+100) -100, Nby) 
           N0z = max(INT((z-dBuffer)/Cell+100)-100, Nmz)  
           N1z = min(INT((z+dBuffer)/Cell+100) -100, Nbz)
           Do i =N0x,N1x
           Do j =N0y,N1y
           Do k =N0z,N1z
              jc = Label(i,j,k)
              Do while (jc.ne.0)                         ! loop over all particles in cell
               If(jc /= ic.and.Vcirc(jc)>=Vcc)Then
                  dd =(Xc(jc)-x)**2 +(Yc(jc)-y)**2 +(Zc(jc)-z)**2

                  If(dd< (Max(Rmc(jc),R0))**2) Then
                       If(Vcirc(jc).gt.Vcc.or.jc>ic)LDistinct(ic) =0        ! subhalo
                        If(Vcc.gt.0.95*Vcirc(jc))LDistinct(ic) =-1  ! Fake
                   EndIf        ! dd<Rmc
               EndIf          ! Rmc => R0
               jc = Lst(jc)
               End Do
           End Do ! k
           End Do ! j
           End Do ! i
         EndDo  ! ic

        write(*,*) 'Done LDistinct'
         nn=0
         nf =0
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i) REDUCTION(+:nn,nf)
        Do i=1,Nmain
           If(LDistinct(i).eq.1)nn=nn+1
           If(LDistinct(i).eq.-1)nf=nf+1
           !if(Vcirc(i)>1600.)write(*,'(i8,i3,f8.2,3x,3f9.3)')i,LDistinct(i),Vcirc(i),Xc(i),Yc(i),Zc(i)
        EndDo
        write(2,*) 'Number of distinct halos =',nn,' Fake=',nf
!        Do i=1,Nmain
!          Write(2,'(i10,i3,1P,2g12.4)')i,LDistinct(i),aMc(i),Vcirc(i)
!         EndDo
      Return
      End SUBROUTINE HaloDistinct
!--------------------------------------------------------------
!                      add buffer of periodically replicated halos
!                      Width of the buffer is dBuffer
!                      Number of particles in the expanded sample is Ntotal
      SUBROUTINE Expand
!--------------------------------------------------------------
        Logical :: inside,insidex,insidey,insidez
        Nsample = Nmain

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (ic,i,j,k,x,y,z,x1,y1,z1, inside,insidex,insidey,insidez) 
        Do ic=1,Nmain
               x = Xc(ic); y = Yc(ic); z = Zc(ic)
             do k=-1,1
                z1 = z+k*Box                
                If(z1.le.Box+dBuffer.and.z1.ge.-dBuffer)Then
                   insidez =.true.
                Else
                   insidez =.false.
                EndIf
             do j=-1,1  
                y1 = y+j*Box
                If(y1.le.Box+dBuffer.and.y1.ge.-dBuffer)Then
                   insidey =.true.
                Else
                   insidey =.false.
                EndIf
             do i=-1,1  
                x1 = x+i*Box
                If(x1.le.Box+dBuffer.and.x1.ge.-dBuffer)Then
                   insidex =.true.
                Else
                   insidex =.false.
                EndIf
                inside = insidex.and.insidey.and.insidez
                If(abs(i) +abs(j) +abs(k) .eq. 0)inside=.false.
                If(inside)Then
!$OMP critical
                Nsample = Nsample +1
                If(Nsample.le.Ncenter)Then
                   Xc(Nsample) = x1
                   Yc(Nsample) = y1
                   Zc(Nsample) = z1
                   VXc(Nsample) = VXc(ic)
                   VYc(Nsample) = VYc(ic)
                   VZc(Nsample) = VZc(ic)
                   Amc(Nsample) =Amc(ic)
                   Rmc(Nsample) =Rmc(ic)
                   Vcirc(Nsample) =Vcirc(ic)
                EndIf
!$OMP end critical
                Endif
        EndDo
        EndDo
        EndDo

        EndDo

        write(*,*) 'Buffer width=',dBuffer,' New number of halos=',Nsample
        write(*,*) ' Maximum number of halos=',Ncenter
        If(Nsample > Ncenter) Stop ' Error: Too many centers. Increase fExpand '
        Ntotal = Nsample
     End  SUBROUTINE Expand
!--------------------------------------------------------------
!                         Mass function of distinct halos
!  
      SUBROUTINE MassDistinct
!--------------------------------------------------------------
Integer, parameter :: Nt =1000
Integer, Dimension(-Nt:Nt) ::        &
                             mHalos,iHalos
Real, Dimension(-Nt:Nt) ::  massbin

mHalos    = 0
massbin = 0.

         nn=0             ! count all halos
        Do i=1,Nmain
                        If(LDistinct(i) == 1)Then
                           nn =nn +1
                     ind =INT(log10(Amc(i))/dMlog+1000)-1000
                     If(ind<Nt)Then
                            mHalos(ind) = mHalos(ind) +1
                            massbin(ind) = massbin(ind) +Amc(i)
                     EndIf
            EndIf
        EndDo
        massbin = massbin/(max(mHalos,1))

        istart =-Nt      ! find first and last non-empty bin
        Do i =-Nt,Nt
           If(mHalos(i)>0)Then
              istart =i ; exit
           EndIf
        EndDo
        ifinish =Nt
        Do i =Nt,-Nt,-1
           If(mHalos(i)>0)Then
              ifinish =i ; exit
           EndIf
        EndDo
        iHalos  = mHalos
        Do i =ifinish-1,istart,-1
           mHalos(i) = mHalos(i) +mHalos(i+1)
        EndDo

        write(2,*) ' Mass(Msunh)    N(>M)       Mass     MdN/dM      Nbin'
        Do i=istart,ifinish         ! print results
           If(mHalos(i)>0)Then
              aMm =10.**(i*dMlog)
              dM    = 10.**((i+1)*dMlog)-aMm
              write(2,'(3x,1P,G11.4,0P,i9,2x,1P,2g11.4,i8)')aMm, mHalos(i),massbin(i), &
                       massbin(i)*iHalos(i)/dM/Box**3,iHalos(i)
           EndIf
        EndDo
      End SUBROUTINE  MassDistinct
!--------------------------------------------------------------
!           Vcirc -Mir relation for Distinct and Subs
!  
      SUBROUTINE VcircMvir
!--------------------------------------------------------------
Integer, Parameter :: M =250000, Nbin =100
Integer, SAVE, Dimension(M) :: indR 
Real, SAVE,    Dimension(M) :: DD,SD 

   write(2,'(2a,i10)')'Ndistinct  <M200>      sigma     V10    Vmax/200median   V90', & 
                               ' Nhalos=',Nmain
                            !  '          Nsubs     <M>       V10       Vmedian      V90', & 

  Do ibin =1,Nbin
   aMin = 1.e11*10.**((ibin-1)*dlog)
   aMax =1.e11*10.**((ibin)*dlog)
   If(aMax>1.e15.or.aMin<1.e10)exit
   iC =0
   ss =0.
   iCb =0
   sb =0.

   Do i=1,Nmain
      If(aMc(i)>=aMin.and.aMc(i)<aMax)Then
         If(LDistinct(i) == 1)Then
            iC =iC +1
            If(iC>M)Stop ' Not enough space for buffer. Increase M'
            DD(iC) = Vcirc(i)
            ss = ss + aMc(i)
        Else
        if(LDistinct(i) == 0)Then
            iCb =iCb +1
            If(iCb>M)Stop ' Not enough space for buffer. Increase M'
            SD(iCb) = Vcirc(i)
            sb = sb + aMc(i)
       EndIf
      EndIf
      EndIf
   EndDo
   write(*,*) ' --> aMin,aMax=',aMin,aMax,iC
      If(iC>8)Then
         CALL R_mrgref (DD,indR,iC)
         M10 = INT(0.1*iC)+1
         M50 = INT(0.5*iC)+1
         M90 = INT(0.9*iC)+1
         V10 = DD(indR(M10))
         V50 = DD(indR(M50))
         V90 = DD(indR(M90))
      EndIf
      If(iCb>8)Then
         CALL R_mrgref (SD,indR,iCb)
         M10 = INT(0.1*iCb)+1
         M50 = INT(0.5*iCb)+1
         M90 = INT(0.9*iCb)+1
         Vs10 = SD(indR(M10))
         Vs50 = SD(indR(M50))
         Vs90 = SD(indR(M90))
         Else
            Vs10 =0.; Vs50 =0.; Vs90 =0.
      EndIf
      Om0    = 0.27
      OL0     = 1. - Om0
      x        = (OL0/Om0)**(1./3.)*AEXPN
      Om      = 1./(1.+x**3)
      OL       = 1.-Om
      Da0   = 2.5*Om0/(Om0**(4./7.)-OL0+(1.+Om0/2.)*(1.+OL0/70.))
      Da     = 2.5*AEXPN*Om/(Om**(4./7.)-OL+(1.+Om/2.)*(1.+OL/70.))/Da0
      y         = 1.e12/(ss/iC)
      sigma = Da*16.9*y**0.41/(1.+1.102*y**0.20+6.22*y**0.333)
      If(iC>10)Then
         !write(*,'(a,2i8,a,1p,6g12.5)') ' Selected halos=',iC,iCb,' Masses=',aMin,aMax,V50,Vs50
         write(2,'(2(i8,1p,5g12.5),6i7)') iC,ss/iC,sigma,V10,V50,V90!,  iCb,sb/max(iCb,1),Vs10,Vs50,Vs90
      EndIf
EndDo
 
      End SUBROUTINE VcircMvir
!!--------------------------------------------------------------
!!                  Find Milky Way candidates       
!!  
!      SUBROUTINE FindMilkyWay
!!--------------------------------------------------------------
!Integer, Dimension(-100:100) ::        &
!                             mHalos,iHalos
!Real, Dimension(-100:100) ::  massbin

!mHalos    = 0
!massbin = 0.

!         nn=0             ! count all halos
!        Do i=1,Nmain
!                        If(LDistinct(i) == 1.and.Amc(i)>1.e12.and.Amc(i)<2.e12)Then
!                           nn =nn +1
!              x = Xc(ic); y = Yc(ic); z = Zc(ic)
!              R0 =Rmc(ic)
!              aM0 =Amc(ic)
!              Vcc =Vcirc(ic)
!              N0x = (x-dBuffer)/Cell ; N1x = (x+dBuffer)/Cell  ! limits for linker list
!              N0y = (y-dBuffer)/Cell ; N1y = (y+dBuffer)/Cell 
!              N0z = (z-dBuffer)/Cell ; N1z = (z+dBuffer)/Cell 
!           Do i =N0x,N1x
!              i1 =mod(i+50*(Nbx+1),Nbx+1)  ! periodical boundaries
!           Do j =N0y,N1y
!              j1 =mod(j+50*(Nby+1),Nby+1)  ! periodical boundaries
!           Do k =N0z,N1z
!              k1 =mod(k+50*(Nbz+1),Nbz+1)  ! periodical boundaries
!              jc = Label(i1,j1,k1)
!              Do while (jc.ne.0)                         ! loop over all particles in cell
 
!               If(jc /= ic.and.LDistinct(jc)==0)Then
!                  ux = Xc(jc)-x
!                  uy = Yc(jc)-y
!                  uz = Zc(jc)-z
                                                    ! periodical boundary conditions
!                  dx = min(abs(ux),abs(ux+Box),abs(ux-Box))
!                  dy = min(abs(uy),abs(uy+Box),abs(uy-Box))
!                  dz = min(abs(uz),abs(uz+Box),abs(uz-Box))

!                  dd =sqrt(dx**2 +dy**2 +dz**2)/R0
!                  If(dd<1.) Then                               ! satellite
!                      ind =MIN(INT(log10(Vcirc(jc))/dlog+100)-100,100)
!!$OMP atomic
!                          iSat(ind) = iSat(ind) +1
!                   EndIf        ! dd< R0
!               EndIf   ! jc /= ic
!               jc = Lst(jc)
!               End Do ! jc
!           End Do ! k
!           End Do ! j
!           End Do ! i

!                     EndIf
!            EndIf
!        EndDo

!      End SUBROUTINE  FindMilkyWay
!--------------------------------------------------------------
!                             Find satellite halos
!  
      SUBROUTINE Satellites
!--------------------------------------------------------------
Integer, Dimension(-100:100) ::  iSat,iHalos,iFake    
Real, Dimension(-100:100) ::  vbin

iSat       = 0 
iHalos   = 0 
vbin       = 0. 
iFake    = 0  
                              ! count satellites
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (ic,jc,x,y,z,R0,aM0,Vcc,dx,ux,dy,uy,dz,uz,dd,ind, &
!$OMP                    N0x,N1x,N0y,N1y,N0z,N1z, i,i1,j,j1,k,k1 ) 
        Do ic=1,Nmain
           If(LDistinct(ic).eq.1)Then   ! distinct halo
              x = Xc(ic); y = Yc(ic); z = Zc(ic)
              R0 =Rmc(ic)
              aM0 =Amc(ic)
              Vcc =Vcirc(ic)
              N0x = INT((x-dBuffer)/Cell+100.)-100 ; N1x = (x+dBuffer)/Cell  ! limits for linker list
              N0y = INT((y-dBuffer)/Cell+100.)-100 ; N1y = (y+dBuffer)/Cell 
              N0z = INT((z-dBuffer)/Cell+100.)-100 ; N1z = (z+dBuffer)/Cell 
           Do i =N0x,N1x
              i1 =mod(i+50*(Nbx+1),Nbx+1)  ! periodical boundaries
           Do j =N0y,N1y
              j1 =mod(j+50*(Nby+1),Nby+1)  ! periodical boundaries
           Do k =N0z,N1z
              k1 =mod(k+50*(Nbz+1),Nbz+1)  ! periodical boundaries
              jc = Label(i1,j1,k1)
              Do while (jc.ne.0)                         ! loop over all particles in cell
 
               If(jc /= ic.and.LDistinct(jc)==0)Then
                  ux = Xc(jc)-x
                  uy = Yc(jc)-y
                  uz = Zc(jc)-z
                                                    ! periodical boundary conditions
                  dx = min(abs(ux),abs(ux+Box),abs(ux-Box))
                  dy = min(abs(uy),abs(uy+Box),abs(uy-Box))
                  dz = min(abs(uz),abs(uz+Box),abs(uz-Box))

                  dd =sqrt(dx**2 +dy**2 +dz**2)/R0
                  If(dd<1.) Then                               ! satellite
                      ind =MIN(INT(log10(Vcirc(jc))/dlog+100)-100,100)
!$OMP atomic
                          iSat(ind) = iSat(ind) +1
                   EndIf        ! dd< R0
               EndIf   ! jc /= ic
               jc = Lst(jc)
               End Do ! jc
           End Do ! k
           End Do ! j
           End Do ! i
           EndIf     ! LDistinct
         EndDo  ! ic

         nn=0             ! count all halos
        Do i=1,Nmain
                     ind =INT(log10(Vcirc(i))/dlog+100)-100
                     If(ind<100)Then
                        If(LDistinct(i)>=0)Then
                            iHalos(ind) = iHalos(ind) +1
                            vbin(ind) = vbin(ind) +Vcirc(i)
                        Else
                           iFake(ind) =iFake(ind) +1
                     EndIf
            EndIf
        EndDo
        vbin = vbin/(max(iHalos,1))

        istart =-100      ! find first and last non-empty bin
        Do i =-100,100
           If(iHalos(i)>0)Then
              istart =i ; exit
           EndIf
        EndDo
        ifinish =100
        Do i =100,-100,-1
           If(iHalos(i)>0)Then
              ifinish =i ; exit
           EndIf
        EndDo
        Do i =ifinish-1,istart,-1
           iHalos(i) = iHalos(i) +iHalos(i+1)
           iSat(i)     = iSat(i)      +iSat(i+1)
           iFake(i)  = iFake(i)   +iFake(i+1)
        EndDo

        write(*,*) ' Vcirc(km/s)   All(>V)  Sat(>V)  Fake(>V) Sat/All'
        Do i=istart,ifinish         ! print results
           If(iHalos(i)>0)Then
              v =10.**(i*dlog)
              write(*,'(3x,f9.2,3x,3i7,f8.3)')v, iHalos(i), iSat(i),iFake(i),float(iSat(i))/iHalos(i)
           EndIf
        EndDo
      Return
      End SUBROUTINE Satellites

!---------------------------------------------------------
!                         Read halos catalog
!                       
      SUBROUTINE ReadHalos(iFlag,vcLimit)
!--------------------------------------------------------------
Character  :: Line*80, Txt*5='Npart', Txt3*8='Coordina',Txt2*8,txt1*4

   Nmain  =0
   Rdmax    =0.
   Do i=1,2
      Read(3,'(a)')Line
      If(iFlag==0)write(*,'(a)') Line
      If(iFlag==1)write(2,'(a)') Line
   EndDo
!  Read(3,'(a)')Line
!  Read(3,'(a)')Line
      Read(3,'(a5,i4)')txt1,ii
      indx =0
      If(iFlag==0)write(*,'(a,i6)') txt1,ii
   Do i=1,7 !21 
     read(3,*,err=20)Txt2
     write(*,*)TRIM(Txt2)
     !If(TRIM(Txt2)==Txt3)Then
     !   If(indx==1)Exit
     !   indx =1
     !   exit
     !EndIf
   EndDo
      ! stop ' ------'
   Do
     read (3,*,iostat=iStat) ii,X0, Y0 , Z0, VvX,VvY,VvZ,   &
                       aM,aMtot,rr,vrmss,vcrc,aMbn,vrmss2,VirRat,dRRvir, AngMom,aLambda,R200,V200,aM200
                                              ! aM,rr,vrmss,vcrc,Npart,Rmax,Conc !,Nb
     !read (3,*,iostat=iStat) ii,X0, Y0 , Z0, &
     !                                          aM0,aM,vnow,vcrc
     !rr = aM0**0.33333/4.82e1      !  virial radius in kpch for Om=0.27 delta=360.8

     !write (*,'(i9,3f10.4,2g12.3,3f9.2)') Nmain, X0, Y0 , Z0,aM !,aM,rr,vnow,vcrc
45          Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4)

        If(iStat.ne.0)Exit
        If(vcrc> vcLimit)Then
        Nmain = Nmain +1
        If(iFlag==1)Then
           if(Nmain>Ncenter)Stop ' Not enough space. Increase Ncenter'
              Xc(Nmain) = X0
              Yc(Nmain) = Y0
              Zc(Nmain) = Z0
              VXc(Nmain) = VvX
              VYc(Nmain) = VvY
              VZc(Nmain) = VvZ
              Amc(Nmain) = aM200
              Rmc(Nmain) = rr/1000.
              Vcirc(Nmain) = max(vcrc,V200)/V200
              !Vrms(Nmain) = vrmss
              Rdmax  = MAX(rr/1000.,Rdmax)
        EndIf
        EndIf     !  test mass
   Enddo
40          Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3)

    If(iFlag==0)write(*,*) 'Read Halos=',Nmain
    If(iFlag==1)write(*,'(a,i8,a,f8.3,a)') ' Read Halos=',Nmain, &
                                            ' Max Radius=',Rdmax,'Mpch'
    rewind(3)
      Return
10    write(*,*) ' Could not read radial bins'
      stop
20    write(*,*) ' Could not read text R(kpc/h)'
      stop
30    write(*,*) ' Unexpected End of file'
      stop

      End SUBROUTINE ReadHalos
 !-------------------------------------------------------------------------
!                    Make linker list of Halos
      SUBROUTINE ListHalos
!------------------------------------------------------------------------
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
             Do i=1,Ncenter
                Lst(i)= -1
             EndDo

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i,j,k)
       Do k=Nmz,Nbz
       Do j=Nmy,Nby
       Do i=Nmx,Nbx
          Label(i,j,k)=0
       EndDo
       EndDo
       EndDo  
      Do jp=1,Nsample
         i=INT(Xc(jp)/Cell+100)-100
         j=INT(Yc(jp)/Cell+100)-100
         k=INT(Zc(jp)/Cell+100)-100
         i=MIN(MAX(Nmx,i),Nbx)
         j=MIN(MAX(Nmy,j),Nby)
         k=MIN(MAX(Nmz,k),Nbz)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo
      write(*,*) ' Done Linker List'
      Return
      End SUBROUTINE ListHalos

!     ------------------------
      real function seconds ()
!     ------------------------
!
!     purpose: returns elapsed time in seconds
      Integer, SAVE :: first=0,rate=0
      Real   , SAVE :: t0=0.
      !real*8 dummy
      !real tarray(2)
      If(first==0)Then
         CALL SYSTEM_CLOCK(i,rate)
         first =1
         t0    =float(i)/float(rate)
         seconds = 0.
      Else
!        seconds = etime(tarray)
         CALL SYSTEM_CLOCK(i)
         seconds = float(i)/float(rate)-t0
      EndIf
      return
      end function seconds


end module SetHalos
!--------------------------------------------------------------------------------------
!
!     
PROGRAM HaloSatellites

use     SetHalos 
Character*120 :: file1,file2,catname,listname,txt*8     

         write(*,'(a,$)')' Enter snapshot number  ='
         read(*,*)iSnap
         write(*,'(a,$)')' Enter expansion parameter  ='
         read(*,*)AEXPN
         write(file1,'(a,i3.3,a)')'CatshortDistinctA200.',iSnap,'.DAT'
         write(file2,'(a,i3.3,a)')'HaloMass.',iSnap,'.DAT'

      Open( 3,file=TRIM(file1),Status='OLD')! short list of halos
      Open( 2,file=TRIM(file2))                         ! output

         !Box =250.
         !Cell = 2.
         !dBuffer = 2.
         !vcLimit =30.
         Box =1000.
         Cell = 5.
          dBuffer = 3.
          vcLimit = 250.
         !write(*,'(a,$)')' Enter limit on circular velocity  ='
         !read(*,*)vcLimit
      CALL ReadHalos(0,vcLimit)                                          ! count halos
      fExpand = (1.+2.*dBuffer/Box)**3*1.1             ! estimate how much more space
                                                                                    ! will be needed 
      Ncenter = Nmain*fExpand                                ! make more centers for periodical conditions
      write(*,*) ' fExpand =',fExpand 
      write(*,*) ' Maximum of halos Ncenter =',Ncenter

      ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter))
      ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter))
      ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter))
      ALLOCATE (Lst(Ncenter))
      Nmx =0 ; Nbx = Box/Cell +1
      Nmy =0 ; Nby = Box/Cell +1
      Nmz =0 ; Nbz = Box/Cell +1
      ALLOCATE (Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz))
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
        Do i=1,Ncenter
           Xc(i) =0. ;Yc(i) =0. ;Zc(i) =0. 
           VXc(i) =0. ; VYc(i) =0. ;  VZc(i) =0.
           Amc(i) =0. ; Rmc(i) =0. ; Vrms(i) =0.
           Vcirc(i) =0.
        EndDo
         write (*,*) 'Allocated Centers'
      CALL ReadHalos(1,vcLimit)                                          ! read halos
      Ncenter = Nmain*fExpand                       ! restore Ncenter

      !Call Expand                                               ! make more halos by adding halos
                                                                          ! due to periodical conditions
      !CALL ListHalos
      Close(3)

       ALLOCATE (LDistinct(Ncenter))
       LDistinct   = 1
      !Call     HaloDistinct    ! make list of distinct halos 
         write (*,*) 'Found Distinct halos'
      !Call     Satellites                                         ! make list of satellites
      !   write (*,*) 'Found Satellites'
      ! Call     MassDistinct                                         ! make list of satellites
         write (*,*) 'Found Mass distribution'
       Call VcircMvir                                            ! Vcirc-Mvir relation for Distinct and Subs

End PROGRAM HaloSatellites
Subroutine R_mrgref (XVALT, IRNGT,N)
!   Ranks array XVALT into index array IRNGT, using merge-sort
! __________________________________________________________
!   This version is not optimized for performance, and is thus
!   not as difficult to read as some other ones.
!   Michel Olagnon - April 2000
! __________________________________________________________
! _________________________________________________________
      Real, Dimension (N), Intent (In) :: XVALT
      Integer, Dimension (N), Intent (Out) :: IRNGT
! __________________________________________________________
!
      Integer, Dimension (:), Allocatable :: JWRKT
      Integer :: LMTNA, LMTNC
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XVALT), SIZE(IRNGT))
      If (NVAL <= 0) Then
         Return
      End If
!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XVALT(IIND-1) < XVALT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo (NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      Allocate (JWRKT(1:NVAL))
      LMTNC = 2
      LMTNA = 2
!
!  Iteration. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
         IWRK = 0
!
!   Loop on merges of A and B into C
!
         Do
            IINDA = IWRKF
            IWRKD = IWRKF + 1
            IWRKF = IINDA + LMTNC
            JINDA = IINDA + LMTNA
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDB = JINDA
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B (no need to do anything)
!
            If (XVALT(IRNGT(JINDA)) <= XVALT(IRNGT(JINDA+1))) Then
               IWRK = IWRKF
               Cycle
            End If
!
!  One steps in the C subset, that we create in the final rank array
!
            Do
               If (IWRK >= IWRKF) Then
!
!  Make a copy of the rank array for next iteration
!
                  IRNGT (IWRKD:IWRKF) = JWRKT (IWRKD:IWRKF)
                  Exit
               End If
!
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (IINDA < JINDA) Then
                  If (IINDB < IWRKF) Then
                     If (XVALT(IRNGT(IINDA+1)) > XVALT(IRNGT(IINDB+1))) &
                    & Then
                        IINDB = IINDB + 1
                        JWRKT (IWRK) = IRNGT (IINDB)
                     Else
                        IINDA = IINDA + 1
                        JWRKT (IWRK) = IRNGT (IINDA)
                     End If
                  Else
!
!  Only A still with unprocessed values
!
                     IINDA = IINDA + 1
                     JWRKT (IWRK) = IRNGT (IINDA)
                  End If
               Else
!
!  Only B still with unprocessed values
!
                  IRNGT (IWRKD:IINDB) = JWRKT (IWRKD:IINDB)
                  IWRK = IWRKF
                  Exit
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!
!  Clean up
!
      Deallocate (JWRKT)
      Return
!
End Subroutine R_mrgref