C***********************************************************************
C
C                             S L U N A R
C
C  PROGRAM:      SLUNAR
C
C  PROGRAMMER:   DR. DAVID G. SIMPSON
C                NASA GODDARD SPACE FLIGHT CENTER
C                GREENBELT, MARYLAND  20771
C
C  DATE:         JULY 29, 2002
C
C  LANGUAGE:     FORTRAN-77
C
C  VERSION:      1.00C
C
C  DESCRIPTION:  THIS PROGRAM IMPLEMENTS A SIMPLE LUNAR EPHEMERIS MODEL.
C                THE RESULTING LUNAR COORDINATES ARE REFERRED TO THE
C                MEAN EQUATOR AND EQUINOX OF EPOCH J2000.  (SEE
C                'AN ALTERNATIVE LUNAR EPHEMERIS MODEL FOR ON-BOARD
C                FLIGHT SOFTWARE USE', D.G. SIMPSON, PROCEEDINGS OF
C                THE 1999 NASA/GSFC FLIGHT MECHANICS SYMPOSIUM,
C                PP. 175-184)
C
C***********************************************************************



      PROGRAM SLUNAR

C-----------------------------------------------------------------------
C
C  VARIABLE DECLARATIONS
C
C PI       PI = 3.14159...
C JDE      INPUT TIME (EPHEMERIS JULIAN DAY)
C T        TIME IN JULIAN CENTURIES FROM J2000
C X,Y,Z    LUNAR J2000 CARTESIAN COORDINATES (KM)
C RA       LUNAR J2000 RIGHT ASCENSION (DEG)
C DEC      LUNAR J2000 DECLINATION (DEG)
C
C-----------------------------------------------------------------------

      DOUBLE PRECISION  PI
      PARAMETER (PI=3.14159265358979323846264338327950288419716939938D0)

      DOUBLE PRECISION  JDE
      DOUBLE PRECISION  T
      DOUBLE PRECISION  X, Y, Z
      DOUBLE PRECISION  RA
      DOUBLE PRECISION  DEC



C-----------------------------------------------------------------------
C  START OF CODE
C-----------------------------------------------------------------------

C
C     READ INPUT TIME (JDE)
C

      WRITE (UNIT=*, FMT='(A)') ' ENTER TIME (JDE): '
      READ (UNIT=*, FMT=*) JDE

C
C     CONVERT TIME TO JULIAN CENTURIES FROM J2000
C

      T = (JDE - 2451545.0D0)/36525.0D0

C
C     COMPUTE LUNAR J2000 CARTESIAN COORDINATES
C

      X =  383.0D3 * SIN( 8399.685D0 * T + 5.381D0)
     $    + 31.5D3 * SIN(   70.990D0 * T + 6.169D0)
     $    + 10.6D3 * SIN(16728.377D0 * T + 1.453D0)
     $    +  6.2D3 * SIN( 1185.622D0 * T + 0.481D0)
     $    +  3.2D3 * SIN( 7143.070D0 * T + 5.017D0)
     $    +  2.3D3 * SIN(15613.745D0 * T + 0.857D0)
     $    +  0.8D3 * SIN( 8467.263D0 * T + 1.010D0)

      Y =  351.0D3 * SIN( 8399.687D0 * T + 3.811D0)
     $    + 28.9D3 * SIN(   70.997D0 * T + 4.596D0)
     $    + 13.7D3 * SIN( 8433.466D0 * T + 4.766D0)
     $    +  9.7D3 * SIN(16728.380D0 * T + 6.165D0)
     $    +  5.7D3 * SIN( 1185.667D0 * T + 5.164D0)
     $    +  2.9D3 * SIN( 7143.058D0 * T + 0.300D0)
     $    +  2.1D3 * SIN(15613.755D0 * T + 5.565D0)

      Z =  153.2D3 * SIN( 8399.672D0 * T + 3.807D0)
     $    + 31.5D3 * SIN( 8433.464D0 * T + 1.629D0)
     $    + 12.5D3 * SIN(   70.996D0 * T + 4.595D0)
     $    +  4.2D3 * SIN(16728.364D0 * T + 6.162D0)
     $    +  2.5D3 * SIN( 1185.645D0 * T + 5.167D0)
     $    +  3.0D3 * SIN(  104.881D0 * T + 2.555D0)
     $    +  1.8D3 * SIN( 8399.116D0 * T + 6.248D0)

C
C     PRINT RESULTS
C

      WRITE (UNIT=*, FMT='(/A/)') ' LUNAR J2000 EQUATORIAL COORDINATES:'
      WRITE (UNIT=*, FMT='(A,F10.1,A)') ' X = ', X, ' KM'
      WRITE (UNIT=*, FMT='(A,F10.1,A)') ' Y = ', Y, ' KM'
      WRITE (UNIT=*, FMT='(A,F10.1,A)') ' Z = ', Z, ' KM'

C
C     COMPUTE RIGHT ASCENSION AND DECLINATION
C

      RA = ATAN2(Y,X) * 180.0D0/PI
      IF (RA .LT. 0.0D0) RA = RA + 360.0D0

      DEC = ASIN(Z/SQRT(X**2+Y**2+Z**2)) * 180.0D0/PI

C
C     PRINT RESULTS
C

      WRITE (UNIT=*, FMT='(A,F7.1,A)') ' RA  = ', RA, ' DEG'
      WRITE (UNIT=*, FMT='(A,F7.1,A)') ' DEC = ', DEC, ' DEG'

      STOP

      END