!----------------------------------------------------------------------------------------
!                    3d properies of halos
!                    
!     Input:  Halo Catalog (long list)
!                 List of distinct halos
!     Output:  short list catalog with V200 and M200
!     Parameters:
!    
!    
MODULE   setHalos
PUBLIC      

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

Integer :: Ncenter,Nmain,Ntotal,Ncurrent
Real      :: Box,Vcmax,xc(3),vv(3),aM,rVir,VcircMvir_unbnd, AEXPN, &
                 Vrms_unbnd,VirRat, dRRvir, AngMomLambda,M200
                              
Contains


!---------------------------------------------------------
!                         Read halos catalog
!                       
      SUBROUTINE GetHalos
!--------------------------------------------------------------
Character  :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4
   Do i=1,11
      Read(3,'(a)')Line
      if(i<10)write(2,'(a)')Line
   EndDo
   write(2,'(2a)')'Npart            Coords(Mpch)              Veloc(km/s)         Mvir         Rvir     Vrms     Vcirc ', &
                          '    Mvir_unbnd Vrms      VirRat   dR/Rvir AngMom Lambda      R200       V200      M200'
   Om0 = 0.27
   Const =  8.64e-4/(Om0+AEXPN**3*(1.-Om0))
   Overdens = 200.

   write (*,*) ' -------'

   Do          !  kk=1,3000
     read (3,*,iostat=iStat) &
                                               Ncurrent,xc,vv,  &
                                               aM,rVir,vrmss,vcrc,aMunb,vrms_u,VirRat,x_off,aJ,aLambda
     !write (*,'(i9,3f10.4,g12.4,2f9.2)') Ncurrent, xc,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
          Do i =1,Nbin
             read(3,*,iostat =iStat)n,r1,r2,am1,d1,vc1,vs1,vr1,vr2,ajx1,ajy1,ajz1,aj 
                Nin(i)           = n
                Rad_out(i) = r1*rVir
                Rad_bin(i)  = r2*rVir
                Mr(i)            =  am1
                Dens(i)       =  d1 
                Vcirc(i)        = vc1
                dd = Mr(i)/Rad_out(i)**3*Const                
             !write (*,'(2i6,f8.3,f8.2,3g12.4)') i,n,Rad_out(i),Vcirc(i),dd,Mr(i) 
                If(r1>0.9999)exit
        EndDo
        ilast =i
        read(3,*,iostat =iStat)n,r1,r2
        If(abs(r1-1.1660)>1.e-3 .and.abs(r2-1.000)>1.e-3)backspace (3)
   Do i=ilast,1,-1
      dd = Mr(i)/Rad_out(i)**3*Const
      If(dd > Overdens)Then
         If(i == ilast)Then
            R200 = Rad_out(i)
            V200 = Vcirc(i)
            M200 = Mr(i)
         Else
            d0      = dd
            d1      = Mr(i+1)/Rad_out(i+1)**3*Const
            R200 =  Rad_out(i+1) - (Overdens -d1)*(Rad_out(i+1)-Rad_out(i))/(d0-d1)
            V200 = Vcirc(i+1) + (Vcirc(i)-Vcirc(i+1))*(Rad_out(i+1)-R200)/(Rad_out(i+1)-Rad_out(i))
            M200 = Mr(i+1) + (Mr(i)-Mr(i+1))*(Rad_out(i+1)-R200)/(Rad_out(i+1)-Rad_out(i))
         EndIf
         !If(M200/aM<3. .and. aM>3.e12)&
         !write(*,'(i4,10g12.4)') i,Mr(i),Rad_out(i),Vcirc(i),dd,R200,V200,M200
         Exit
       EndIf
    EndDo
       !  If(M200/aM<1./3. .and. aM>3.e12)Then
       !  write(*,'(10g12.4)') aM,rVir,vcrc,M200,R200,V200
   !Do i=1,ilast
   !   dd = Mr(i)/Rad_out(i)**3*Const
   !   write(*,'(i3,7g12.4)') i,Rad_out(i),Mr(i),dd,Vcirc(i)
   !EndDo
   !      EndIf

       write(2,'(i7,3f9.4,3f9.2,2(g12.4,3f9.2),f9.3,3x,f9.4,f8.2,1p2g12.4)')Ncurrent,xc,vv,  &
                                               aM,rVir,vrmss,vcrc,aMunb,vrms_u,VirRat,x_off,aJ,aLambda, &
                                               R200,V200,M200
   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 

      write(*,'(a,$)')' Enter snapshot number     = '
      read(*,*) iSnap
      write(*,'(a,$)')' Enter expansion prameter = '
      read(*,*) AEXPN
      write(file1,'(a,i3.3,a)') 'CatshortDistinctA200.',iSnap,'.DAT'
      Open( 2,file=file1)                                 ! output file
      write(file1,'(a,i3.3,a)') 'CatalogDistinctA.',iSnap,'.DAT'
      Open( 3,file=file1,Status='OLD')    !  list of halos

 CALL GetHalos

End PROGRAM HaloProfiles