PROGRAM XEPHEM
C
C     EXERCISE DRIVER FOR SUBROUTINE EPHEM
C
C     B. KNAPP, 1994-06-14, 2000-10-12
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
      REAL*8 DAY, SEC, JDE, JD, ECL_A(6), EQU_A(6), HC_A(6),
     &   ECL_T(6), EQU_T(6), HC_T(6)
      INTEGER*4 P, YR, MON, HOUR, MINUTE, STATUS1, STATUS2, STATUS3,
     &   STATUS4, STATUS5, STATUS6
      CHARACTER*8 NAME(0:10)
     &   / 'Sun'    , 'Mercury', 'Venus'  , 'Mars'   , 'Earth'  ,
     &     'Jupiter', 'Saturn' , 'Uranus' , 'Neptune', 'Pluto'  ,
     &     'Moon'   /
      CHARACTER*6 DIST_UNIT
      REAL*8 C, AU, JM, AU2JM
      PARAMETER (C=299792.458D0, AU=149597870.66D0, JM=365250.D0,
     &   AU2JM = AU/C/86400.D0/JM)
      INTEGER*4 TRUE, APPARENT, EQUATORIAL, ECLIPTIC, HELIOCENTRIC,
     &   SPHERICAL
      PARAMETER (TRUE=1, APPARENT=2, EQUATORIAL=1, ECLIPTIC=2,
     &   HELIOCENTRIC=3, SPHERICAL=1)
C
C     EXTERNAL FUNCTIONS, SUBROUTINES
      REAL*8 TD2UT
      EXTERNAL TD2UT, YMD2JD, EPHEM
C
C
    1 CONTINUE
      WRITE(*,'(/A$)') ' Year, month, day, hour, min, sec (UT)? '
      READ(*,*) YR, MON, DAY, HOUR, MINUTE, SEC
      DAY = DAY + HOUR/24.D0 + MINUTE/1440.D0 + SEC/86400.D0
      CALL YMD2JD( YR, MON, DAY, JDE )
      JD = TD2UT(JDE)
C
    2 CONTINUE
      WRITE(*,'( A$)') ' Body (0,...,9,10, for Sun,...,Pluto,Moon)? '
      READ(*,*) P
      IF ( P .LT. 0 .OR. 10 .LT. P ) GOTO 2
C
      CALL EPHEM( JD, P, APPARENT, EQUATORIAL, SPHERICAL, 
     &   EQU_A, STATUS1 )
      CALL EPHEM( JD, P, APPARENT, ECLIPTIC, SPHERICAL,
     &   ECL_A, STATUS2 )
      CALL EPHEM( JD, P, APPARENT, HELIOCENTRIC, SPHERICAL,
     &   HC_A,  STATUS3 )
      CALL EPHEM( JD, P, TRUE, EQUATORIAL, SPHERICAL,
     &   EQU_T, STATUS4 )
      CALL EPHEM( JD, P, TRUE, ECLIPTIC, SPHERICAL,
     &   ECL_T, STATUS5 )
      CALL EPHEM( JD, P, TRUE, HELIOCENTRIC, SPHERICAL,
     &   HC_T,  STATUS6 )
C
C     PRINT THE RESULTS
      IF (P .EQ. 10) THEN
         DIST_UNIT = '(km)  '
         ECL_A(3) = ECL_A(3)*AU
         ECL_A(6) = ECL_A(6)*AU/86400.D0
         EQU_A(3) = EQU_A(3)*AU
         EQU_A(6) = EQU_A(6)*AU/86400.D0
         ECL_T(3) = ECL_T(3)*AU
         ECL_T(6) = ECL_T(6)*AU/86400.D0
         EQU_T(3) = EQU_T(3)*AU
         EQU_T(6) = EQU_T(6)*AU/86400.D0
      ELSE
         DIST_UNIT = '(A.U.)'
      ENDIF
C
      WRITE(*,5) STATUS1, STATUS2, STATUS3, STATUS4, STATUS5, STATUS6
 5    FORMAT(' Return statuses (1-6): ',6i6)
C
      WRITE(*,10) NAME(P), JD, ECL_A(1), HC_A(1),
     &                         ECL_A(2), HC_A(2),
     &              DIST_UNIT, ECL_A(3), HC_A(3),
     &                         EQU_A(1),
     &                         EQU_A(2)
   10 FORMAT(
     &   /' Ephemeris for ',A8,', JD =',F16.6/
     &   /'    (Apparent geocentric position, corrected for light ',
     &    'travel time,',
     &   /'     annual aberration, nutation, and FK5 system)'/
     &   /'    Ecliptic longitude (deg) = ',2F16.6,
     &   /'    Ecliptic latitude (deg)  = ',2F16.6,
     &   /'    Apparent distance ',A6,' = ',2F16.6,
     &   /'    Right ascension (deg)    = ',1F16.6,
     &   /'    Declination (deg)        = ',1F16.6// )
C
      WRITE(*,20) NAME(P), JD, ECL_T(1), HC_T(1),
     &                         ECL_T(2), HC_T(2),
     &              DIST_UNIT, ECL_T(3), HC_T(3),
     &                         EQU_T(1),
     &                         EQU_T(2)
   20 FORMAT(
     &   /' Ephemeris for ',A8,', JD =',F16.6/
     &   /'    (True geocentric position, corrected for',
     &   /'     annual aberration, nutation, and FK5 system)'/
     &   /'    Ecliptic longitude (deg) = ',2F16.6,
     &   /'    Ecliptic latitude (deg)  = ',2F16.6,
     &   /'    True distance ',A6,'     = ',2F16.6,
     &   /'    Right ascension (deg)    = ',1F16.6,
     &   /'    Declination (deg)        = ',1F16.6// )
C
      GOTO 1
      END