!***********************************************************************************************************************************
!
!                                                          S L U N A R
!
!  Program:      SLUNAR
!
!  Programmer:   Dr. David G. Simpson
!                NASA Goddard Space Flight Center
!                Greenbelt, Maryland  20771
!
!  Date:         July 29, 2002
!
!  Language:     Fortran-90
!
!  Version:      1.00c
!
!  Description:  This program implements a simple lunar ephemeris model.  The resulting lunar coordinates are referred to the
!                mean equator and equinox of epoch J2000.  (See 'An Alternative Lunar Ephemeris Model for On-Board Flight Software
!                Use', D.G. Simpson, Proceedings of the 1999 NASA/GSFC Flight Mechanics Symposium, pp. 175-184)
!
!***********************************************************************************************************************************



      PROGRAM SLUNAR

      IMPLICIT NONE

!-----------------------------------------------------------------------------------------------------------------------------------
!  Variable declarations
!-----------------------------------------------------------------------------------------------------------------------------------

      DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0

      DOUBLE PRECISION :: JD                                                        ! input time (ephemeris Julian day)
      DOUBLE PRECISION :: T                                                         ! time in Julian centuries from J2000
      DOUBLE PRECISION :: X, Y, Z                                                   ! lunar J2000 cartesian coordinates (km)
      DOUBLE PRECISION :: RA                                                        ! lunar J2000 right ascension (deg)
      DOUBLE PRECISION :: DEC                                                       ! lunar J2000 declination (deg)



!-----------------------------------------------------------------------------------------------------------------------------------
!  Start of code
!-----------------------------------------------------------------------------------------------------------------------------------

!
!     Read input time.
!

      WRITE (UNIT=*, FMT='(A)') ' Enter time (JD): '                                ! read desired time (JD)
      READ (UNIT=*, FMT=*) JD

      T = (JD - 2451545.0D0)/36525.0D0                                              ! convert time to Julian centuries from J2000

!
!     Compute lunar J2000 cartesian coordinates.
!

      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)

!
!     Print results.
!

      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'

!
!     Compute right ascension and declination.
!

      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

!
!     Print results.
!

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

      STOP

      END PROGRAM SLUNAR