subroutine neighbor_consistency( NSets, Epoch, Incl, Node, E,
     &   Omega, M, N, bStar, j, cR, cV )
!
!     Compute the consistency in position and velocity of TLE set j
!     with its neighbors.
!
!     B. Knapp, 2000-05-04
!
      implicit none
!
!     Input
      integer*4 NSets, j
      real*8 Epoch(NSets), Incl(NSets), Node(NSets), E(NSets),
     &   Omega(NSets), M(NSets), N(NSets), bStar(NSets)
!
!     Output
      real*8 cR, cV
!
!     Local
      real*8 dT12, dR12, dV12, dRdT12, dVdT12,
     &       dT23, dR23, dV23, dRdT23, dVdT23,
     &       dT13, dR13, dV13, dRdT13, dVdT13
!
!     Externals
      external Pair_Difference
!
!     Can't do 1st or last
      if ( j .eq. 1 .or. j .eq. NSets ) then
         cR = 0.
         cV = 0.
      else
!
!        Get 3 distances (in both dRdT and dVdt)
         call Pair_Difference( NSets, Epoch, Incl, Node, E, Omega,
     &      M, N, bStar, j-1, j,   dT12, dR12, dV12, dRdT12, dVdT12 )
         call Pair_Difference( NSets, Epoch, Incl, Node, E, Omega,
     &      M, N, bStar, j  , j+1, dT23, dR23, dV23, dRdT23, dVdT23 )
         call Pair_Difference( NSets, Epoch, Incl, Node, E, Omega,
     &      M, N, bStar, j-1, j+1, dT13, dR13, dV13, dRdT13, dVdT13 )
!
!        Set j is inconsistent if sets j-1 and j+1 agree with each
!        other, but not with j:
         cR = (dRdT12 + dRdT23 + dRdT13)/(3.*dRdT13)
         cV = (dVdT12 + dVdT23 + dVdT13)/(3.*dVdT13)
      endif
      return
      end