C-------------------------------------------------
c
C     
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'
      Real             INPUT

      Np_min  =INPUT('Enter Minimum number of particles in a halo => ')
      V_c_min=INPUT('Enter Minimum circular velocity             => ')  
      Radius  =INPUT('Enter comoving radius(Mpc/h) for spheres    => ')
      Nsph     =INPUT('Enter number of  random spheres           => ')
      Fract   =INPUT('Enter fraction of DM particles (1/4,1/2,1)=> ')
      Fract   =Min(Max(Fract,0.),1.0)                             
      Box     =INPUT('Enter Comoving Box size(Mpc/h)            => ')
      write (*,*)
C     ------------------------------------------ Open files
      Open( 2,file='Vresults.DAT')   ! short list of halos
      Open( 3,file='Catshort.DAT')
C     ------------------------------------------ Read data
      CALL RDTAPE
            write (*,*) ' RDTAPE is done'
            write (*,*) ' Nnu=Number of Species=',Nspecies
         WRITE (2,100) HEADER,
     +               AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +               EKIN,NROW,NGRID,Om0,Oml0,hubble
100      FORMAT(1X,'Header=>',A45,/
     +         ' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +         ' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,
     +         ' Nrow=',I4,' Ngrid=',I4,/' Omega_0=',f6.3,
     +         ' Omega_L=',f6.3,' h=',f6.3)
C     ------------------------------------------ Constants ---------
           hsmall= hubble            ! Hubble/100.
           Box   = Box /hsmall       ! Scale to real Mpc          
           Radius= Radius/hsmall  ! Comoving Search radius in real Mpc
           Xscale= Box/Ngrid        ! Scale for comoving coordinates
           Vscale= 100.*hsmall*Xscale/AEXPN ! Scale for velocities
           Dscale= 2.746e+11*hsmall**2*(Box/NROW)**3 ! mass scale
C                            Dscale is a factor used to convert particles to mass 
C                                       =  mass of a particle for Omega=1 model
C                            Xscale is a factor to convert PM coordinates into Mpc
C                                       = Box / Ngrid = length of a PM cell in Mpc
C                           Vscale is a factor to convert PM pomenta to V (km/sec)
C                                       = the Hubble velocity at the distance of 1 PM cell
C                              the factor AEXPN is needed to go from momenta to pec.velocity
           write (*,*) ' Vscale=',Vscale, ' X=',Xscale,' Box(Mpc)=',Box
           Radius=Radius/Xscale    ! comoving radius in cell units
           RadSr1 = Radius**2
           Cell  =2.0              ! Cell-size for Linker List in grid units
           Box   =NGRID            ! all internal variables are in Grid units
           imin  =INT(-Radius/Cell)! Make recommendations. Code works 
           imax  =INT((Box+Radius)/Cell) ! if the limits are not optimal
                   If(imin.lt.Nm)Write(*,*)
     .               ' Wrong low limit:',Nm,' better be',imin
                   If(imax.gt.Nb)Write(*,*)
     .               ' Wrong upper limit:',Nb,' better be',imax
           write (*,*) 
     .        ' Box(cells)=',Box,' Cell=',Cell,' Rad(cells)=',Radius
c         DtoNscale = Om0*Dscale/Fract ! scale dark matter particles
C 
C      CALL ReadCent(Nsph,Ncentr)                     ! random centers          
      CALL ReadMax(N,Xscale,Np_min,V_c_min)            ! read halos
      Open( 3,file='Catshort.DAT')
      Ncentr = N
          Do i=1,N
                  Rmc(i) =Radius
                  Amc(i) =1.
                  Xm(i)   =Xp(i)
                  Ym(i)   =Yp(i)
                  Zm(i)   =Zp(i)
         EndDo
          
      CALL ReadPnt(N,Fract,Vscale)                         ! read particles
      Ndm  = N
c----------------------------------   DM particles
      CALL Points(Box,N,Ncp)            ! Make more points periodically
          write (*,*) ' Nparticles in periodic set=',Ncp,
     .              '( Maximim =',Np,')'
          write (2,*) ' Ndm in periodic set=',Ncp,
     .              '( Maximum =',Np,')'
      Call List(Box,Ncp)                ! make linker List, Label

          write (*,*) ' Nparticles=',N,' Ncenters=',Ncentr
      CALL SPairs(Box,N,Ncentr,1)  ! Accumulate  statistics for particles 
 
C---------------------------------- Halos
      CALL ReadMax(N,Xscale,Np_min,V_c_min)            ! read halos
      Nhalos = N
      CALL Points(Box,N,Ncp)            ! Make more points periodically

          write (*,*) ' N_halos=',Nhalos,'  in periodic set=',Ncp,
     .              '( Maximum =',Np,')',' N_dm min=',Np_min,
     .              ' V_circ_min=',V_c_min
          write (2,*) ' N_halos=',Nhalos,'  in periodic set=',Ncp,
     .              '( Maximum =',Np,')'
          write (2,*) ' N_dm min=',Np_min,
     .              ' V_circ_min=',V_c_min
      Call List(Box,Ncp)                ! make linker List, Label

          write (*,*) ' Nparticles=',N,' Ncenters=',Ncentr
      CALL SPairs(Box,N,Ncentr,2)  ! Accumulate  statistics for halos
C------------------------------- Print Results
      write ( 2,13)
 13   Format(7x,'N_dm Contrast N_halos V_dm rmsV_dm',
     &                     ' V_halos rmsV_halos')
      AMmean =Ndm/FLOAT(NGRID**3)*4.188*Radius**3 
      Do i=1,Ncentr          ! Scale all to physical units
         Xc    =Xm(i)  *Xscale
         Yc    =Ym(i)  *Xscale
         Zc    =Zm(i)  *Xscale
         aM   =Nbc(i,1)
         Contrast = aM/AMmean
c         If(Contrast.ge.1.)Then
         VelocityD  = sqrt(max(1.e-10,
     &        Wxc(i,1)**2 +Wyc(i,1)**2 +Wzc(i,1)**2))/
     &         max(Nbc(i,1),1)
         VelocityH  = sqrt(max(1.e-10,
     &        Wxc(i,2)**2 +Wyc(i,2)**2 +Wzc(i,2)**2))/
     &         max(Nbc(i,2),1)
         rmsD=sqrt(max(1.e-10,Vrmc(i,1)/max(Nbc(i,1),1) -
     &                  VelocityD**2))   
         rmsH=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) -
     &                  VelocityD**2))
         rmsH2=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) -
     &                  VelocityH**2))
         rmsH3=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) ))
c         write (*,200) i,Nbc(i,1),Contrast ,Nbc(i,2),
c     &                      VelocityD,rmsD, VelocityH, rmsH,rmsH2,rmsH3
         write (2,200) i,Nbc(i,1),Contrast ,Nbc(i,2),
     &                      VelocityD,rmsD, VelocityH, rmsH,rmsH2,rmsH3,
     &                      Xc,Yc,Zc
 200     format(i5,i7,f8.3,i5,2(f8.2,f8.2,2x),2f8.2,3f7.3)
c         Endif 
        Enddo 
      Stop
      End
C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE List(Box,Ncp)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'
             Do i=1,Ncp
                Lst(i)=-1
             EndDo

             Do k=Nm,Nb
             Do j=Nm,Nb
             Do i=Nm,Nb
                Label(i,j,k)=0
             EndDo
             EndDo
             EndDo  
      Do jp=1,Ncp
         i=INT(Xp(jp)/Cell)
         j=INT(Yp(jp)/Cell)
         k=INT(Zp(jp)/Cell)
         i=MIN(MAX(Nm,i),Nb)
         j=MIN(MAX(Nm,j),Nb)
         k=MIN(MAX(Nm,k),Nb)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      Return
      End
C--------------------------------------------------------------
C                                Find all neibhours for a center
C                                Xc,Yc,Zc - center; 
C                                ic = center, js = 1 (dm), =2 (halos)    
      SUBROUTINE SNeib(ic,js)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'

         Radmax =Radius

c                      limits for Label
         i1=INT((Xc-Radmax)/Cell)
         j1=INT((Yc-Radmax)/Cell)
         k1=INT((Zc-Radmax)/Cell)
            i1=MIN(MAX(Nm,i1),Nb)
            j1=MIN(MAX(Nm,j1),Nb)
            k1=MIN(MAX(Nm,k1),Nb)
         i2=INT((Xc+Radmax)/Cell)
         j2=INT((Yc+Radmax)/Cell)
         k2=INT((Zc+Radmax)/Cell)
            i2=MIN(MAX(Nm,i2),Nb)
            j2=MIN(MAX(Nm,j2),Nb)
            k2=MIN(MAX(Nm,k2),Nb)
C                                        Look for neibhours
         Do k=k1,k2
         Do j=j1,j2
         Do i=i1,i2
            jp =Label(i,j,k)   !  Dark matter
 10         If(jp.ne.0)Then
               dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
                  If(dd.lt.RadSr1)Then
                    Nbc(ic,js)= Nbc(ic,js)+ 1
                    Rrc(ic,js)= Rrc(ic,js)+ sqrt(dd)
                    Vrmc(ic,js) =Vrmc(ic,js) +
     &                                VXp(jp)**2+VYp(jp)**2+VZp(jp)**2
                    Wxc(ic,js) =Wxc(ic,js) +VXp(jp)  
                    Wyc(ic,js) =Wyc(ic,js) +VYp(jp)
                    Wzc(ic,js) =Wzc(ic,js) +VZp(jp)
                  EndIf
               jp =Lst(jp)
               GoTo10
            EndIf
         EndDo
         EndDo
         EndDo
      Return
      End
C--------------------------------------------------------------
C                       Add points around the Box to take into account periodical conditions
C  It replicates points periodically and keeps those which are at  distance
C  less than Radius from the Box. 
C                        N     = number of points inside the Box
C                        Ncp = total number of points
      SUBROUTINE Points(Box,N,Ncp)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'
      Logical    Inside
         B1  =-Radius
         B2  =Box+Radius
         B3  =Box-Radius
         Ncp =N
      Do i=1,N               ! Add dark matter particles
         x0 =xp(i)
         y0 =yp(i)
         z0 =zp(i)
         Inside =            (x0.gt.Radius).and.(x0.lt.B3)
         Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3))
         Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3))
         If(.not.Inside)THEN         
         ux =VXp(i)
         uy =VYp(i)
         uz =VZp(i)
            Do k1 = -1,1
               zr =z0+k1*Box
               If(zr.GT.B1.AND.zr.LT.B2)Then
               kk =k1**2
               Do j1 = -1,1
                  yr =y0+j1*Box
                  If(yr.GT.B1.AND.yr.LT.B2)Then
                     jj = kk +j1**2
                     Do i1 = -1,1
                        xr =x0+i1*Box
                        ii =jj +i1**2
                        If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)
     .                     Then
                             Ncp =Ncp +1
                             If(Ncp.le.Np)Then
                               xp(Ncp)  =xr
                               yp(Ncp)  =yr
                               zp(Ncp)  =zr
                               VXp(Ncp)  =ux
                               VYp(Ncp)  =uy
                               VZp(Ncp)  =uz
                             EndIf
                        EndIf
                     EndDo
                  EndIf
               EndDo
               EndIf
            EndDo
         EndIf
      EndDo
      If(Ncp.Gt.Np)Then
         write(*,*)' Too many points:',Ncp,' Maximum is =', Np
         STOP
      EndIf
      Return
      End

C-------------------------------- Update Statistics of pairs 
C                   
      SUBROUTINE SPairs(Box,N,Ncentr,js)
C---------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'

      If(Ncentr.gt.Ncc)Then
         write (*,*) ' Number of centers is too large:'
         write (*,*) '         Ncentr =',Ncentr,' Max=Ncc=',Ncc,' STOP'
         STOP
      EndIf
         Do i=1,Ncentr
                  Wxc(i,js) = 0.
                  Wyc(i,js) = 0.
                  Wzc(i,js) = 0.
                  Nbc(i,js) = 0
                  Rrc(i,js) = 0.
                  Vrmc(i,js)= 0.
         EndDo
      Do ic=1,Ncentr
         Xc =Xm(ic)
         Yc =Ym(ic)
         Zc =Zm(ic)
         Call SNeib(ic,js)
      EndDo
      Return
      End

C---------------------------------------------------------
C                                                  Select random particles
c                                                   as initial seeds of centers
      SUBROUTINE ReadCent(N,Ncentr)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'
      
      Nseed =120731
          Do i=1,N
                  Rmc(i) =Radius
                  Amc(i) =1.
                  Xm(i)   =NGRID*RANDd(Nseed)
                  Ym(i)   =NGRID*RANDd(Nseed)
                  Zm(i)   =NGRID*RANDd(Nseed)
         EndDo
         Ncentr =N
      Return
      End

C---------------------------------------------------------
C                                        Read list of maxima
c                                     Xnew =X-1: coord. now are 0-NGRID       
      SUBROUTINE ReadMax(Ncentr,Xscale,Np_min,V_c_min)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'
      Character*79 Text
      Do i=1,16
         read(3,'(A)') Text
         write(*,'(A)') Text
      enddo 
       Ncentr =0
      xxm =1000.
      xxn =-1000.
 10   read (3,*,err=5,end=5) X0, Y0, Z0, VvX,VvY,VvZ,
     &                                     amh,rh,wvel,Vcirc,ihh
           If(ihh.ge.Np_min.and.Vcirc.ge.V_c_min)Then
            Ncentr=Ncentr+1
            If(Ncentr.le.Ncc)Then
            Xp(Ncentr)  =X0/Xscale/hubble
            Yp(Ncentr)  =Y0/Xscale/hubble
            Zp(Ncentr)  =Z0/Xscale/hubble
            VXp(Ncentr)=VvX
            VYp(Ncentr)=VvY
            VZp(Ncentr)=VvZ
            xxm    =min(xxm,Xp(Ncentr))
            xxn     =max(xxn,Xp(Ncentr))
            EndIf
            EndIf
            If(Ncentr.lt.Ncc)GoTo 10
 5    write(*,*) ' Read Halos: X{Min,Max}=',xxm,xxn,' N=',Ncentr
      CLOSE (3)
      If(Ncentr.GT.Ncc)Then
         write (*,*) ' Too many particles: Max=',Ncc
         Stop
      EndIf
      Return
      End
C---------------------------------------------------------
C                                        Read particles
      SUBROUTINE ReadPnt(N,Fract,Vscale)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMvlists.h'

      N =0
      Iz=1
      Iy=1
      Ff=1.
      If(Fract.le.0.666)Then
         Iz=2
         Ff=0.5
         If(Fract.le.0.333)Then
            Iy=2
            Ff=0.25
         Endif
      EndIf
      Fract =Ff
C				       Loop over particles
      DO 10 IROW=1,NROW,Iz
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,NPAGE,Iy
            N   = N +1
            Xp(N)  =XPAR(IN)-1.
            Yp(N)  =YPAR(IN)-1.
            Zp(N)  =ZPAR(IN)-1.
            VXp(N)=VX(IN)*Vscale
            VYp(N)=VY(IN)*Vscale
            VZp(N)=VZ(IN)*Vscale
         ENDDO
 10   CONTINUE
      Return
      End

