C ______________________________	Select particles, find density
C		             
C		                     
C          Scales internal PM coordinates and velocities 
C          to  coordinates in Mpc/h and km/s

      INCLUDE 'PMparameters.h'
      INCLUDE 'PMlists.h'
      REAL        AMprevious(Np),AMprevious2(Np)
      REAL        INPUT
      Character   FileASCII*50
C...................................................................
C			Read data and open files
      WRITE (*,'(A,$)') 'Enter File Name for  processed data = '
      READ  (*,'(A)') FileASCII
      OPEN(17,FILE=FileASCII,STATUS='UNKNOWN')
c      OPEN(17,FILE=FileASCII,STATUS='UNKNOWN',FORM='UNFORMATTED')
      CALL RDTAPE

      Cell   = 0.15     ! Cell size in Mpc/h
      write (*,*) ' RDTAPE is done'
      Fraction = 1.
c      Fraction =INPUT(' Enter fraction of particles     =')
      If(Fraction.le.0.)Stop
      Ifraction = INT(1./Fraction+0.1)
      write (*,*) ' Every ',Ifraction,' particle will be selected'
      write (*,*)
     +     ' Enter  Center and Radius XYZ R (Mpc/h)=' 
      read (*,*) Xc,Yc,Zc,Radius
c      Rsearch =0.001
      Rsearch =INPUT(' Enter radius to look for neighbors (Mpc/h)=')
       WRITE (17,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,Xc,Yc,Zc,Radius,Rsearch
100   FORMAT(1X,'Header=>',A45,/
     +            1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3,
     +                 ' Center=',3f8.3,' Rad=',f7.2,' Rsearch=',f7.3)
      Box  =INPUT(' Enter box size in comoving Mpc/h =')

      BoxV =Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
      ScaleC = Box/NGRID         ! scale factor for Coordinates
      xshift =0. 
      yshift =0. 
      zshift =0.
c      CALL GetCentre(ScaleC,ScaleV,Xc,Yc,Zc,
c     &          xshift, yshift ,  zshift  )
c      write (*,*) ' Center=',Xc,Yc,Zc
c      Radius =0.5
      CALL ReadPnts(N,Nsm,ScaleC,ScaleV,xmin,ymin,zmin,
     &             xshift,yshift,zshift,Xc,Yc,Zc,Radius,Ifraction)
      CALL List(Box,N,xmin,ymin,zmin)
      ScaleInc = 4.
c      ScaleInc =INPUT(' Enter factor for larger radius or 1 for none=')
c       WRITE (17) AEXPN,ISTEP,Box,Rsearch,ScaleInc
      write (*,*) ' Center=',Xc,Yc,Zc
      If(ScaleInc.lt.1.1)ScaleInc=1.  ! background scale
      Radius =Rsearch *ScaleInc
      If(ScaleInc.gt.1.1)Then
        write (*,*) ' Search neighbors inside large radius=',Radius
        Do i=1,Nsm
c           if(i.eq.100000)STOP
          if(i/100000*100000.eq.i)write (*,*) '   particle=',i,' of ',N
            X =Xp(i)
            Y =Yp(i)
            Z =Zp(i)
           Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Radius)
           AMprevious(i)=AMnew
c            Rr =Radius *ScaleInc
c           Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Rr)
           AMprevious2(i)=Amnew
        Enddo 
      Else
      write (*,*) ' Center=',Xc,Yc,Zc
         Do i=1,Nsm
           AMprevious(i)=0. 
           AMprevious2(i)=0. 
         EndDo
      Endif 
      Radius =Rsearch
      write (*,*) ' Search neighbors inside small radius=',Radius
      Do i=1,Nsm
          if(i/100000*100000.eq.i)write (*,*) '   particle=',i,' of ',N
         X =Xp(i)
         Y =Yp(i)
         Z =Zp(i)
         Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Radius)
         If(ScaleInc.lt.1.1)AMprevious(i)=0.
         Amextr =Amnew -AMprevious(i)/ScaleInc**3
         If(AMprevious(i).lt.1.05*Amnew)Amextr =0.001
         If(Amextr.lt.0.01)Amextr=0.001
         NpMax = INT(MAX(Amnew,AMprevious(i)))
c         If(iWeight(i).lt.2.)
c     &         write (17) X-Xc,Y-Yc,Z-Zc,Int(Amnew),
c     &                 INT(AMprevious(i)),INT(AMprevious2(i)) !,AMextr,

         If(NpMax.gt.1.and.iWeight(i).lt.2.)
     &      write (17,200) X-Xc,Y-Yc,Z-Zc,Int(Amnew),
     &            INT(AMprevious(i)),AMextr,
     &            INT(iWeight(i))
 200       Format(3F9.4,i6,i7,g11.3,i3)
      Enddo 
      write (*,*) ' Scaled Coordinates were written to:',
     &            FileASCII
      write (*,*) ' Number of particles          =',N

      END

C------------------------------------------
C                                                      Read particles
C
      SUBROUTINE GetCentre(ScaleC,ScaleV,Xc,Yc,Zc,
     &          xshift, yshift ,  zshift  )
C------------------------------------------
       INCLUDE 'PMparameters.h'      
      INCLUDE 'PMlists.h'
      Dimension ListPart(10000)
      Character*90 Line
      Logical Tcen
      Open(18,file='CenterParticles.dat')

      Tcen =.false.
      Box =ScaleC*NGRID
         ifile      =1
      If(Nspecies.eq.0)Then ! old case: 1 complete set
         Npages =NROW          ! Number of pages
         N_in_last=NPAGE       ! Number of particles in last page
      Else                             ! multiple masses
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
C                      --------------------------------------
      If(Tcen)Then                ! assign the center
         Xc = 0.988 + xshift
         Yc = 0.923 + yshift
         Zc = 11.638 + zshift
         Vxcc = -98.643
         Vycc = 110.17
         Vzcc = 35.71
         Radius = 0.002
c                                  find particles close to the center
         Radius2 = Radius**2                       ! square of radius
         write (*,*) '    Center =',Xc,Yc,Zc,' shift=',xshift,
     &                   yshift, zshift,' Radius=',  Radius  
         write (18,*) '    Center =',Xc,Yc,Zc,' shift=',xshift,
     &                   yshift, zshift,' Radius=',  Radius  
         Vxcen =0.
         Vycen =0.
         Vzcen =0.
         Ncen  =0

C				       Loop over particles
       DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            write (*,*)' Read page=',IROW,' file=',ifile,' N=',N
            iL = NPAGE*(IROW-1)
	      CALL GETROW(IROW,ifile)
	      DO IN=1,In_page
	        X  =ScaleC* (XPAR(IN)-1.)
	        Y  =ScaleC* (YPAR(IN)-1.)
	        Z  =ScaleC* (ZPAR(IN)-1.)
                  X =X+xshift 
                  Y =Y+yshift 
                  Z =Z+zshift
                  IF(X.lt.0.)X=X +Box 
                  IF(Y.lt.0.)Y=Y +Box 
                  IF(Z.lt.0.)Z=Z +Box
                  IF(X.GE.Box)X=X -Box 
                  IF(Y.GE.Box)Y=Y -Box 
                  IF(Z.GE.Box)Z=Z -Box
                Vxs=ScaleV* VX(IN)
                Vys=ScaleV* VY(IN)
                Vzs=ScaleV* VZ(IN)
                If(Nspecies.eq.0)Then   ! find weight of particle
                   W =PARTW
                Else
                   Ipart =IN+iL             ! current particle number
                   Do i =1,Nspecies
                      if(Ipart.le.lspecies(i))Then
                         W =wspecies(i)
                         goto 50
                      endif 
                   EndDo 
                EndIf
 50         Ipart  =IN+iL                     ! current particle number
                  iWeigh = INT(W/wspecies(1)+0.1) ! W  ! ! W  ! particles weight   iWeight(Ipart)   
                  If(iWeigh.lt.1.001)Then
                    DD =(X-Xc)**2+(Y-(Yc))**2
     &                     +(Z-(Zc))**2
                    If(DD.lt.Radius2)Then !  inside the limits -- take it
                       Ncen =Ncen +1
                       Vxcen = Vxcen +Vxs
                       Vycen = Vycen +Vys
                       Vzcen = Vzcen +Vzs
                       Vv  = SQRT((Vxs-Vxcc)**2+
     &                       (Vys-Vycc)**2+(Vzs-Vzcc)**2)
        If(Vv.lt.30.)   write (18,'(i8,3f9.4,4f8.2)') Ipart,X,Y,Z,
     &                Vxs-Vxcc,Vys-Vycc,Vzs-Vzcc,Vv
                   EndIf
                 Endif
	      ENDDO
        ENDDO 
 300  write (*,*) '  Number of particle at the center=',Ncen
        write (*,*) '  Radius      =',Radius,' Centre=',Xc,Yc,Zc
        write (*,*) '  Vmean       =',
     &          Vxcen/Ncen,Vycen/Ncen,Vzcen/Ncen
        Close (18)
        Return
c
      ELSE                      ! read the list of particles and find the center
C                      --------------------------------------
c
         read(18,'(A)') Line
         read(18,'(A)') Line
         Ncen =0
         Xc = 0.
         Yc = 0.
         Zc = 0.
 400     read(18,*,err=410,end=410) In,X,Y,Z
               Ncen =Ncen +1
               ListPart(Ncen) =In
         goto 400
c                                           find coords of the particles
 410     Ncurrent =1
         IpCurrent =ListPart(Ncurrent)
C				       Loop over particles
          DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            write (*,*)' Read page=',IROW,' file=',ifile,' Nc=',Ncurrent
            iL = NPAGE*(IROW-1)
	      CALL GETROW(IROW,ifile)
	      DO IN=1,In_page
                Ipart =IN+iL             ! current particle number
                If(Ipart.eq.IpCurrent)Then
	            X  =ScaleC* (XPAR(IN)-1.)
	            Y  =ScaleC* (YPAR(IN)-1.)
	            Z  =ScaleC* (ZPAR(IN)-1.)
                  X =X+xshift 
                  Y =Y+yshift 
                  Z =Z+zshift
                  IF(X.lt.0.)X=X +Box 
                  IF(Y.lt.0.)Y=Y +Box 
                  IF(Z.lt.0.)Z=Z +Box
                  IF(X.GE.Box)X=X -Box 
                  IF(Y.GE.Box)Y=Y -Box 
                  IF(Z.GE.Box)Z=Z -Box
                  Xc =Xc +X
                  Yc =Yc +Y
                  Zc =Zc +Z
                  If(Ncurrent.ge.Ncen) Goto 500
                  Ncurrent =Ncurrent +1
                  IpCurrent =ListPart(Ncurrent)
                EndIf
         EndDo
         EndDo
 500     write (*,*) ' Number of centers found=',Ncurrent
         write(*,*) '    it should be equal to =',Ncen
         Xc =Xc/Ncurrent
         Yc =Yc/Ncurrent
         Zc =Zc/Ncurrent
         write (*,*) ' Center=',Xc,Yc,Zc
      ENDIF   
        Return
        End
C------------------------------------------
C                                                      Read particles
C
      SUBROUTINE ReadPnts(N,Nsm,ScaleC,ScaleV,xmin,ymin,zmin,
     &          xshift,yshift ,zshift ,Xc,Yc,Zc,Radius,Ifraction)
C------------------------------------------
       INCLUDE 'PMparameters.h'      
      INCLUDE 'PMlists.h'
      Integer Xdist(0:25),Ydist(0:25),Zdist(0:25)

      Box =ScaleC*NGRID
      Radius2 =Radius**2
      N =0                                             ! Number of particles
      Nsm = 0
      xmax =-1.e+9
      xmin = 1.e+9
      ymax =-1.e+9
      ymin = 1.e+9
      zmax =-1.e+9
      zmin = 1.e+9
      vmax =-1.e+9
      vmin = 1.e+9
      Do i=0,25
         Xdist(i) =0
         Ydist(i) =0
         Zdist(i) =0
      Enddo 
      Icount =0
      Icount_p =Icount-1
      ifile      =1
      If(Nspecies.eq.0)Then ! old case: 1 complete set
         Npages =NROW          ! Number of pages
         N_in_last=NPAGE       ! Number of particles in last page
      Else                             ! multiple masses
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
      write (*,*) ' N_in_last=',N_in_last
C				       Loop over particles
       DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            write (*,*)' Read page=',IROW,' file=',ifile,' N=',N
            iL = NPAGE*(IROW-1)
            if(Icount_p.eq.Icount.and.Icount_p.gt.100)goto300  ! stop if number does not change
	      CALL GETROW(IROW,ifile)
                  Icount_p=Icount
	      DO IN=1,In_page
	        X  =ScaleC* (XPAR(IN)-1.)
	        Y  =ScaleC* (YPAR(IN)-1.)
	        Z  =ScaleC* (ZPAR(IN)-1.)
C           write (*,'(i8,10g11.4)') IN+iL,
c     &              X,Y,Z,XPAR(IN),YPAR(IN),ZPAR(IN)
                  X =X+xshift 
                  Y =Y+yshift 
                  Z =Z+zshift
                  IF(X.lt.0.)X=X +Box 
                  IF(Y.lt.0.)Y=Y +Box 
                  IF(Z.lt.0.)Z=Z +Box
                  IF(X.GE.Box)X=X -Box 
                  IF(Y.GE.Box)Y=Y -Box 
                  IF(Z.GE.Box)Z=Z -Box
              DD =(X-Xc)**2+(Y-(Yc))**2
     &                     +(Z-(Zc))**2
              If(DD.lt.Radius2)Then !  inside the limits -- take it
                Vxs=ScaleV* VX(IN)
                Vys=ScaleV* VY(IN)
                Vzs=ScaleV* VZ(IN)
                Icount =Icount +1
                vmax =MAX(vmax,Vxs,Vys,Vzs)
                vmin =MIN(vmin,Vxs,Vys,Vzs)
                If(Nspecies.eq.0)Then   ! find weight of particle
                   W =PARTW
                Else
                   Ipart =IN+iL             ! current particle number
                   Do i =1,Nspecies
                      if(Ipart.le.lspecies(i))Then
                         W =wspecies(i)
                         goto 50
                      endif 
                   EndDo 
                EndIf
 50         If(Icount/Ifraction*Ifraction.eq.Icount)Then
                  N   = N +1
                  If(Ipart.le.lspecies(1))Nsm=Nsm +1
                  IF(N.gt.Np)Then    ! too many particles
                           write (*,*) ' Too many particles. Max=',Np
                           STOP
                  Endif 
                  Xp(N)  =X   
                  Yp(N)  =Y   
                  Zp(N)  =Z   
                  Ipart  =IN+iL                     ! current particle number
                  iWeight(N) = INT(W/wspecies(1)+0.1) ! W  ! ! W  ! particles weight   iWeight(Ipart)   
                  If(iWeight(N).lt.1.001)Then
                     xmax =MAX(xmax,Xp(N))
                     xmin  =MIN(xmin,Xp(N))
                     ymax =MAX(ymax,Yp(N))
                     ymin  =MIN(ymin,Yp(N))
                     zmax =MAX(zmax,Zp(N))
                     zmin  =MIN(zmin,Zp(N))
                     ii =Min(INT(X/Box*25.),25)
                     jj =Min(INT(Y/Box*25.),25)
                     kk =Min(INT(Z/Box*25.),25)
                     Xdist(ii) =Xdist(ii) +1
                     Ydist(jj) =Ydist(jj) +1
                     Zdist(kk) =Zdist(kk) +1
                   EndIf
                 Endif
              Endif      ! end test for center
	      ENDDO
        ENDDO 
 300    write (*,*) ' Number of particles selected:=',N,' Small=',Nsm
        write (*,*) ' Limits                              =',xmin,xmax       
        write (*,*) '                                     =',ymin,ymax
        write (*,*) '                                     =',zmin,zmax
        write (*,*) '  Radius      =',Radius,' Centre=',Xc,Yc,Zc
        write (*,'(i3,2x,i9,2x,i9,2x,i9)')
     &             (i,Xdist(i),Ydist(i),Zdist(i),i=0,25)
        Return
        End
C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE List(Box,N,xmin,ymin,zmin)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMlists.h'
             Do i=1,N
                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,N
         i=INT((Xp(jp)-xmin)/Cell)
         j=INT((Yp(jp)-ymin)/Cell)
         k=INT((Zp(jp)-zmin)/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; ww-its weight
C
      SUBROUTINE Neib(Amnew,xmin,ymin,zmin,Xc,Yc,Zc,Radius)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMlists.h'
C                                       Initiate counters
  	  Amnew =0.
        Radius2=Radius**2
c                      limits for Label
         i1=INT((Xc-xmin-Radius)/Cell)
         j1=INT((Yc-ymin-Radius)/Cell)
         k1=INT((Zc-zmin-Radius)/Cell)
            i1=MIN(MAX(Nm,i1),Nb)
            j1=MIN(MAX(Nm,j1),Nb)
            k1=MIN(MAX(Nm,k1),Nb)
         i2=INT((Xc-xmin+Radius)/Cell)
         j2=INT((Yc-ymin+Radius)/Cell)
         k2=INT((Zc-zmin+Radius)/Cell)
            i2=MIN(MAX(Nm,i2),Nb)
            j2=MIN(MAX(Nm,j2),Nb)
            k2=MIN(MAX(Nm,k2),Nb)
C                                        Look for neibhours
c            nn =0
         Do k=k1,k2
         Do j=j1,j2
         Do i=i1,i2
            jp =Label(i,j,k)
 10         If(jp.ne.0)Then
               dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
               If(dd.lt.Radius2)Then
c                  nn =nn +1
	             Amnew=Amnew + iWeight(jp)
               EndIf  
                  jp =Lst(jp)
               GoTo10
            EndIf
         EndDo
         EndDo
         EndDo

c         write (*,100) nn,i1,i2,j1,j2,k1,k2,Xc,Yc,Zc
c 100     format(' Nn=',i6,3(2i4,2x),3f8.3)
      Return
      End
