!---------------------------------------------------------
!               Find spherical voids
!               Find MW candidates in a catshort
Module Structures

Real*4, parameter :: Box = 200., &
                     MassMax = 3.4e12

Real*4 ::     dx,amw
Integer*4 :: Nhalos
Real*4, Allocatable, Dimension(:) :: Xh,Yh,Zh,Vx,Vy,Vz,vrms,Rh,Mh,Vh
Real*4, Dimension(100) :: Xc,Yc,Zc,Rc,Mc,Vc

Contains
!---------------------------------------------------------
  SUBROUTINE GridSearch(N,xl,yl,zl,xr,yr,zr,Radius)

   Radius = 11.
   dx     = 0.5
   xl     = Radius;     yl  = Radius;      zl = Radius
   xr     = Box-Radius; yr  = Box-Radius; zr = Box-Radius

   !xl     = 19.;     yl  = 22.;      zl = 46.
   !xr     = 23.; yr  = 26.; zr = 50
   write(2,*) ' Radius    =',Radius

   N = INT((xr-xl)/dx)

!$OMP parallel do default(shared) private(ii,i,x,y,z,N1,N2,N3,N4,dd,iFlag,aMtot,an)
   Do ii =0,N
      x = xl+ii*dx
      !write(*,*)' x=',x,Radius
      y      = Radius
      Do while(y<yr)
         z      = Radius
         Do while(z<zr)
            N1 = 0; N2 =0; N3 =0; N4 =0
            aMtot = 0. ; iFlag =0
            !write(*,*)x,y,z
            Do i=1,Nhalos
               dd = sqrt((Xh(i)-x)**2 +(Yh(i)-y)**2 +(Zh(i)-z)**2)-Rh(i)/1000.
               If(dd<Radius)Then
                  If(Mh(i)>MassMax)Then
                   iFlag =1 ; exit
                  End If
                  N1 = N1 + 1
                  aMtot = aMtot + Mh(i)
                  If(Vh(i)>250.)N2 =N2+1
                  If(Vh(i)>200.)N3 =N3+1
                  If(Vh(i)>150.)N4 =N4+1
               EndIf
            EndDo
            an = N3/(4.188*Radius**3)
            If(iFlag == 0 .and. an>0.7*amw .and. an<1.3*amw)Then
!$OMP critical
               write(*,'(3f10.4,g12.4,2x,4i5)')x,y,z,aMtot,N2,N3,N4,N1
!$OMP end critical
            End If
            z = z + dx
         EndDo
            y = y + dx
      End Do
   End Do

  end SUBROUTINE GridSearch
!---------------------------------------------------------
  Subroutine MWsearch

   Radius = 10.

   write(2,*) ' Radius    =',Radius

!$OMP parallel do default(shared) private(ic,i,x,y,z,N1,N2,N3,N4,dd,iFlag,aMtot,an)
    Do ic=1,Nhalos
       If(mod(ic,10000) ==0)write(*,*) ' ic =',ic
       If(Vh(ic)>150. .and. Vh(ic)<200.)Then
          x = Xh(ic) ; y = Yh(ic) ; z = Zh(ic) 
            N1 = 0; N2 =0; N3 =0; N4 =0
            aMtot = 0. ; iFlag =0
            !write(*,*)x,y,z
            Do i=1,Nhalos
               dd = sqrt((Xh(i)-x)**2 +(Yh(i)-y)**2 +(Zh(i)-z)**2)-Rh(i)/1000.
               If(dd<Radius)Then
                  If(Mh(i)>MassMax)Then
                   iFlag =1 ; exit
                  End If
                  N1 = N1 + 1
                  aMtot = aMtot + Mh(i)
                  If(Vh(i)>250.)N2 =N2+1
                  If(Vh(i)>200.)N3 =N3+1
                  If(Vh(i)>100.)N4 =N4+1
               EndIf
            EndDo
            an = N4/(4.188*Radius**3)
            If(iFlag == 0 .and. an>0.7*amw .and. an<1.3*amw)Then
!$OMP critical
               write(2,'(3f10.4,g12.4,f9.3,g12.4,2x,4i5)')x,y,z,Mh(ic),Vh(ic),aMtot,N2,N3,N4,N1
               write(3,'(3f9.4,2x,3f9.3,1p,g12.5,0p,4f9.2)')    &
                   Xh(ic),Yh(ic),Zh(ic),vx(ic),vy(ic),vz(ic),Mh(ic),Radius*1000.,vrms(ic),Vh(ic),Rh(ic) 
!$OMP end critical
            End If
      end If
   End Do

  end Subroutine MWsearch
!---------------------------------------------------------
  SUBROUTINE FindSelected

   i = 1
   Do 
      read(3,*,iostat=ierr)Xc(i),Yc(i),Zc(i),v1,v2,v3,Mc(i),Rc(i),vr,Vc(i)
      If(ierr /= 0)Exit
      i = i+1
   Enddo
   Ns = i
   write(2,*) ' Read N centers =',Ns


            N1 = 0; N2 =0; N3 =0; N4 =0; aMtot =0.
!$OMP parallel do default(shared) private(ic,i,x,y,z,dd) REDUCTION(+:N1,N2,N3,N4,aMtot)
    Do ic=1,Nhalos
       If(mod(ic,10000) ==0)write(*,*) ' ic =',ic
          x = Xh(ic) ; y = Yh(ic) ; z = Zh(ic) 
            !write(*,*)x,y,z
            Do i=1,Ns
               dd = sqrt((Xc(i)-x)**2 +(Yc(i)-y)**2 +(Zc(i)-z)**2)
               If(dd<Rc(i)/1000.)Then
                  N1 = N1 + 1
                  aMtot = aMtot + Mh(ic)
                  If(Vh(ic)>250.)N2 =N2+1
                  If(Vh(ic)>200.)N3 =N3+1
                  If(Vh(ic)>100.)N4 =N4+1
!$OMP critical
               write(2,'(3f9.4,2x,3f9.3,1p,g12.5,0p,4f9.2)')    &
                   Xh(ic),Yh(ic),Zh(ic),vx(ic),vy(ic),vz(ic),Mh(ic),vrms(ic),Vh(ic),Rh(ic) 
!$OMP end critical
                  exit
               EndIf
            EndDo
   End Do
write(*,*) ' Nhalos selected =',N1
write(*,*) ' Total   mass    =',aMTot
write(*,*) ' N               =',N2,N3,N4


  End SUBROUTINE FindSelected
End Module Structures
!---------------------------------------------------------
Program FindSpheres
  use Structures
  Character*80 :: String

   Open(1,file='CatshortV.0330.DAT')
   Open(2,file='listOut.dat')
   Open(3,file='CatSelectV.dat')

   Nhalos = 0
   Do i=1,17
      read(1,'(a)') String
      write(*,'(a)') String
      write(2,'(a)') String
     ! write(3,'(a)') String
      read(3,'(a)') String
   EndDo
   Do 
      read(1,*,iostat=ierr)x,y,z
      If(ierr /= 0)Exit
      Nhalos = Nhalos +1
   Enddo
   write(*,*) ' Number of halos= ',Nhalos
   write(2,*) ' Number of halos= ',Nhalos
   Allocate(Xh(Nhalos),Yh(Nhalos),Zh(Nhalos),Rh(Nhalos),Mh(Nhalos),Vh(Nhalos))
   Allocate(Vx(Nhalos),Vy(Nhalos),Vz(Nhalos),vrms(Nhalos))
   rewind(1)
   Do i=1,17
      read(1,'(a)') String
      !write(*,'(a)') String
   EndDo
   N1 = 0 ; N2 =0
   Do i=1,Nhalos
      read(1,*,iostat=ierr)Xh(i),Yh(i),Zh(i),vx(i),vy(i),vz(i),Mh(i),Rh(i),vrms(i),Vh(i)
      If(ierr /= 0)Exit
      If(Vh(i)>150..and.Vh(i)<200.)N1 =N1+1
      If(Vh(i)>100.)N2 =N2+1
   Enddo
   amw = N2/Box**3
   write(2,*) ' Number of MW halos = ',N1
   write(2,*) ' NumDen of    halos = ',amw

   !Call MWsearch
   Call FindSelected

end Program FindSpheres