C ______________________________	
       INCLUDE 'PMbox.h'
       REAL     INPUT
       Character*90 Line
       PI =3.1415926535
C...................................................................
C			Read data and open files
       
      Box =160.      
      BoxV =Box*100.    ! Box size in km/s
      write(*,*) ' Enter: Center(xyz) Radius N_part_limit V_circ_limit '
      read (*,*) Xc,Yc,Zc,Rsphere,N_plim,V_clim
c      Open(1,file='Catshort.999.small.DAT')
      Open(1,file='Catshort.DAT')
      open(18,file='Box.dat')
      CALL ReadHalos(N_plim,V_clim)

c      CALL CONE(ux,uy,uz,ACone,ScaleC)
      CALL SPHERE(Rsphere,Xc,Yc,Zc)
c      CALL DumpParticles(z0,dz0,ScaleC)
      STOP
      END


C---------------------------------------
C                                                       Read Particles 
      SUBROUTINE ReadHalos(N_plim,V_clim)
C---------------------------------------
      INCLUDE 'PMbox.h'
      Character*90 Line
      Lines_header =20
      Do i=1,Lines_header
         Read(1,'(A)') Line
         If(i.le.4)write(18,'(A)') Line
         write(*,'(A)') Line
      ENDDO 
      Nh =0
 10   read(1,*,err=20,end=20) a,x,y,z,vvx,vvy,vvz,am,ar,v1,v2,n1,n2
      If(n2.ge.N_plim.and.v2.ge.V_clim) Then    ! take the halo
         Nh =Nh + 1
         Xm(Nh) =x
         Ym(Nh) =y
         Zm(Nh) =z
         Vx(Nh) =vvx
         Vy(Nh) =vvy
         Vz(Nh) =vvz
         Amc(Nh) =am
         Rmc(Nh) =ar
         Vrm(Nh) =v2
      Endif
      goto 10
 20   write (*,*)'Number of halos read =',Nh,' Np>',N_plim,' Vc>',V_clim
        RETURN
      END
C---------------------------------------
C                              Particles in a sphere
      SUBROUTINE Sphere(Rsphere,Xc,Yc,Zc)
C---------------------------------------
C                 Find points close to the center
      INCLUDE 'PMbox.h'
      write (*,*) ' Center=',Xc,Yc,Zc,' R=',Rsphere

      V  =0.
      V2 =0.
      Vm =0.
      Vm2=0.
      Nn  = 0
      R2 = Rsphere**2
      Ic =0
      Vvx =0.
      Vvy =0.
      Vvz =0.
      Vcx =-215.1
      Vcy =72.3
      Vcz =-5.6
	   DO   IN=1,Nh              ! find central halo
	        dX=(Xm(IN)-Xc)            
	        dY=(Ym(IN)-Yc)
	        dZ=(Zm(IN)-Zc)
           dd =dX**2 +dY**2 +dZ**2
c           write (*,50)IN,Xm(IN),Ym(IN),Zm(IN),dX,dY,dZ
c 50        format(i5,3f8.2,3g11.3)
           If(dd.lt.R2)Then
                Nn =Nn +1
                Vvx =Vvx +Vx(IN)
                Vvy =Vvy +Vy(IN)
                Vvz =Vvz +Vz(IN)
c             If(dd.lt.1.e-4)Then 
c                Ic  =IN
c                Vcx =Vx(Ic)
c                Vcy =Vy(Ic)
c                Vcz =Vz(Ic)
c            Endif
            EndIf
         ENDDO 
         If(Nn.le.1)Then
            write(*,*) '    Too few halos inside the radius: ',Nn
            stop
         Endif
         Vvx =Vvx/Nn
         Vvy =Vvy/Nn
         Vvz =Vvz/Nn
         If(Ic.eq.0)Then     ! no center halo was found
            Vcx =Vvx
            Vcy =Vvy
            Vcz =Vvz
         endif 

         Nc = 0
	   DO   IN=1,Nh            ! Main statistics
	        dX=(Xm(IN)-Xc)            
	        dY=(Ym(IN)-Yc)
	        dZ=(Zm(IN)-Zc)
             dd =dX**2 +dY**2 +dZ**2
             If(dd.lt.R2.and.dd.gt.1.e-4)Then
                d   = max(sqrt(dd+1.e-10),1.e-10)
                Nc =Nc +1
                wx =Vx(IN)-Vcx
                wy =Vy(IN)-Vcy
                wz =Vz(IN)-Vcz
                vr  =(dX*wx+dY*wy+dZ*wz)/d
                ww =sqrt(wx**2+wy**2+wz**2)
c                Xn =dX*(1.+vr/100./d)+Xc
c                Yn =dY*(1.+vr/100./d)+Yc
c                Zn =dZ*(1.+vr/100./d)+Zc
                
                write (*,90)d,ww,vr,100.*d+vr,Amc(IN),Vrm(IN)
 90           format(f8.3,f9.2,2f8.2,
     &              g11.3,f7.2)
                V2 = V2   + ww**2
                Vm2= Vm2  + vr**2
             endif 
      Enddo 

      write (*,*) ' There are ',Nc,' particles inside R=',Rsphere
      write (*,*)
     &               ' Center(Mpc/h) =',Xc,Yc,Zc
      write (*,*)
     &               ' V(km/s)=', Vcx, Vcy, Vcz
      write (*,*) ' Rms Radial velocity=',sqrt(Vm2/Nc),
     &               ' Rms total velocity=',sqrt(V2/Nc)
      Return
      End
C---------------------------------------
C                            write particles in a slice
      SUBROUTINE DumpParticles(z0,dz0,ScaleC)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
C                 Find points close to the center
      Xc =NGRID/2.+1.
      Yc =NGRID/2.+1.
      Zc =NGRID/2.+1.
      Rc =1.
 10   R2 =Rc**2
C				       Loop over particles
      Nc =0
      Vxc =0.
      Vyc =0.
      Vzc =0.
	   DO   IN=1,Np
	        dX=(Xp(IN)-Xc)            !coords
	        dY=(Yp(IN)-Yc)
	        dZ=(Zp(IN)-Zc)
           dd =dX**2 +dY**2 +dZ**2
             If(dd.lt.R2)Then
                Nc =Nc+1
                Vxc =Vxc +Vxp(IN)
                Vyc =Vyc +Vyp(IN)
                Vzc =Vzc +Vzp(IN)
             endif 
      Enddo 
      If(Nc.eq.0)Then
         Rcold =Rc
         Rc       =Rc*1.2
         write (*,*) ' Too small radius=',Rcold,' Increase it to ',Rc
         goto 10
      endif 
      Vxc =Vxc/Nc
      Vyc =Vyc/Nc
      Vzc =Vzc/Nc
      Xc   =ScaleC*(Xc-1.)
      Yc   =ScaleC*(Yc-1.)
      Zc   =ScaleC*(Zc-1.)

      write (*,*) ' There are ',Nc,' particles inside Rc=',Rc*ScaleC,
     &               ' Center(Mpc/h) =',Xc,
     &               Yc,Zc,
     &               ' V(km/s)=', Vxc, Vyc, Vzc
C				       Loop over particles
	   DO   IN=1,Np
	        X=ScaleC*(Xp(IN)-1.)            !coords
	        Y=ScaleC*(Yp(IN)-1.)
	        Z=ScaleC*(Zp(IN)-1.)
             If(Z.gt.z0.and.Z.lt.z0+dz0)Then
	        dX=(X-Xc)            !coords
	        dY=(Y-Yc)
	        dZ=(Z-Zc)
              dd =dX**2 +dY**2 +dZ**2
              d   =max(sqrt(dd+1.e-10),1.e-10)
                ux =Vxp(IN)-Vxc
                uy =Vyp(IN)-Vyc
                uz =Vzp(IN)-Vxc
                vr  =(dX*ux+dY*uy+dZ*uz)/d
                Xn =dX*(1.+vr/100./d)+Xc
                Yn =dY*(1.+vr/100./d)+Yc
                Zn =dZ*(1.+vr/100./d)+Zc

                write (18,100) X,Y,Z,ux,uy,uz,Wp(IN),vr,
     &                 Xn,Yn,Zn
 100            format(3f8.3,3f8.1,g11.3,4f8.3)
             endif 
      Enddo 
      Return
      End
C---------------------------------------
C                          Rotate around Y-axis by angle alpha
      SUBROUTINE RotateYParticles(alpha)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      XN =NGRID+1
c      XC =(NGRID+1.)/2.     ! center of rotation
c      ZC =(NGRID+1.)/2.
      cosA =cos(alpha*3.1415926/180.)
      sinA =sin(alpha*3.1415926/180.)
C				       Loop over particles
	   DO   IN=1,Np
	        X=Xp(IN) -XC            !coords
	        Y=Yp(IN)
	        Z=Zp(IN) -ZC
           Xnew =X*cosA -Z*sinA + XC
           Znew =Z*cosA +X*sinA + ZC
           Vxp(IN) =Vxp(IN)*cosA-Vzp(IN)*sinA
           Vyp(IN) =Vyp(IN) 
           Vzp(IN) =Vzp(IN)*cosA+Vxp(IN)*sinA
           
           X =Xnew
              If(X.ge.XN)X=X-NGRID
              If(X.lt.1.)X=X+NGRID
           Z =Znew
              If(Z.ge.XN)Z=Z-NGRID
              If(Z.lt.1.)Z=Z+NGRID
           Xp(IN) =X
           Yp(IN) =Y
           Zp(IN) =Z

      Enddo 
      Return
      End
C---------------------------------------
C                          Shift particles with periodical conditions
      SUBROUTINE ShiftParticles(xsh,ysh,zsh)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      XN =NGRID+1
C				       Loop over particles
	   DO   IN=1,Np
	        X=Xp(IN)             !coords
	        Y=Yp(IN)
	        Z=Zp(IN)
           X =X+xsh
              If(X.ge.XN)X=X-NGRID
              If(X.lt.1.)X=X+NGRID
           Y =Y+ysh
              If(Y.ge.XN)Y=Y-NGRID
              If(Y.lt.1.)Y=Y+NGRID
           Z =Z+zsh
              If(Z.ge.XN)Z=Z-NGRID
              If(Z.lt.1.)Z=Z+NGRID
           Xp(IN) =X
           Yp(IN) =Y
           Zp(IN) =Z
c           Vxpb(IN) =Vxp(IN)
c           Vypb(IN) =Vyp(IN) 
c           Vzpb(IN) =Vzp(IN)

      Enddo 
      Return
      End
C---------------------------------- Read in variables      
      REAL FUNCTION INPUT(text)
C------------------------------------------------
      Character text*(*)
          write (*,'(A,$)')text
          read (*,*) x
          INPUT =x
      Return
      End
