real x,y,z,r(15000),vx,vy,vz,v2,vt2,vr,vt(15000)
      integer rindex(15000)

      real rin,rout,sumvt,sumvr,meanvt,meanr,numdens,pi
      integer i,ntot,nbin,ninbin,ibin,ipart,inum,partid

      ninbin=100     ! number of particles per radial bin
      pi=3.1415927

C open data file; read coords, calculate r, v_tangential

      open (unit=10,file='data.dat')
      i=0
 1    read (10,*,end=2) x,y,z,vx,vy,vz
      i=i+1
      r(i)=sqrt(x*x+y*y)   ! cylindrical R
      vr=(vx*x+vy*y)/r(i)  ! radial velocity
      v2=vx*vx+vy*vy       ! velocity squared
      vt2=v2-vr*vr         ! v_tangential squared
      if (vt2.gt.0) then   ! if we can, calculate vt
         vt(i)=sqrt(vt2)
      else                 ! if not, set to zero
         vt(i)=0.0
      endif
      go to 1

 2    close (10)
      ntot=i                    ! total particle number


C sort particles in R using an index variable rindex (see num rec)

      call indexx(ntot,r,rindex)

C now bin particles and calculate meanr, meanvt, and numdens

      nbin=ntot/ninbin       ! calculate radial bins of ninbin particles each
      rin=0                  ! inner radius of 1st bin is 0
      open (unit=10,file='diskprop.dat')
      do ibin=1,nbin         ! for all bins...
         sumvt=0             ! initialize the sum of vt's
         sumr=0              ! initialize the sum of r's
         do ipart=1,ninbin   ! for each of the ninbin particles in the bin...
            inum=(ibin-1)*ninbin+ipart  ! calculate where it ranks in r
            partid=rindex(inum)         ! find the particle id 
            sumvt=sumvt+vt(partid)      ! add its vt to the bin sum
            sumr=sumr+r(partid)         ! add its r to the bin sum
         enddo        
         meanvt=sumvt/ninbin  ! calculate average v_tangential for bin
         meanr=sumr/ninbin    ! calculate mean r for bin
         rout=r(partid)       ! find outermost particle in bin
         numdens=ninbin/(pi*(rout*rout-rin*rin))   ! density in the bin
         write (10,*) meanr, numdens, meanvt       ! write out the data
         rin=rout             ! for next bin set inner radius
      enddo
      close (10)

      end

C==========================================================================

      SUBROUTINE INDEXX(N,ARRIN,INDX)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END