      PROGRAM VENUS_EARTH
C
C     GET DAILY HELIOCENTRIC COORDINATES OF VENUS AND EARTH OVER A
C     GIVEN RANGE OF TIME, USING VSOP87 THEORY.
C
C     B. KNAPP, 1993-05-28, 1999-06-25, 2001-01-03, 2001-04-25
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
      INTEGER*4 YD, J, IERRE, IERRV
      REAL*8 YR, DAY, YD_BEG, YD_END, JD, JDE, T, UT2TD,
     &   JD_BEG, PV_V(6), PV_E(6)
      REAL*8 D2R
      PARAMETER (D2R = 3.14159265358979324D0/180.D0)
C
C     CONSTANTS
      INTEGER*4 VENUS, EARTH
      PARAMETER (VENUS=2, EARTH=3)
C
C     EXTERNAL FUNCTIONS, SUBROUTINES
      REAL*8 DYD, YD2
      EXTERNAL UT2TD, YD2JD, DYD, YD2, HCPVD
C
C
      WRITE(*,'(/A$)') ' Beginning Year, day? '
      READ(*,*) YR, DAY
      YD_BEG = 1000.D0*YR+DAY
      CALL YD2JD( YD_BEG, JD_BEG )
C
      WRITE(*,'(A$)')  '    Ending Year, day? '
      READ(*,*) YR, DAY
      YD_END = 1000.D0*YR+DAY
C
C     LOOP OVER GIVEN TIME RANGE
      DO J = 0,NINT( DYD( YD_BEG, YD_END ) )

         YD = NINT( YD2( YD_BEG, DBLE( J ) ) )
         JD = JD_BEG + DBLE( J )
         JDE = UT2TD( JD )
C
C        GET DYNAMIC TIME (IN JULIAN MILLENIA SINCE J2000)
         T = (JDE-2451545.D0)/365250.D0
C
C        GET APPARENT HELIOCENTRIC ECLIPTIC LONGITUDE AND LATITUDE
         CALL HCPVD( VENUS, T, PV_V, IERRV )
         CALL HCPVD( EARTH, T, PV_E, IERRE )
C
C        PRINT THE RESULTS
         PV_V(1) = MOD( PV_V(1)/D2R, 360D0 )
         IF ( PV_V(1) .LT. 0. ) PV_V(1) = PV_V(1) + 360.
         PV_E(1) = MOD( PV_E(1)/D2R, 360D0 )
         IF ( PV_E(1) .LT. 0. ) PV_E(1) = PV_E(1) + 360.
         PV_V(2) = PV_V(2)/D2R
         PV_E(2) = PV_E(2)/D2R
C
         WRITE(*,10)
     &      YD, PV_V(1), PV_V(2), PV_V(3), PV_E(1), PV_E(2), PV_E(3)
   10    FORMAT( I10, F9.2, F8.2, F9.5, 2X, F9.2, F8.2, F9.5 )
      ENDDO
C
      END
