!-------------------------------------------------------------------- 
!             Find isolated  halos in a given mass range in Catalog.DAT
!             Construct their averaged density and vel profile
!-------------------------------------------------------------------- 
PROGRAM HaloProfiles
      include 'haloprofiles.h' 
      CHARACTER  Name*80,sName*2

!      Open(1,file='Catalog.unbound.DAT')
      Open(1,file='Halolist.a1.0003.DAT')
      Vc1 =250.
      Vc2 =300.
      write (*,'(A,$)')' Enter 2 characters  for file name=>'
      read (*,*) sName
         write(Name,'(3a,i3.3,a)')'results.',sName,'.',INT(Vc2),'.dat'
      Open(2,file=Name)
 
!                                                   Write headers
      write (*,'(T10,"Hubble=",f8.2,T25,"Box=",f8.2,"Mpc/h", &
                   T35,"Vc limits=",3f8.2,"km/s")')         &
                   hubble,Box,Vc1,Vc2,Vclimit 

      Lines =7                       ! skip header of catalog file
      Do i=1,Lines
         read(1,'(A)')Line
            if(i.le.40)write(*,*)Line
      EndDo 

      Nh = 1                           ! read halo catalog
 10   read(1,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh),       &
                     vxc(Nh),vyc(Nh),vzc(Nh),aMass(Nh),   &
                     Rvir(Nh),Vrms(Nh),vc(Nh),nn,Rmax,conc  !,Nbins(Nh)
       ! write(*,'(i7,3f9.3,i4)') Nh,xc(Nh),yc(Nh),zc(Nh)    !,Nbins(Nh)
         Nbins(Nh) =100
        Do j=1,Nbins(Nh)
          read(1,*)rad(j,Nh),nn,aM(j,Nh),overd,Vdisp(j,Nh), &
          Vcirc(j,Nh),Dens(j,Nh),i1,i2,i3,Vrad(j,Nh)
        EndDo
      if(vc(Nh).lt.Vclimit)goto 10  ! do not take halo if too small
      if(Nh.ge.Nmax)goto 100
      Nh =Nh +1
      goto 10

 100  close(1)
      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 ExpandSample
      CALL Isolated3D(Vc1,Vc2)
      CALL GetProfiles

      stop
end PROGRAM HaloProfiles
!-------------------------------------------------------------------- 
!                        make sample larger by replicating
!                        some halos using periodical
!                        boundary conditions
!          dR is the width of the buffer in Mpc/h
!          Nhext = new number of halos
      SUBROUTINE ExpandSample
!-------------------------------------------------------------------- 
      include 'haloprofiles.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
      Nhext = 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
                     Nhext =Nhext +1
                     if(Nhext.gt.Nmax) STOP 'too many halos'
                     xc(Nhext)    = x+dx
                     yc(Nhext)    = y+dy
                     zc(Nhext)    = z+dz
                     vxc(Nhext)   = vxc(ip)
                     vyc(Nhext)   = vyc(ip)
                     vzc(Nhext)   = vzc(ip)
                     aMass(Nhext) = aMass(ip)
                     vc(Nhext)    = vc(ip)
                     Rvir(Nhext)  = Rvir(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=',Nhext  
      write (2,'("  Nhalos in extended sample=",T30,i8)')Nhext  
      RETURN 
END   SUBROUTINE ExpandSample

!-------------------------------------------------------------------- 
!                  select isolated halos
!                  range of circular velocities Vc1<Vc< Vc2
!                  do it in 3D, not in projection
!                  Liso = 1 if halo is isolated
      SUBROUTINE Isolated3D(Vc1,Vc2)
!-------------------------------------------------------------------- 
      include 'haloprofiles.h' 

      Nx = 0
      aMav  =0.
      rvirav=0.
      vcav  =0.
      R2     = (RadLmax)**2  ! square of max distance in Mpc
 
      Do ip =1,Nh            ! Initialize counters
        vv =vc(ip)
        If(vv.gt.Vc1.and.vv.lt.Vc2)Then
          Liso(ip) = 1
          Nx       = Nx +1
        Else
          Liso(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(Liso(ip).eq.1)Then
            x =xc(ip)
            y =yc(ip)
            z =zc(ip)
            Do jp =1,Nhext
               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)Liso(ip) =0
              EndIf                  ! jp=ip
            EndDo   ! jp loop
         EndIf         ! Liso(ip) =1
      EndDo      ! ip loop

      Do ip =1,Nh
         If(Liso(ip).eq.1)Then
            Nx    =Nx     +1
            aMav  =aMav   +aMass(ip)
            rvirav=rvirav +Rvir(ip)
            vcav  =vcav   +vc(ip)
         EndIf 
      EndDo
      write(*,'(" Nisolated=",3i7)') Nx
      write(2,'("  Nisolated halos =",T20,i4," Vc limits=",2f8.2,  &
                3x,"<Mvir>=",g12.4," rvir=",f8.3," <Vc>=",f8.3)')  &
                   Nx,Vc1,Vc2, aMav/Nx, rvirav/Nx, vcav/Nx
      RETURN 
End SUBROUTINE Isolated3D
!-------------------------------------------------------------------- 
!                  Cosnstruct averaged profiles of isolated halos
!
SUBROUTINE GetProfiles
!-------------------------------------------------------------------- 
    include 'haloprofiles.h' 
    DIMENSION Ss(Nshells),Slope(Nshells)
    Do i=1,Nshells
        rdsh(i) =0.
        aMsh(i) =0.
        Vdsh(i) =0.
        Vcsh(i) =0.
        Dnsh(i) =0.
        Vradsh(i) =0.
        Ss(i)     =0.
        Slope(i)  =0.
    EndDo
    Nx     =0           ! count halos. Find number of bins Nb
    aMav   =0.
    rvirav =0.
    vcav   =0.
    Do ip =1,Nh
       If(Liso(ip).eq.1)Then
          Nx    =Nx     +1
          aMav  =aMav   +aMass(ip)
          rvirav=rvirav +Rvir(ip)
          vcav  =vcav   +vc(ip)
          If(Nx.eq.1)Then 
             Nb = Nbins(ip)
          Else   
            If(Nbins(ip).ne.Nb)Then
              write(*,*)' Error:Number of radial bins not constant'
              stop
            EndIf
          EndIf 
       EndIf 
    EndDo
    If(Nx.eq.0)Then
        write(*,*) ' No isolated halos. Stop'
        stop
    EndIf
    aMav   =aMav  /Nx
    rvirav =rvirav/Nx
    vcav   =vcav  /Nx

    Do ip =1,Nh         ! get profiles
      If(Liso(ip).eq.1)Then
        Do j=1,Nb   
          rdsh(j)  =rdsh(j) + rad(j,ip)
          aMsh(j)  =aMsh(j) + aM(j,ip)
          Dnsh(j)  =Dnsh(j) + Dens(j,ip)
          Vdsh(j)  =Vdsh(j) + sqrt(max(Vdisp(j,ip)**2-Vrad(j,ip)**2,1.e-10))
          Vradsh(j) =Vradsh(j) +(Vrad(j,ip)+100.*rad(j,p)/1000.)
        EndDo
      EndIf 
    EndDo

   Do j=1,Nb   
     rdsh(j)  =rdsh(j)/Nx
     aMsh(j) =aMsh(j)/Nx
     Dnsh(j)  =Dnsh(j)/Nx
     Vdsh(j)  =Vdsh(j)/Nx
     Vradsh(j)=Vradsh(j)/Nx
   EndDo

   Ss =Dnsh                ! smooth density profile
   Do j=2,Nb-1   
     Ss(j)  =(Dnsh(j-1)+Dnsh(j)+Dnsh(j+1))/3.
   EndDo                   

   Do j=2,Nb-1             ! find slope
     Slope(j)  =log10(Ss(j+1)/Ss(j-1))/log10(rdsh(j+1)/rdsh(j-1))
   EndDo
    Slope(1) =Slope(2)
    Slope(Nb) =Slope(Nb-1)

   Ss =Slope                ! smooth slope
   Do j=2,Nb-1   
     Ss(j)  =(Slope(j-1)+Slope(j)+Slope(j+1))/3.
   EndDo                   ! find slope
   Slope = Ss

   write(2,'(3x,a,T15,a,T30,a,T45,a,T60,a,T75,a)')  &
       'Rad(kpch)','Mass','OverDens','Vrms','Vrad','Slope'
   Do j=1,Nb 
     write(2,'(f8.3,T15,G12.3,T30,G12.3,T45,f8.2,T60,f8.2,T75,f8.2)')&
       rdsh(j),aMsh(j),Dnsh(j), Vdsh(j),Vradsh(j),Slope(j)
   EndDo

    RETURN 
End  SUBROUTINE GetProfiles