!---------------------------------
!
!    Find Halo
!
!---------------------------------
Module DATA
Integer(kind=8) :: Nhalos
Integer(kind=4), parameter :: Ntab =60
Real(kind=8), dimension(0:Ntab) :: Dprof, Vrad, Vrad2, Vtan, Vrms
Integer(kind=8), dimension(0:Ntab) :: Nprof
Real(kind=8) :: xc,yc,zc,VXc,Vyc,Vzc, Rad, x0, y0, z0
Real(kind=4), ALLOCATABLE,DIMENSION(:) :: XPAR,YPAR,Zpar,Dens, &
     VX, VY, VZ
Real(kind=4) :: redshift
end Module DATA
!------------------------------------
!
Program FindHalo
  use DATA
  !character :: txt1*10,filename*20

  Rad = 0.050
Call ReadData
Call FindCenter
!xc = 0.; yc = 0. ; zc =0.
Call FindProperties

end Program FindHalo
!------------------------------------
!
SUBROUTINE ReadData
  use DATA
  character :: txt1*10,filename*20

  write(*,'(a)',ADVANCE='NO')'Enter file name: '
  read(*,*)filename
  Open(1,file=TRIM(filename))
  Read(1,*) txt1,redshift
  print *,'Redshift= ',redshift
  dmax = 0.
  x0 = 0.; y0 =0.; z0 =0.
  Nhalos =0
Do
   Read(1,*,iostat=ierr) x
   If(ierr==0)Then
      Nhalos = Nhalos +1
   else
      exit
   end If
end Do
print *,'Number of particles = ',Nhalos
ALLOCATE(Xpar(Nhalos),Ypar(Nhalos),Zpar(Nhalos),Dens(Nhalos))
ALLOCATE(VX(Nhalos),VY(Nhalos),VZ(Nhalos))
rewind(1)
read(1,*) txt1
ii =0
Do i=1,Nhalos
   Read(1,*) Xpar(i),Ypar(i),Zpar(i),vx(i),vy(i),vz(i),Dens(i)
   If(Dens(i)>dmax)Then
      dmax = Dens(i)
      x0 = Xpar(i) ; y0 =Ypar(i); z0 =Zpar(i)
      ii = i
   end If
end Do
close(1)
print '(a,3f9.3,a,f9.2,a,i12)',' Center=',x0,y0,z0, &
     ' Max.density= ',dmax,' particle= ',ii
Xpar(:) = Xpar(:) - x0
Ypar(:) = Ypar(:) - y0
Zpar(:) = Zpar(:) - z0  
end SUBROUTINE ReadData
!------------------------------------
!
SUBROUTINE FindCenter
  use DATA
  Real(kind=8) :: xx, yy, zz, aM

  xc = 0. ; yc =0. ; zc =0.
  R2 = Rad**2
  print *,Rad,R2
  do j=1,20
     Nn  =0
     xx = 0.; yy =0.; zz=0. ; aM = 0.
     Do i=1,Nhalos
        dd = (Xpar(i)-xc)**2 + (Ypar(i)-yc)**2 + (Zpar(i)-zc)**2
        If(dd<R2)Then
           aM = aM + dens(i)
           Nn = Nn + 1
           xx = xx + Xpar(i)*dens(i)
           yy = yy + Ypar(i)*dens(i)
           zz = zz + Zpar(i)*dens(i)
        end If
     end Do
     xc = xx/max(aM,1._8)
     yc = yy/max(aM,1._8)
     zc = zz/max(aM,1._8)
     if(sqrt(R2)>1.5e-3)R2 = R2*0.97
     print '(a,i8,a,3f9.5,a,f9.5)',' particles= ',Nn,' Center= ',xc,yc,zc,' Search Radius= ',(R2)
  end do
  xx =0. ; ii = 0
     Do i=1,Nhalos
        dd = (Xpar(i)-xc)**2 + (Ypar(i)-yc)**2 + (Zpar(i)-zc)**2
        If(dd<R2.and.dens(i)>xx)Then
           xx = dens(i)
           ii = i
        end If
     end Do
     xc = Xpar(ii)
     yc = Ypar(ii)
     zc = Zpar(ii)
     print '(a,f9.2,a,3f9.5,a,f9.5)',' Densest ',dens(ii),' Center= ',xc,yc,zc,' Search Radius= ',(R2)  
end SUBROUTINE FindCenter
!------------------------------------
!
SUBROUTINE FindProperties
  use DATA
  real(kind=8) ::  dx,dy,dz, dvx,dvy,dvz, vv,vvr,vvt 
  
  Rmax  = 0.100
  R2 = Rmax**2
  d  = Rmax/Ntab
  Nn = 0
  VXc = 0.;  VYc = 0.;  VZc = 0   !----- find bulk velocity
  Do i=1,Nhalos
     dd = (Xpar(i)-xc)**2 + (Ypar(i)-yc)**2 + (Zpar(i)-zc)**2
     If(dd<R2)Then
        VXc = VXc + VX(i)
        VYc = VYc + VY(i)
        VZc = VZc + VZ(i)
        Nn  = Nn + 1
     end If
  end Do
  VXc = VXc/max(Nn,1)
  VYc = VYc/max(Nn,1)
  VZc = VZc/max(Nn,1)
  print '(a,i8,a,3f9.3)',' particles= ',Nn,' Velocity Center= ',VXc,VYc,VZc

  
  Dprof(:) = 0. ; Vrad(:)=0. ; Vrad2(:)=0. ; Vtan(:)=0. ; Vrms(:) =0.
  Nprof(:) = 0
     Do i=1,Nhalos
        dd = (Xpar(i)-xc)**2 + (Ypar(i)-yc)**2 + (Zpar(i)-zc)**2
        If(dd<R2)Then
           dd = sqrt(dd+1.e-8)
           ind = INT(dd/d)
           If(ind>Ntab)Stop 'error: radius is to big'
           Nprof(ind) = Nprof(ind) +1
           Dprof(ind) = Dprof(ind) +dd

           dx = Xpar(i)-xc ; dy = Ypar(i)-yc ; dz = Zpar(i)-zc
           dvx = VX(i)-VXc ; dvy = VY(i)-VYc ; dvz = VZ(i)-VZc
           wvx = dvx*dx/dd;  wvy = dvy*dy/dd; wvz = dvz*dz/dd;
           vvr = wvx + wvy + wvz
           vv  = dvx**2 +dvy**2 +dvz**2
           vvt = vv - vvr**2
           !print '(4es12.4,3x,3f9.3,5x,6f9.2)',sqrt(vv),vvr,vvt,dd,dx,dy,dz,dvx,dvy,dvz,wvx,wvy,wvz
           Vrms(ind)  = Vrms(ind)  + vv
           Vrad(ind)  = Vrad(ind)  + vvr
           Vrad2(ind) = Vrad2(ind) + vvr**2
           Vtan(ind)  = Vtan(ind)  + vvt
        end If
     end Do
     open(10,file='HaloProperties.dat',position='APPEND')
     write(10,'(a,f8.3,a,3f9.3)')'Z=',redshift,' Center= ',xc+x0,yc+y0,zc+z0
     write(10,'(a)') ' r_in     <r>    r_out     Np     Overden Ovd(<r_out)  Mass     Vrms    Vradial    Vtangential  beta'
     aNum = 3200.**3/20.**3   !--- number-density of particles
     print *,' number-density= ',aNum
     Ovd = 0.
     aMass  =0.
     aM1 = 2.11e4
     Do ind=0,Ntab
        If(Nprof(ind)/=0)Then
           Dprof(ind) = Dprof(ind)/Nprof(ind)        ! radius
           Vrms(ind)  = sqrt(Vrms(ind)/Nprof(ind))   ! rms velocity
           Vrad(ind)  = Vrad(ind)/Nprof(ind)         ! radial velocity
           Vrad2(ind) = Vrad2(ind)/Nprof(ind)         ! radial velocity
           Vtan(ind)  = sqrt(Vtan(ind)/Nprof(ind))   ! rms tangential velocity
           R1 = ind*d ; R2 = (ind+1)*d
           Volume = 4.188*(R2**3-R1**3)
           Av = aNum*Volume
           beta = 1.- Vtan(ind)**2/2./Vrad2(ind)
           Ovd = Ovd + Nprof(ind)
           aMass = Ovd*aM1
           write(10,'(3f8.4,i8,2f10.3,es12.4,3f9.2,f9.3)')R1,Dprof(ind),R2,Nprof(ind), &
                Nprof(ind)/Av,Ovd/(4.188*R2**3*aNum),aMass,              &
                Vrms(ind),Vrad(ind),Vtan(ind),beta
        end If
     end Do
     
end SUBROUTINE FindProperties