C ______________________________	 Selects a halo from a large Catalog.DAT
C	
C
C          Scales internal PM coordinates and velocities 
C          to  coordinates in Mpc/h and km/s
C          In PM the coordinates are in the range 1 - (NGRID+1)
C                     velocities are P = a_expansion*V_pec/(x_0H_0)
C                     where     x_0 = comoving cell_size=Box/Ngrid
C                                    H_0 = Hubble at z=0
C                   
C	     NROW = number of particles in 1D
C	     NGRID= number of cells        in 1D
      INCLUDE 'PMselect.h'
      REAL     INPUT
      Character*50  FileASCII
      Logical    Inside

C...................................................................
C			Read data and open files
      WRITE (*,'(A,$)') 'Enter Name for output file = '
      READ  (*,'(A)') FileASCII
      OPEN(17,FILE=FileASCII,STATUS='UNKNOWN')
      WRITE (*,'(A,$)') 'Enter Name of input file = '
      READ  (*,'(A)') FileASCII
      OPEN(20,FILE=FileASCII,STATUS='UNKNOWN')
      WRITE (*,'(A,$)') 'Enter Name of catalog with',
     &     ' corrected concentrations file = '
      READ  (*,'(A)') FileASCII
      OPEN(21,FILE=FileASCII,STATUS='UNKNOWN')

      Box  =INPUT(' Enter box size in comoving Mpc/h =')
      dBox = 10.  ! expand box by dBox Mpc/h
      write (*,'(A,$)')' Enter selection Limits: ',
     +     'Mass dM X Y Z dR C dC =>'
      read  (*,*) aMass, dM, X0, Y0, Z0, dR, Cc, dCc

      CALL ReadCatalog(Nh,Nhex)
      CALL ExpandSample(Nh,Nhex,dBox,Box)
      write(17,30)aMass-dM,aMass+dM,Cc-dCC,Cc+dCC,aMiso,Riso
c      CALL PrintSelected(Nh,Nhex)
      CALL AverageSelected(Nh,Nhex)
 30   format(' Mass range=',2g11.4,' Concentrations=',2g11.3,
     &  /' Vcirc_neighbor/Vcirc < ',f8.2,' Rneighbor/Rvir> ',f8.2)
      END
c-------------------------------------------------------------------- 
      SUBROUTINE AverageSelected(Nh,Nhex)
c-------------------------------------------------------------------- 
      INCLUDE 'PMselect.h'
      Logical inside,isolated
      COMMON /AVER/ Rapr(Nrmax),Napr(Nrmax),
     +             Vacpr(Nrmax),RHOapr(Nrmax),Varpr(Nrmax),
     +             bMa(Nrmax),Dapr(Nrmax),Vadpr(Nrmax),
     &             
      
      write (*,*)  ' Concentration limits=',Cc-dCc,Cc+dCc
      nn = 0
      aMav =0.
      rav  =0.
      Do i=1,Nh
         amm = aM(i)
         Cnc = Conc(i)
         x   = Xp(i)
         y   = Yp(i)
         z   = Zp(i)
         inside = (amm.gt.aMass-dM.and. amm.lt.aMass+dM )
         inside = (Cnc.gt.Cc-dCc.and.Cnc.lt.Cc+dCc.and.inside)
         inside = (x.gt.X0-dR.and. x.lt.X0+dR .and.inside)
         inside = (y.gt.Y0-dR.and. y.lt.Y0+dR .and.inside)
         inside = (z.gt.Z0-dR.and. z.lt.Z0+dR .and.inside)
         inside = (errFit(i).lt.errLimit.and.inside)
         If(inside)Then
          CALL CheckIsolated(i,Nh,Nhex,isolated)
c          isolated =.true.
          If(isolated)Then  
           write(*,50) Xp(i),Yp(i),Zp(i),
     &         VXp(i),VYp(i),VZp(i),
     &         aM(i),Radh(i),Vrms(i),Vcirc(i),
     &         Npart(i),Rmax(i),Conc(i),Nlines(i)
            nn = nn +1
            aMav = aMav + aM(i)
            rav  = rav  + Radh(i)
             If(nn.eq.1)Then
                ifirst = i
               Do j=1,Nlines(i)
                  Rapr(j)  = Rpr(j,i)/ Radh(i)
                  Vacpr(j) = Vcpr(j,i)
                  RHOapr(j)= RHOpr(j,i)
                  bMa(j)   = RHOpr(j,i)**2
                  Varpr(j) = Vrpr(j,i)
               EndDo
             else 
                ilast = i
               Do j=1,Nlines(i)
                  Rapr(j)  = Rapr(j)  +Rpr(j,i)/ Radh(i)
                  Vacpr(j) = Vacpr(j) +Vcpr(j,i)
                  RHOapr(j)= RHOapr(j)+RHOpr(j,i)
                  bMa(j)   = bMa(j)   +RHOpr(j,i)**2
                  Varpr(j) = Varpr(j) +Vrpr(j,i)
               EndDo
             EndIf 
          EndIf 
         EndIf 
      EndDo
      If(nn.eq.0)Then
         write (*,*)  ' no halos were selected'
         return
      EndIf 
      write (*,*)  nn,' halos were selected'
               Do j=1,Nlines(ifirst)
                  Rapr(j)  = Rapr(j)/nn
                  Vacpr(j) = Vacpr(j)/nn
                  RHOapr(j)= RHOapr(j)/nn
                  bMa(j)   = sqrt(bMa(j)/nn -RHOapr(j)**2)
                  Varpr(j) = Varpr(j)/nn
               EndDo
            aMav = aMav/nn
            rav  = rav/nn

          write(17,'(" Nselected=",i6," Mass=",g11.3," Rvir=",g11.3)')
     &              nn,aMav,rav
          write(17,50) Xp(ilast),Yp(ilast),Zp(ilast),
     &         VXp(ilast),VYp(ilast),VZp(ilast),
     &         aM(ilast),Radh(ilast),Vrms(ilast),Vcirc(ilast),
     &         Npart(ilast),Rmax(ilast),Conc(ilast),Nlines(ilast)
          Do j=1,Nlines(ilast)
            write(17,60) Rapr(j),Npr(j,ilast),
     &                bM(j,ilast),Dpr(j,ilast),Vdpr(j,ilast),
     &          Vacpr(j),RHOapr(j)*Dpr(1,1)/RHOpr(1,1),
     &          Npr(j,ilast),0,0,Varpr(j),bMa(j)*Dpr(1,1)/RHOpr(1,1)
          EndDo 

 40   Format(2i7,3g11.3,g12.4,f9.2)
 50   Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4)
 60   format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,2g11.3)
      Return
      End
c-------------------------------------------------------------------- 
      SUBROUTINE PrintSelected(Nh,Nhex)
c-------------------------------------------------------------------- 
      INCLUDE 'PMselect.h'
      Logical inside,isolated

      nn = 0
      Do i=1,Nh
         amm = aM(i)
         Cnc = Conc(i)
         x   = Xp(i)
         y   = Yp(i)
         z   = Zp(i)
         inside = (amm.gt.aMass-dM.and. amm.lt.aMass+dM )
         inside = (Cnc.gt.Cc-dCc.and.Cnc.lt.Cc+dCc.and.inside)
         inside = (x.gt.X0-dR.and. x.lt.X0+dR .and.inside)
         inside = (y.gt.Y0-dR.and. y.lt.Y0+dR .and.inside)
         inside = (z.gt.Z0-dR.and. z.lt.Z0+dR .and.inside)
         inside = (errFit(i).lt.errLimit.and.inside)
         If(inside)Then
          CALL CheckIsolated(i,Nh,Nhex,isolated)
c          isolated =.true.
          If(isolated)Then  
             nn = nn +1
          write(*,40) nn,i,Conc(i),aM(i),Vcirc(i)
          write(17,50) Xp(i),Yp(i),Zp(i),
     &         VXp(i),VYp(i),VZp(i),
     &         aM(i),Radh(i),Vrms(i),Vcirc(i),
     &         Npart(i),Rmax(i),Conc(i),Nlines(i)
c          Do j=1,Nlines(i)
c            write(17,60) Rpr(j,i)/Radh(i),Npr(j,i),
c     &                bM(j,i),Dpr(j,i),Vdpr(j,i),
c     &          Vcpr(j,i),RHOpr(j,i)*Dpr(1,i)/RHOpr(1,i),
c     &          Npr(j,i),0,0,Vrpr(j,i)
c          EndDo 
         EndIf 
         EndIf 
      EndDo 
      write (*,*)  nn,' halos were selected'
 40   Format(2i7,3g11.3,g12.4,f9.2)
 50   Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4)
 60   format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3)
      Return
      End
c-------------------------------------------------------------------- 
      SUBROUTINE CheckIsolated(i,Nh,Nhex,isolated)
c-------------------------------------------------------------------- 
      INCLUDE 'PMselect.h'
      Logical isolated

      isolated = .true.
      amm      = Vcirc(i) ! aM(i)
c      amm      = aM(i)
      Rhalo    = Radh(i)
         x   = Xp(i)
         y   = Yp(i)
         z   = Zp(i)

      Do j=1,Nhex
         If(i.ne.j.and.amm.lt.Vcirc(j)*aMiso)Then
            dd = (Xp(j)-x)**2+(Yp(j)-y)**2+(Zp(j)-z)**2
            diso = Riso*max(Rhalo,Radh(j))/1.e3  ! put it in Mpc
c            write (*,'(2i7,10g11.3)')i,j,sqrt(dd),diso,Rhalo,
c     &            Radh(j),Vcirc(i),Vcirc(j)
            If(dd.lt.diso**2)Then
               isolated = .false.
               return
            EndIf 
         EndIf 
      EndDo 
      If(Vrms(i).gt.1.4*Vcirc(i))write (*,*)' Unbound halo:',
     &                           i,Vrms(i),Vcirc(i)
      Return
      End
c-------------------------------------------------------------------- 
      SUBROUTINE ReadCatalog(Nh)
c-------------------------------------------------------------------- 
      INCLUDE 'PMselect.h'
      Character*120 Line

      Lines = 27

      Do i=1,Lines
         read(20,'(A)')Line
         read(21,'(A)')Line
         write(*,'(A)')Line
         write(17,'(A)')Line
      EndDo 
      Nh = 1

  10   read(20,*,end=100,err=100) Xp(Nh),Yp(Nh),Zp(Nh),
     &         VXp(Nh),VYp(Nh),VZp(Nh),
     &         aM(Nh),Radh(Nh),Vrms(Nh),Vcirc(Nh),
     &         Npart(Nh),Rmax(Nh),Conc(Nh),Nlines(Nh)

       read(21,*,end=100,err=100) Xp0,Yp0,Zp0,
     &         VXp0,VYp0,VZp0,
     &         aM0,Radh0,Vrms0,Vcirc0,
     &         Npart0,Rmax0,Conc0,Nlines0,errFit(Nh)

       If(Nh.lt.10.or.Nh/5000*5000.eq.Nh)
     &         write (*,'(" header=",i7,g11.3,i6,3g12.4)')
     &         Nh,aM(Nh),Nlines(Nh),Xp(Nh),Yp(Nh),Zp(Nh)
       If(Nh.lt.10.or.Nh/5000*5000.eq.Nh)
     &         write (*,'(" header=",i7,g11.3,i6,3g12.4)')
     &         Nh,aM0,Nlines0,Xp0,Yp0,Zp0

       If(abs(Xp(Nh)-Xp0).gt.1.e-2.or.abs(Yp(Nh)-Yp0).gt.1.e-2
     &    .or.abs(Zp(Nh)-Zp0).gt.1.e-2)stop 'Mismatch between catalogs'
       If(abs(aM(Nh)-aM0).gt.1.e-2*aM0.or.Nlines(Nh).ne.Nlines0)
     &       stop 'Mismatch between catalogs'
       Conc(Nh) =Conc0
       If(Nlines(Nh).gt.Nrmax) STOP ' wrong number of profile lines'
       If(Nlines(Nh).le.0) STOP ' wrong number of profile lines'
       Do i=1,Nlines(Nh)
          read (20,*) Rpr(i,Nh),Npr(i,Nh),
     &                bM(i,Nh),Dpr(i,Nh),Vdpr(i,Nh),
     &                Vcpr(i,Nh),RHOpr(i,Nh),
     &                nsm,nmed,nlg,Vrpr(i,Nh)
       EndDo

      if(Nh.ge.Np)goto100
      Nh =Nh +1
      goto 10

 100  Nh =Nh -1

      write (*,'("  Number of halos read=",T30,i8)')  Nh

      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(Nh,Nhex,dBox,Box)
c-------------------------------------------------------------------- 
      INCLUDE 'PMselect.h'
      Logical inx,iny,inz

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

      Do ip =1,Nh    ! go over all halos
         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
                     Next =Next +1
                     if(Next.gt.Np) STOP 'too many halos'
                     Xp(Next) =x+dx
                     Yp(Next) =y+dy
                     Zp(Next) =z+dz
                     VXp(Next) =VXp(ip)
                     VYp(Next) =VYp(ip)
                     VZp(Next) =VZp(ip)
                     aM(Next)  =aM(ip)
                     Radh(Next) =Radh(ip)
                     Vrms(Next) =Vrms(ip)
                     Vcirc(Next)=Vcirc(ip)
                     Conc(Next) =Conc(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 
      Nhex =Next
      RETURN 
      END

C---------------------------------- Read in variables      
      REAL FUNCTION INPUT(text)
C------------------------------------------------
      Character text*(*)
          write (*,'(A,$)')text
          read (*,*) x
          INPUT =x
      Return
      End
