C ______________________________	  
C		                    
C	     NROW = number of particles per row
C	     NGRID= number of cells on a side (normal grid)
       INCLUDE 'PMparameters.h'
       INCLUDE 'PMfields.h'
       REAL     INPUT
       PI =3.1415926535
C...................................................................
C			Read data and open files
      CALL RDTAPE

      write (*,*) ' RDTAPE is done'
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,/
     +            1x,' Number of part.1specie=',i8,' Nspecies=',i3)
c      Box =INPUT(' Enter box size in Mpc/h   =')
      Box =160.
c      alpha =INPUT(' Enter angle of rotation   =')
      
      BoxV =Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
      ScaleC = Box/NGRID         ! scale factor for Coordinates

C
      CALL ReadParticles(ScaleV)
      xsh = 3.2/ScaleC                   ! shift coordinates
      ysh = -12.3/ScaleC
      zsh = 0.8/ScaleC
c      Xc = -0.72                            !center of the VIRGO
c      Yc =14.43
c      Zc =0.83
      Xc = 0.                           !center of the cone
      Yc = 0.
      Zc = 0.
      Rsphere =5.
      Rc =40.
      ACone=20.
      ux =-1.                           ! direction of the cone
      uy =13.68
      uz =1.
c      ux =INPUT(' Enter x-direction   =')

      CALL ShiftParticles(xsh,ysh,zsh)
c      CALL RotateYParticles(alpha)

      open(18,file='Cone.dat')
      WRITE (18,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1),
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,lspecies(1),Nspecies
      WRITE (18,*) ' Angle=',alpha,' Cone=',Acone,'deg R=',
     +            Rc,' direction=',ux,uy,uz
      ACone=ACone*PI/180.
      Xc =(Xc/Box +0.5)*NGRID+0.5
      Yc =(Yc/Box +0.5)*NGRID+0.5
      Zc =(Zc/Box +0.5)*NGRID+0.5
      Rc =Rc/Box*NGRID
      CALL CONE(ux,uy,uz,ACone,ScaleC)
c      CALL SPHERE(Rsphere,ScaleC)
c      CALL DumpParticles(z0,dz0,ScaleC)
      STOP
      END


C---------------------------------------
C                                                       Read Particles 
      SUBROUTINE ReadParticles(ScaleV)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)

      W     = (FLOAT(NGRID)/FLOAT(NROW))**3
      PARTW = W  
      WPAR = wspecies(1)
      N =0
      xmin =1.e+5
      xmax =-1.e+5
      ymin =1.e+5
      ymax =-1.e+5
      zmin =1.e+5
      zmax =-1.e+5
      vxmin =1.e+15
      vxmax =-1.e+15
      vymin =1.e+15
      vymax =-1.e+15
      vzmin =1.e+15
      vzmax =-1.e+15

C				       Loop over particles
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      write (*,*) Nspecies,
     + lspecies(Nspecies),wspecies(Nspecies)
      Nsp_current =1
      WPAR     = wspecies(1)
      Lspec_next=lspecies(1)+1
      DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            If(IROW/50*50.eq.IROW)
     +        write (*,*) ' Read page=',IROW,' N=',N,
     +             ' Npages=',Npages,In_page
            iL = NPAGE*(IROW-1)
C				       Loop over particles
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,In_page
            N   = N +1
            If(N.ge.Lspec_next)Then   ! next specie
               write (*,*) ' Another specie. Current=',Nsp_current
               write (*,*) '     Particle=',N,WPAR,
     &                         ' lspecies=',lspecies(Nsp_current),
     &                         ' next=',lspecies(Nsp_current+1)
               Nsp_current=Nsp_current+1
               Lspec_next =lspecies(Nsp_current)+1
               WPAR     = wspecies(Nsp_current)
            EndIf
	        X=XPAR(IN)
	        Y=YPAR(IN)
	        Z=ZPAR(IN)
           If(N.gt.Nmaxp)Then
              write (*,*) ' Error: not enough space for particles:',N
              stop
           Endif 
           Xp(N) =X                             ! coordinates in 1-NGRID
           Yp(N) =Y
           Zp(N) =Z
           Vxp(N)=ScaleV* VX(IN)    ! velocities in km/s
           Vyp(N)=ScaleV* VY(IN)
           Vzp(N)=ScaleV* VZ(IN)
           Wp(N) =WPAR                   ! particle weight

           xmin =MIN(xmin,X)
           ymin =MIN(ymin,Y)
           zmin =MIN(zmin,Z)
           xmax =MAX(xmax,X)
           ymax =MAX(ymax,Y)
           zmax =MAX(zmax,Z)
           vxmin =MIN(vxmin,Vxp(N))
           vymin =MIN(vymin,Vyp(N))
           vzmin =MIN(vzmin,Vzp(N))
           vxmax =MAX(vxmax,Vxp(N))
           vymax =MAX(vymax,Vyp(N))
           vzmax =MAX(vzmax,Vzp(N))
	        I=INT(X)
	        J=INT(Y)
	        K=INT(Z)
            If(I.le.0)write (*,*) ' X:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(J.le.0)write (*,*) ' Y:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(K.le.0)write (*,*) ' Z:',X,Y,Z,I,' Irow=',IROW,IN,ifile

	      ENDDO
        ENDDO 

        Np = N
      write (*,*) ' Number of particles read=', N
      write (*,*) ' X:',xmin,xmax,' V=',vxmin,vxmax
      write (*,*) ' Y:',ymin,ymax,' V=',vymin,vymax
      write (*,*) ' Z:',zmin,zmax,' V=',vzmin,vzmax
        RETURN
      END
C---------------------------------------
C                              Particles in a sphere
      SUBROUTINE Sphere(Rsphere,ScaleC)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'

      Dimension Rline(0:100),Vline(0:100),Mline(0:100),
     &                V2line(0:100),Wline(0:100)
      Dimension Xq(250000),Yq(250000),Zq(250000),
     &               VXq(250000),VYq(250000),VZq(250000)
C                 Find points close to the center
      write (*,*) ' Center=',Xc,Yc,Zc,' R=',Rsphere
      write (*,*) ' ScaleC=',ScaleC
      dRline =0.15
      Do i=0,100
         Rline(i) =0.
         Vline(i) =0.
         Mline(i) =0
         V2line(i) =0.
         Wline(i) =0.
      Enddo 
      R =Rsphere /ScaleC/10.    ! Find velocity of the center
 10   R2 =R**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
                 dX =dX*ScaleC
                 dY =dY*ScaleC
                 dZ =dZ*ScaleC
c           write (*,90)IN,dX,dY,dZ,Vxp(IN),Vyp(IN),Vzp(IN)
c 90        format(' i=',i8,' x=',3f8.2,' V=',3f9.2)
                Nc =Nc+1
                Vxc =Vxc +Vxp(IN)
                Vyc =Vyc +Vyp(IN)
                Vzc =Vzc +Vzp(IN)
             endif 
      Enddo 
      If(Nc.lt.5)Then
         Rold =R
         R       =R*1.2
         write (*,*) ' Too small radius=',Rold,' Increase it to ',R
         goto 10
      endif 
      Vxc =Vxc/Nc                  ! This is the center
      Vyc =Vyc/Nc
      Vzc =Vzc/Nc
      Xc0   =ScaleC*(Xc-1.)
      Yc0   =ScaleC*(Yc-1.)
      Zc0   =ScaleC*(Zc-1.)

      write (*,*) ' There are ',Nc,' particles inside R=',R*ScaleC,
     &               ' Center(Mpc/h) =',Xc0,
     &               Yc0,Zc0,
     &               ' V(km/s)=', Vxc, Vyc, Vzc
C				    Find radial velocities of particles in the Sphere
      Nc =0
      Rsph =Rsphere/ScaleC
      R3    =Rsph**2 
	   DO   IN=1,Np             ! Loop over particles
	        X=Xp(IN)            !coords
	        Y=Yp(IN)
	        Z=Zp(IN)
	        dX=(X-Xc)            !coords
	        dY=(Y-Yc)
	        dZ=(Z-Zc)
              dd =dX**2 +dY**2 +dZ**2
              If(dd.lt.R3)Then
                 d     = max(sqrt(dd+1.d-10),1.d-10)
                    Nc =Nc +1
                wx =Vxp(IN)-Vxc
                wy =Vyp(IN)-Vyc
                wz =Vzp(IN)-Vzc
                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
                d =d*ScaleC
                dX =dX*ScaleC
                 dY =dY*ScaleC
                 dZ =dZ*ScaleC
              write (18,100) d,vr,ww,dX,dY,dZ,wx,wy,wz,Wp(IN)
 100            format(f8.2,2f9.2,3f8.3,3f8.1,g11.3,4f8.3)
                iL =max(min(int(d/dRline),100),0)
                Rline(iL) = Rline(iL) +d*Wp(IN)
                Wline(iL) = Wline(iL) +Wp(IN)
                Mline(iL) = Mline(iL) +1
                Vline(iL) = Vline(iL) +vr*Wp(IN)
                V2line(iL) = V2line(iL) +vr**2*Wp(IN)
c                write (*,*) d,dRline,
                If(Nc.gt.250000)Then
                   write (*,*) ' Not enough space'
                   stop
                endif 
                Xq(Nc) =dX
                Yq(Nc) =dY
                Zq(Nc) =dZ
                VXq(Nc) =wx
                VYq(Nc) =wy
                VZq(Nc) =wz
             endif 
      Enddo 
      write (*,*) ' N_inside_sphere=',Nc,' R=',Rsphere
      open(30,file='Sphere.dat')
      do i=0,99
         If(Mline(i).gt.3)Then
            Rline(i) =Rline(i)/Wline(i)
            Vline(i) = Vline(i)/Wline(i)
            V2line(i) =sqrt(V2line(i)/Wline(i)-Vline(i)**2)
            write(30,80) i,(i+0.5)*dRline,Rline(i),Vline(i), 
     &                        V2line(i),Mline(i),Wline(i)
 80         format(i3,f8.3,f8.3,2f9.2,i6,g11.3)
         EndIf
      Enddo 
      Rsmooth =0.15
      R2 =Rsmooth**2
      Do i=1,Nc
         x =Xq(i)
         y =Yq(i)
         z =Zq(i)
         Nn =0
         wx =0.
         wy =0.
         wz =0.
         Do j=1,Nc
            dd =(x-Xq(j))**2+(y-Yq(j))**2+(z-Zq(j))**2
            if(dd.lt.R2)Then
               Nn =Nn +1
               wx =wx+Vxq(j)
               wy =wy+Vyq(j)
               wz =wz+Vzq(j)
            endif 
         Enddo 
         If(Nn.eq.0)write(*,*) ' Error Nn'
               wx =wx/Nn
               wy =wy/Nn
               wz =wz/Nn
               d   =sqrt(x**2+y**2+z**2)
                vr  =(x*wx+y*wy+z*wz)/d
                ww =sqrt(wx**2+wy**2+wz**2)
                vra  =(x*Vxq(i)+y*Vyq(i)+z*Vzq(i))/d
                wwa=sqrt(Vxq(i)**2+Vyq(i)**2+Vzq(i)**2)
              write (30,110) d,vr,vra,ww,wwa,x,y,z,wx,wy,wz,Nn
 110            format(f8.2,4f9.2,3f8.3,3f8.1,i4)
          Enddo 
      Return
      End
C---------------------------------------
C                              Particles in a cone
      SUBROUTINE CONE(ux,uy,uz,ACone,ScaleC)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      Dimension Rline(0:100),Vline(0:100),Mline(0:100),
     &                V2line(0:100),Wline(0:100)
C                 Find points close to the center
      write (*,*) ' Center=',Xc,Yc,Zc,' R=',Rc
      write (*,*) ' ScaleC=',ScaleC
      dRline =1.
      Do i=0,100
         Rline(i) =0.
         Vline(i) =0.
         Mline(i) =0
         V2line(i) =0.
         Wline(i) =0.
      Enddo 
      R =1.5
 10   R2 =R**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
                 dX =dX*ScaleC
                 dY =dY*ScaleC
                 dZ =dZ*ScaleC
c           write (*,90)IN,dX,dY,dZ,Vxp(IN),Vyp(IN),Vzp(IN)
 90        format(' i=',i8,' x=',3f8.2,' V=',3f9.2)
                Nc =Nc+1
                Vxc =Vxc +Vxp(IN)
                Vyc =Vyc +Vyp(IN)
                Vzc =Vzc +Vzp(IN)
             endif 
      Enddo 
      If(Nc.lt.5)Then
         Rold =R
         R       =R*1.2
         write (*,*) ' Too small radius=',Rold,' Increase it to ',R
         goto 10
      endif 
      Vxc =Vxc/Nc
      Vyc =Vyc/Nc
      Vzc =Vzc/Nc
      Xc0   =ScaleC*(Xc-1.)
      Yc0   =ScaleC*(Yc-1.)
      Zc0   =ScaleC*(Zc-1.)

      write (*,*) ' There are ',Nc,' particles inside R=',R*ScaleC,
     &               ' Center(Mpc/h) =',Xc0,
     &               Yc0,Zc0,
     &               ' V(km/s)=', Vxc, Vyc, Vzc
C				    Find radial velocities of particles in the Cone
      Nc =0
      R2 =Rc**2
      u = sqrt(ux**2 +uy**2 +uz**2)
      ux =ux/u                             ! unit vector for cone direction
      uy =uy/u
      uz =uz/u
      CosCone =COS(ACone)
      Rcyl =1.5/ScaleC
      R3    =Rcyl**2 
	   DO   IN=1,Np             ! Loop over particles
	        X=Xp(IN)            !coords
	        Y=Yp(IN)
	        Z=Zp(IN)
	        dX=(X-Xc)            !coords
	        dY=(Y-Yc)
	        dZ=(Z-Zc)
              dd =dX**2 +dY**2 +dZ**2
              If(dd.lt.R2.and.Wp(IN).lt.5.)Then
                 d     = max(sqrt(dd+1.d-10),1.d-10)
                 CosTh = (dX*ux+dY*uy+dZ*uz)/d
                 d2 =dd*(1.-CosTh**2)
                 If(d2.lt. R3)Then  ! inside the cylinder
c                 If(CosTh.gt. CosCone)Then  ! inside the cone
                    Nc =Nc +1
                wx =Vxp(IN)-Vxc
                wy =Vyp(IN)-Vyc
                wz =Vzp(IN)-Vxc
                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
                d =d*ScaleC
                dX =dX*ScaleC
                 dY =dY*ScaleC
                 dZ =dZ*ScaleC
              write (18,100) d,vr,ww,dX,dY,dZ,wx,wy,wz,Wp(IN)
 100            format(f8.2,2f9.2,3f8.3,3f8.1,g11.3,4f8.3)
                iL =max(min(int(d/dRline),100),0)
                Rline(iL) = Rline(iL) +d*Wp(IN)
                Wline(iL) = Wline(iL) +Wp(IN)
                Mline(iL) = Mline(iL) +1
                Vline(iL) = Vline(iL) +vr*Wp(IN)
                V2line(iL) = V2line(iL) +vr**2*Wp(IN)
c                write (*,*) d,dRline,
             endif 
             endif 
      Enddo 
      write (*,*) ' N_inside_cone=',Nc
      open(30,file='Line.dat')
      do i=0,99
         If(Mline(i).gt.3)Then
            Rline(i) =Rline(i)/Wline(i)
            Vline(i) = Vline(i)/Wline(i)
            V2line(i) =sqrt(V2line(i)/Wline(i)-Vline(i)**2)
            write(30,80) i,(i+0.5)*dRline,Rline(i),Vline(i), 
     &                        V2line(i),Mline(i),Wline(i)
 80         format(i3,f8.3,f8.3,2f9.2,i6,g11.3)
         EndIf
      Enddo 
      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.d-10),1.d-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