subroutine Spacecraft_PV( SpacecraftID, pathlen, path, JD,
     &   PV, dPV, f, Status )
!
!     Given NORAD (US SpaceCom) spacecraft ID and Julian Date (and
!     fraction), return the position and velocity (x, y, z, x', y',
!     z' and 1 s.d. uncertainties) of the spacecraft in the Earth-
!     centered inertial frame, using SGP4 and the NORAD two-line
!     element sets.  Also, return f, the interpolation fraction,
!     which defines the weighting used to construct the result from
!     the two nearest NORAD TLE sets: result = f*PV(1) + (1-f)*PV(2).
!     (F will be in the range 0 to 1 if the requested JD is bracketed
!     by the epochs of the two TLE sets used.)
!
!     B. Knapp, 1999-04-20, 2000-05-02, 2000-05-26, 2000-06-02
!               2003-11-20, Make path an input variable
!
      implicit none

!     Input
      integer*4 SpacecraftID, pathlen
      character*(1024) path
      real*8 JD
!
!     Output
      real*8 PV(6), dPV(6), f
      integer Status
!
!     Static
      integer*4 MAX_SETS, NSets, NGood, prevID/-1/, jl, ju
      parameter (MAX_SETS=8192)
      integer*4 EltSet(MAX_SETS), Orbit(MAX_SETS)
      real*8 EpochJD(MAX_SETS), dNdT(MAX_SETS), d2NdT2(MAX_SETS),
     &   bStar(MAX_SETS), Incl(MAX_SETS), Node(MAX_SETS), E(MAX_SETS),
     &   Omega(MAX_SETS), M(MAX_SETS), N(MAX_SETS)
      character*8 IntlNo(MAX_SETS)
      character*1 ClsLvl(MAX_SETS), EphTyp(MAX_SETS)
      save EltSet, Orbit, EpochJD, dNdT, d2NdT2, bStar, Incl, Node, E,
     &   Omega, M, N, IntlNo, ClsLvl, EphTyp, NSets, NGood, prevID,
     &   jl, ju
!
!     Local
      integer LUN, fStatus, j, jm, j0, step, SatID
      real*8 EpochYD
      character*12 fname
!
!     Functions, subroutines
      external Get_LUN, Free_LUN, NextTle, yd2jd, TLE_Subset,
     &   PV_from_Elts
!
!     Initialize TLE sets?
      if ( SpacecraftID .ne. prevID ) then
!
         write(fname,10) SpacecraftID
 10      format(i8.8,'.tle')
         call Get_LUN( LUN, fStatus )
         if ( fStatus .ne. 0 ) then
            Status=fStatus+8000
            return
         endif
!         print *, path(1:pathlen)//fname
         open( LUN, file=path(1:pathlen)//fname, status='old', 
     &      form='formatted', iostat=fStatus )
         if ( fStatus .ne. 0 ) then
            Status=fStatus+9000
            return
         endif
!
         j = 1
         Status = 0
         do while( j .lt. MAX_SETS .and. Status .ge. 0 )
            call NextTle( lun, SatID, ClsLvl(j), IntlNo(j), EltSet(j),
     &         EphTyp(j), Orbit(j), EpochYD, dNdT(j), d2NdT2(j),
     &         bStar(j), Incl(j), Node(j), E(j), Omega(j), M(j), N(j),
     &         Status )
            if ( Status .eq. 0 ) then
               N(j) = N(j)/1440.d0
               if ( EpochYD/1000.d0 .ge. 57.d0 ) then
                  EpochYD = 1900000.d0+EpochYD
               else
                  EpochYD = 2000000.d0+EpochYD
               endif
               call yd2jd( EpochYD, EpochJD(j) )
               NSets = j
               j = j+1
            endif
         enddo
         close( LUN )
         call Free_LUN( LUN, fStatus )
!
!        How did we end?
         if ( Status .eq. 0 ) then
            Status = 99  !Too many element sets
            return
         else if ( Status .gt. 0 ) then
            return       !Problem reading elements
         endif
!
!        Select a consistent subset
         call TLE_Subset( NSets, NGood, ClsLvl, IntlNo, EltSet, EphTyp,
     &      Orbit, EpochJD, dNdT, d2NdT2, bStar, Incl, Node, E,
     &      Omega, M, N, Status )
         if ( Status .ne. 0 ) return

!        Initialization successful!
         jl = 1
         ju = NGood
         prevID = SpacecraftID
      endif
!
!     Binary search for j s.t. EpochJD(j) <= JD < EpochJD(j+1)
!
!     Expand search interval forward?
      step = 1
      do while ( ju .lt. NGood .and. JD .ge. EpochJD(ju) )
         ju = min( NGood, ju+step )
         step = 2*step
      enddo
!
!     Expand backward?
      step = 1
      do while ( jl .gt. 1 .and. JD .lt. EpochJd(jl) )
         jl = max( jl-step, 1 )
         step = 2*step
      enddo
!
!     Shrink interval?
      do while ( ju-jl .gt. 1 )
         jm = (ju+jl)/2
         if ( JD .ge. EpochJD(jm) ) then
            jl = jm
         else
            ju = jm
         endif
      enddo
!
!     Pass four sets (jl-1, jl, jl+1, jl+2) to the orbit propagator
!     (but take care if we are near one end of the element set arrays)
      j0 = max( 1, min( max( 1, jl-1 ), NGood-3 ) )
      call PV_from_Elts( JD, EpochJD(j0), bStar(j0), Incl(j0), Node(j0),
     &   E(j0), Omega(j0), M(j0), N(j0), PV, dPV, f, Status )
!
      return
      end