!---------------------------------------------------------------------------------
!     Read HaloTrace files
!     Select Halos with given coords and dump their tracks 
!--------------------------------------------------------------------------------
MODULE   setHalos
PUBLIC      

Integer, parameter  :: Ncats =170          !  Maximum number of  snapshots to analyze

Real, parameter :: Box=250.,   &
                                 Omega0 = 0.27, &  ! Omega_matter
                                 hubble = 0.70
Real    ::  Xcenter,Ycenter,Zcenter,Radius,aSnap

TYPE :: hal
  Real*4         ::   aexpn                      ! expansion parameter
  Integer*4    ::   InCat                        ! pointer to halo in Catshort
  Integer*4    ::   Parent                      ! if not 0, pointer to parent halo in Catshort
  Real         ::   x,y,z,vx,vy,vz        !  coordinates (Mpch units), peculiar velocities(km/s)
  Real         ::  aM,Rm,circV            ! Mass(Msunh), Vir radius(kpch), max circ Velocity
End TYPE hal

TYPE  :: haloEvolution
   Integer :: N                                   ! number of moments for which halo was traced
   Type(hal) :: h(Ncats)
End TYPE haloEvolution

Type(haloEvolution) ::     halo
   
Contains
!---------------------------------------------------------
!                         Read 'iC' halo catalog 
!
      SUBROUTINE ReadHalos(iC)
!--------------------------------------------------------------
Character*120   :: listname,  Txt*5, Txt2*50
Type(hal)                      :: halo2
Real   ::           buff(Ncats)
                              !                                                                                 Open file
      write(listname,'(a,i3.3,a)') 'HaloTree/HaloTrace.',iC,'.dat'
      Open( iC+100,file=TRIM(listname),Status='UNKNOWN',form='unformatted')

   read(iC+100)      ! first line is empty
   nHalos = 0
   Do                                                   ! Loop over all halos in the file
      read(iC+100,iostat=iStat) ih
      If(iStat.ne.0)exit
      Line =0
      Do
         read(iC+100,iostat=iStat) halo2                  ! there is a bug in reading these data
         backspace(iC+100)                                       ! read it twice to get around the bug
         read(iC+100,iostat=iStat) halo2%aexpn   !     read aexpn again
         If(iStat.ne.0)exit
         Line =Line +1
         halo%N = ih
         halo%h(Line) = halo2
         !write(*,'(i9,i8,f9.4,2i9,3f9.4,3f9.2,1pg12.4,0p2f8.2)')ih,nHalos, halo2
      EndDo
      nHalos = nHalos +1

         da =10.
         iClose =1
      Do i =1,Line
         aa = halo%h(i)%aexpn
            If(abs(aa-aSnap).lt.da)then
               da = abs(aa-aSnap) ; iClose = i
            EndIf
      EndDo
         Dist = sqrt((halo%h(iClose)%x -Xcenter)**2 +  &
                             (halo%h(iClose)%y- Ycenter)**2 +  &
                             (halo%h(iClose)%z - Zcenter)**2 )   ! distance to prime halo in kpch
         If(mod(nHalos,10000)==0)&
      write(*,'(a,i8,a,i3,a,i3,4f9.4,1pg12.3)') ' Halo=',nHalos,' iC=',iC,' Lines=',Line, &
                        halo%h(1)%aexpn,halo%h(1)%x,halo%h(1)%y,halo%h(1)%z,halo%h(1)%aM
         If(Dist< Radius)  &
            write(2,'(/a,2i9/(f9.4,2i9,3f9.4,3f9.2,p1g11.3,0p,2f9.2))') ' ----  Halo=',halo%N,Line, &
            (halo%h(j),j=1,Line)

   EndDo

      write(*,'(a,i9,a,i4)') ' Read Halos=',nHalos,' File =',iC
      close(100+iC)
      Return

      End SUBROUTINE ReadHalos

    End module setHalos

!-------------------------------------------------
!            
!            
PROGRAM HaloCrossIdentity

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

   Call CPU_time(t0); call system_clock(iT0,iTc,ib) 
   write(*,*)
         write(*,'(a,$)')' Enter   center (Mpc) and  radius(kpc)  ='
         read(*,*)  Xcenter,Ycenter,Zcenter,Radius
         write(*,'(a,$)')' Enter   Aexpn  ='
         read(*,*)  aSnap

      Open( 2,file='ListSelected.DAT')! updated short list of halos
      !write(2,'(a,T20,a,T38,a,5(2x,a))')'  Halo','Coords(Mpch)',                &
      !                            'Mass(z=0)','MassMax','Vcirc(z=0)','Vcirc(max)','Sub(z=0)','SubEver'


   Call CPU_time(t1); call system_clock(iT2,iTc,ib) ; write(*,*) '          time for Init =',t1-t0,float(iT2-iT0)/iTc

                              ! ---------------------------------------------   Main loop 
Do iC =1, 100
   Call ReadHalos(iC)
   Call CPU_time(t1); call system_clock(iT2,iTc,ib) ; write(*,*) '          time  =',t1-t0,float(iT2-iT0)/iTc
EndDo

stop
end Program HaloCrossIdentity