!----------------------------------------------------------------------------------------
!                    3d properies of halos
!                    
!     Input:  Halo Catalog (long list)
!                 List of distinct halos
!     Output:  average profiles of halos in a range of Vcirc 
!     Parameters:
!    
!    
MODULE   setHalos
PUBLIC      

Real,   PARAMETER    ::                                        &
                         dLog              = 0.1,                    & ! bin size in log(r)
                         dMlog            =0.15
Integer, PARAMETER ::     Nbin = 27
  Real,    DIMENSION (Nbin)     :: Nin,Rad_out,Rad_bin,Mr,Dens, &
                                    Vcirc,Vrms,Vrad,Vrad2,Jx,Jy,Jz,Jtot

Integer :: Ncenter,Nmain,Ntotal
Real      :: Box,  Risomax, Vcmin,Vcmax,Vfrac, dBuffer
                              
Contains


!---------------------------------------------------------
!                         Read halos catalog
!                       
      SUBROUTINE GetHalos(Vlow,Vup)
!--------------------------------------------------------------
Character  :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4
Logical :: take
   Ncenter  = 0;   Rvir_a     = 0.;    aM_a      =  0.
   Vc_a       = 0.;  aL_a        = 0.;    x_a          = 0.
   vir_a         = 0.
   Nin     = 0; Rad_out =0.; Rad_bin = 0.; Mr =0.
   Dens = 0.; Vcirc =  0.; Vrms = 0.; Vrad = 0.; Vrad2 = 0.
   Jx       = 0.; Jy      = 0.; Jz       = 0.; Jtot = 0.
   Do i=1,11
      Read(3,'(a)')Line
      If(i<10)write(2,'(a)')Line
   EndDo

   Do                 !kk=1,100000
     read (3,*,iostat=iStat) j,X0, Y0 , Z0, VvX,VvY,VvZ,   &
                                               aM,rr,vrmss,vcrc,aMunb,vrms_u,VirRat,x_off,aJ,aLambda
     !write (*,'(2i9,3f10.4,g12.4,2f9.2)') kk,Ncenter, X0, Y0 , Z0,aM,vcrc,aLambda
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 > Vlow .and. vcrc < Vup.and.abs(VirRat)<0.5)Then
        If(vcrc > Vlow .and. vcrc < Vup)Then
          take = .true. 
          Ncenter = Ncenter +1
          aM_a = aM_a + aMunb
          Vc_a  = Vc_a + vcrc
          aL_a   = aL_a + aLambda
          x_a     = x_a    + x_off
          vir_a   = vir_a  + VirRat
        Else
           take = .false.
        EndIf
          Do i =1,Nbin
             read(3,*,iostat =iStat)n,r1,r2,am1,d1,vc1,vs1,vr1,vr2,ajx1,ajy1,ajz1,aj
             If(take)Then
                Nin(i)           = Nin(i)          + n
                Rad_out(i) = Rad_out(i) + r1
                Rad_bin(i) = Rad_bin(i)  + r2
                Mr(i) = Mr(i)  + am1
                Dens(i) = Dens(i)  + d1 
                Vcirc(i) = Vcirc(i)    + vc1
                Vrms(i) = Vrms(i)  + vs1**2 
                Vrad(i) = Vrad(i)   + vr1
                Vrad2(i) = Vrad2(i)  + vr2 **2
                Jtot(i) = Jtot(i)  + aj**2 
             EndIf
        EndDo
   Enddo

   Do j=0,2,2
    write(j,'(a,i6,a,2f8.1)') ' Read Halos =',Ncenter,' Vc   =',Vlow,Vup
    write(j,'(a,g13.4)')      &
                       ' Mass         =',aM_a/Ncenter,      &
                       ' Vcirc        =',Vc_a/Ncenter,                &
                       ' Lambda       =',aL_a/Ncenter,       &
                       ' X_offset     =', x_a/Ncenter,  &
                       '  2-Epot/Ekin =',vir_a/Ncenter
    EndDo
    write(2,'(T4,a,T12,a,T25,a,T36,a,T44,a,T59,a,T67,a,T76,a,T85,a,T96,a)')    &
                   'Npart','R_out/Rvir','R_av/Rvir','M(<r)','Dens(Mpc/kpc3)',     &
                   'Vcirc','Vrms_3d','Vradial','Vrms_rad','AngMom'      
    n = Ncenter
    Do i =1,Nbin
       write(2,'(f9.1,2f11.4,1p,2g12.4,0p,8f9.3)')(Nin(i))/float(n),Rad_out(i)/n,Rad_bin(i)/n,     &
                          Mr(i)/n,Dens(i)/n,Vcirc(i)/n, sqrt(Vrms(i)/n),                                 &
                          Vrad(i)/n,sqrt(Vrad2(i)/n),sqrt(Jtot(i)/n)
    EndDo

      End SUBROUTINE GetHalos

!     ------------------------
      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 HaloProfiles

use     SetHalos 
Character*120 :: file1 

      Open( 3,file='CatalogDistinctB.416.DAT',Status='OLD')    !  list of halos
      write(*,'(a,$)') ' Enter limits on circular velocity:'
      read(*,*) Vlow,Vup
      write(file1,'(a,i4.4,a,i4.4,a)') 'HaloProfileB.',INT(Vlow),'.',INT(Vup),'.dat'
      Open( 2,file=file1)                                 ! output file

 CALL GetHalos(Vlow,Vup)

End PROGRAM HaloProfiles