subroutine ExtractPositionn8(ITP,iYr,Doy8,RVect)

            parameter               (NMAX = 7000)
            integer                 ITP
            integer                 iYr
            double precision        Doy8
            real                    RVect(3)

       	real                    VVect(5)

            real                    error(3) ! error estimate for interpolation
            real                    dates(2,NMAX)
            real                    coords(3,NMAX)

            double precision        dLngSun
            double precision        dLatSun
            double precision        dDisSun
            
            character(len=100)      :: parker_file,bepi_file,solar_file,stereoa_file
            integer                 :: NTPS, io

        if (ITP .eq. 3) then            ! Earth
            call SunNewcomb8(0,iYr,Doy8,dLngSun,dLatSun,dDisSun)
            RVect(1) =  dLngSun+180
            RVect(2) = -dLatSun/3600    ! arcsec -> degrees
            RVect(3) =  dDisSun

        else if (ITP .eq. 10) then      ! Ulysses
            Doy = sngl(Doy8)
            call UlyssesOrbit(iYr,Doy,RVect,VVect)

        else if (ITP .eq. 11 .or.       ! Helios 1
     &           ITP .eq. 12) then      ! Helios 2
            Doy = sngl(Doy8)
            call HOSOrbit(ITP-10,iYr,Doy,RVect(3),RVect(1),VR,VT)
            RVect(2) = 0

        else if (ITP .eq. 13) then      ! STEREO A
            Doy = sngl(Doy8)
            if (iYr.lt.2021) then
              call StereoAOrbit(iYr,Doy,RVect,VVect)
            else            
              stereoa_file = '/home/soft/dat/scposition/stereoa_trajectory.txt'
              stereoa_file = trim(stereoa_file)            
              call read_locations(NMAX,stereoa_file,NPTS,dates,coords)
              call get_orbit(ITP, iYr, Doy, NPTS, NMAX, dates, coords, RVect, error)
            end if
        else if (ITP .eq. 14) then      ! STEREO B
            Doy = sngl(Doy8)
            call StereoBOrbit(iYr,Doy,RVect,VVect)

        else if (ITP .eq. 15) then      ! MESSENGER
            Doy = sngl(Doy8)
            call Messengerorbit(iYr,Doy,RVect,VVect)

        else if (ITP .eq. 16) then      ! Comet 67P/GC
            Doy = sngl(Doy8)
            call comet67PCGorbit(iYr,Doy,RVect,VVect)
            
        else if (ITP .eq. 17) then
            Doy = sngl(Doy8)
            parker_file = '/home/soft/dat/scposition/psp_trajectory.txt'
            parker_file = trim(parker_file)            
            call read_locations(NMAX,parker_file,NPTS,dates,coords)
            call get_orbit(ITP, iYr, Doy, NPTS, NMAX, dates, coords, RVect, error)

        else if (ITP .eq. 18) then
            Doy = sngl(Doy8)
            bepi_file = '/home/soft/dat/scposition/bc_trajectory.txt'
            bepi_file = trim(bepi_file)
            call read_locations(NMAX,bepi_file,NPTS,dates,coords)
            call get_orbit(ITP, iYr, Doy, NPTS, NMAX, dates, coords, RVect, error)

        else if (ITP .eq. 19) then
            Doy = sngl(Doy8)
            solar_file = '/home/soft/dat/scposition/so_trajectory.txt'
            solar_file = trim(solar_file)
            call read_locations(NMAX,solar_file,NPTS,dates,coords)
            call get_orbit(ITP, iYr, Doy, NPTS, NMAX, dates, coords, RVect, error)
            
        else
            Doy = sngl(Doy8)
            call PlanetOrbit(ITP,iYr,Doy,RVect,VVect)

        end if
C        write(*,'(A,3F9.3)'), 'In extractpositionn8.f Rvect(1), Rvect(2), Rvect(3)',RVect(1), RVect(2), RVect(3)
                                        ! Convert to heliographic coordinates
        call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,RVect(1),RVect(2))

        return
        end