C$Procedure                         LUNORB
C
      SUBROUTINE LUNORB ( DATE , MOON )
C
C
C*******************************************************************************
C
C       Copyright (C) 1993, California Institute of Technology.  U.S.
C       Government Sponsorhip under NASA Contract NAS7-918 is
C       acknowledged.
C
C*******************************************************************************
C
C$ Purpose
C
C  LUNORB computes the mean orbital elements of the Moon in EME and equinox
C  of date coordinate system.
C
C$ Input_Arguments
C
C  DATE   - epoch at which to determine orbital elements of Sun
C
C$ Output_Arguments
C
C  MOON   - mean orbital elements of the Moon (a,e,i,LAN,w,M), km,deg
C
C$ Log
C
C  21-Mar-1989 - Eric Cannell     - creation
C
C$ Library_Links
C
C  TIMETRANS
C
C$ Parameters
C
      CHARACTER*25         J1899
      PARAMETER          ( J1899  = '31-DEC-1899 12:00:00.0000' )

      DOUBLE PRECISION     PI
      PARAMETER          ( PI     = 3. 14159 26535 89793 23846 D0 )

      DOUBLE PRECISION     D2R
      PARAMETER          ( D2R    = PI / 180.0D0 )
C
C$ Declarations_of_Input_and_Output_Arguments
C
      CHARACTER*(*)        DATE
      DOUBLE PRECISION     MOON   ( 6 )
C
C$ Declarations_of_Local_Variables
C
      DOUBLE PRECISION     ALPHA
      DOUBLE PRECISION     APM
      DOUBLE PRECISION     COSALF
      DOUBLE PRECISION     COSB
      DOUBLE PRECISION     COSOM
      DOUBLE PRECISION     D
      DOUBLE PRECISION     DF                            
      DOUBLE PRECISION     EPM
      LOGICAL              ERROR
      DOUBLE PRECISION     IPM
      DOUBLE PRECISION     JD1899
      DOUBLE PRECISION     JDEPOC
      DOUBLE PRECISION     MPM
      DOUBLE PRECISION     OOE
      DOUBLE PRECISION     OPM
      DOUBLE PRECISION     SINB
      DOUBLE PRECISION     SINOM
      DOUBLE PRECISION     TE
      DOUBLE PRECISION     WPM
C           
C$ External_Statements
C
      DOUBLE PRECISION     CH2JD
      EXTERNAL             CH2JD
C
C$ Method
C-&

C1    Convert reference epoch and orbit epoch to Julian date.

      JD1899 = CH2JD ( J1899 , .TRUE. , ERROR )
      JDEPOC = CH2JD ( DATE  , .TRUE. , ERROR )

C1    Compute number of ephemeris days from reference epoch and fraction.

      D  = JDEPOC - JD1899
      DF = D / 10000.0D0

C1    Compute orbital elements of Moon (a, e, i, LAN , w, M) in Earth
C1    mean ecliptic and equinox of date.

      APM = 384 400.0D0
      EPM = 0.054900489D0
      IPM = 5.1453964D0 
      OPM = 259.183275D0 - 0.0529539222*D + .0001557*DF**2 +
     &      5.D-8*DF**3
      WPM = 334.329556D0 + .1114040803D0*D - .0007739D0*DF**2 - 
     &      2.6D-7*DF**3 - OPM
      MPM = 270.434164D0 + 13.1763965268*D - 8.5D-5*DF**2 +
     &      3.9D-8*DF**3 - OPM - WPM

      OPM = DMOD( OPM , 360.0D0 )
      WPM = DMOD( WPM , 360.0D0 )
      MPM = DMOD( MPM , 360.0D0 )

      IF ( OPM .LT. 0.0D0 ) OPM = OPM + 360.0D0
      IF ( WPM .LT. 0.0D0 ) WPM = WPM + 360.0D0
      IF ( MPM .LT. 0.0D0 ) MPM = MPM + 360.0D0

C1    Compute the mean obliquity of the ecliptic.

      TE  = ( D + 1.0D0 ) / 36525.0D0

      OOE = 23.45229444D0 - .0130125D0 * TE - .16388889D-5*TE**2 +
     &      .50277778D-6*TE**3

      OOE = OOE * D2R

C1    Copy a, e, and M into return vector.

      MOON( 1 ) = APM
      MOON( 2 ) = EPM
      MOON( 6 ) = MPM

C1    Convert inclination from ecliptic to equator.

      COSALF =   DSIN(IPM*D2R)*DSIN(OOE)*DCOS(OPM*D2R)
     &         - DCOS(OOE)*DCOS(IPM*D2R)
      ALPHA  = DACOS( COSALF )

      MOON( 3 ) = ( PI - ALPHA ) / D2R

C1    Convert LAN from ecliptic to equator.

      SINOM = DSIN(OPM*D2R) * DSIN(IPM*D2R) / DSIN(ALPHA)
      COSOM =   ( DCOS(IPM*D2R) + DCOS(ALPHA) * DCOS(OOE) ) 
     &        / DSIN(ALPHA) / DSIN(OOE)

      IF ( SINOM .EQ. 0.0D0 .AND. COSOM .EQ. 0.0D0 ) THEN
         MOON( 4 ) = 0.0
      ELSE
         MOON( 4 ) = DATAN2( SINOM , COSOM ) / D2R
         IF ( MOON( 4 ) .LT. 0.0D0 ) MOON( 4 ) = MOON( 4 ) + 360.0D0
      END IF

C1    Convert w from ecliptic to equator.

      SINB = DSIN(OPM*D2R) * DSIN(OOE) / DSIN(ALPHA)
      COSB =   ( DCOS(OOE) + DCOS(ALPHA) * DCOS(IPM*D2R) ) 
     &        / DSIN(ALPHA) / DSIN(IPM*D2R)

      IF ( SINB .EQ. 0.0D0 .AND. COSB .EQ. 0.0D0 ) THEN
         MOON( 5 ) = WPM
      ELSE
         MOON( 5 ) = WPM + DATAN2( SINB , COSB ) / D2R
         IF ( MOON( 5 ) .LT. 0.0D0 ) MOON( 5 ) = MOON( 5 ) + 360.0D0
      END IF

      RETURN
      END