subroutine PV_from_Elts( JD, EpochJD,
     &   bStar, Incl, Node, E, Omega, M, N, PV, dPV, f, Status )
!
!     Given four (consecutive) sets of NORAD (USSPACECOM) two-line
!     orbit elements and a Julian Day Number, select the pair which
!     yields the mimimum error in propagating with SGP4, and return
!     the corrected position and velocity, and 1-s.d. uncertainties.
!
!     B. Knapp, 2000-05-26
!
      implicit none
!
!     Input
      real*8 JD, EpochJD(4), bStar(4), Incl(4), Node(4), E(4), Omega(4),
     &   M(4), N(4)
!
!     Output
      real*8 PV(6), dPV(6), f
!
!     Externals
      real*8 ut2td
      external SGP4, ut2td, Nutate_Matrix
!
!     Static: none
!
!     Local
      integer*4 i, j, IFlag, Status
      real*8 fj, x, w, Tl, Tu, PVl(6), PVu(6), dPmag, dPmin, td,
     &   p(3), v(3), nu(3,3)
!
!     To estimate the uncertainty in the result, we will propagate to
!     the given JD from the two "best" TLE sets, defined as the two
!     with adjacent epochs that minimize the error, the rationale 
!     being that these are two most "consistent" TLE sets.
!     
!     The error is a function of the difference Delta between these two
!     propagations, and the distance in time from the nearest TLE set.
!     The best case is that JD coincides exactly with one of the two TLE
!     epochs, in which case we estimate the uncertainty to be  Delta/2
!     (the minimum).  If JD is equidistant from two bracketing TLE sets,
!     we estimate the uncertainty to be Delta (a local maximum). These 
!     conditions define a quartic: dy = (1-4x^2(1-2x^2))*Delta, where x
!     is the ratio of the distance in time from the midpoint of the two
!     TLE epoch times to the difference in those epoch times, and y is
!     one of the 6 position, velocity magnitudes.
!
!     Initialize output
      do i=1,6
         PV(i) = 0.d0
      enddo
!
!     Search for best pair of TLE sets
      dPmin = 1.d38
      do j=1,3
         fj = (JD-EpochJD(j))/(EpochJD(j+1)-EpochJD(j))
         x = ((fj-(1.d0-fj))**2)/2.d0
         w = 1.-2.d0*x*(1.d0-x)
!
         Tl = (JD-EpochJD(j))*1440.d0
         Tu = (JD-EpochJD(j+1))*1440.d0
         IFlag = 1
         call SGP4( IFlag, Tl, M(j), Node(j), Omega(j), E(j),
     &      Incl(j), N(j), bStar(j), PVl, Status )
         IFlag = 1
         call SGP4( IFlag, Tu, M(j+1), Node(j+1), Omega(j+1), E(j+1),
     &      Incl(j+1), N(j+1), bStar(j+1), PVu, Status )
         dPmag = 0.d0
         do i=1,3
            dPmag = dPmag + ((PVu(i)-PVl(i))*w)**2
         enddo
!
         if ( dPmag .lt. dPmin ) then
            do i=1,6
               PV(i) = PVl(i)*(1.d0-fj)+PVu(i)*fj
               dPV(i) = abs( PVu(i)-PVl(i) )*w
            enddo
            f = fj
            dPmin = dPmag
         endif
      enddo
!
!     Make nutation correction
      td = (ut2td( JD )-2451545.d0)/365250.d0
      call Nutate_Matrix( td, nu )
      do j=1,3
         p(j) = PV(j)
         v(j) = PV(j+3)
      enddo
      do i=1,3
         PV(i) = 0.d0
         PV(i+3) = 0.d0
         do j=1,3
            PV(j) = PV(j)+nu(i,j)*p(j)
            PV(j+3) = PV(j+3)+nu(i,j)*v(j)
         enddo
      enddo
!
      return
      end