      SUBROUTINE MOON (PHIP, LAMBDAP, YR, YDAY, UT, ALPHAP, DELTAP,
     &                 ALTITUDE, AZIMUTH)
*
*	Subroutine MOON accepts observer's coordinates and a time argument, 
*	and uses low-precision formulae (error < .3 deg) to calculate 
*	the apparent lunar position, including right ascension, declination,
*	altitude and azimuth, for those coordinates at that time. The 
*     	formulae come from The Astronomical Almanac, 1984, p. D46, with the 
*	exception of the formula for the time argument, T, which is borrowed 
*	from the reference cited in SUN.FOR, and the formulae for elevation
*	and azimuth, which come from The Observer's Handbook, 1986.
*
      REAL     YR, YDAY, UT
      REAL * 8 PHIP 	! observer's geographic N latitude
      REAL * 8 LAMBDAP  ! observer's geographic E longitude
      REAL * 8 FDAY 	! fraction of elapsed day
      REAL * 8 JD 	! Julian day - 2451545.0
      REAL * 8 T 	! number of Julian centuries from J2000.0
      REAL * 8 LAMBDA 	! ecliptic longitude (degrees)
      REAL * 8 BETA 	! ecliptic latitude (degrees)
      REAL * 8 PI 	! horizontal parallax (degrees)
      REAL * 8 R 	! distance (Earth radii)
      REAL * 8 L, M, N 	! geocentric direction cosines
      REAL * 8 X, Y, Z 	! geocentric rectangular coordinates
      REAL * 8 XP, YP, ZP 	! topocentric rectangular coordinates 
      REAL * 8 THETA0 	! local sidereal time
      REAL * 8 RP 	! topocentric distance
      REAL * 8 ALPHAP 	! topocentric right ascension
      REAL * 8 DELTAP 	! topocentric declination
      REAL * 8 HANGLE 	! local lunar hour angle
      REAL * 8 ALTITUDE ! local lunar elevation angle
      REAL * 8 AZIMUTH	! local lunar azimuth
*                 
      FDAY = UT / 24.0
      IYR = INT (YR)
      JD = DBLE (REAL (365 * (IYR - 1900) + (IYR - 1901) / 4)) + 
     &                                         DBLE (YDAY) + FDAY - 0.5
      T = (JD / 36525.0) - 1.0
*
      LAMBDA = + 218.32 + 481267.883 * T 
     &         + 6.29 * DSIND (DBLE (134.9 + 477198.85 * T))
     &         - 1.27 * DSIND (DBLE (259.2 - 413335.38 * T))
     &         + 0.66 * DSIND (DBLE (235.7 + 890534.23 * T))
     &         + 0.21 * DSIND (DBLE (269.9 + 954397.70 * T))
     &         - 0.19 * DSIND (DBLE (357.5 +  35999.05 * T))
     &         - 0.11 * DSIND (DBLE (186.6 + 966404.05 * T))
*
      BETA =   + 5.13 * DSIND (DBLE ( 93.3 + 483202.03 * T))
     &         + 0.28 * DSIND (DBLE (228.2 + 960400.87 * T))
     &         - 0.28 * DSIND (DBLE (318.3 +   6003.18 * T))
     &         - 0.17 * DSIND (DBLE (217.6 - 407332.20 * T))
*
      PI =   + 0.9508
     &       + 0.0518 * DCOSD (DBLE (134.9 + 477198.85 * T))
     &       + 0.0095 * DCOSD (DBLE (259.2 - 413335.38 * T))
     &       + 0.0078 * DCOSD (DBLE (235.7 + 890534.23 * T))
     &       + 0.0028 * DCOSD (DBLE (269.9 + 954397.70 * T))
*
      R = 1.0 / DSIND (PI)
*
      L = DCOSD (BETA) * DCOSD (LAMBDA)
      M = 0.9175 * DCOSD (BETA) * DSIND (LAMBDA) - 0.3978 * DSIND (BETA)
      N = 0.3978 * DCOSD (BETA) * DSIND (LAMBDA) + 0.9175 * DSIND (BETA)
*
      X = R * L
      Y = R * M
      Z = R * N
*          
      THETA0 = 100.46 + 36000.77 * T + LAMBDAP + 15.0 * UT
*
      XP = X - DCOSD (PHIP) * DCOSD (THETA0)
      YP = Y - DCOSD (PHIP) * DSIND (THETA0)
      ZP = Z - DSIND (PHIP)
*
      RP = DSQRT (XP * XP + YP * YP + ZP * ZP)
      ALPHAP = DATAN2D (YP, XP)
      IF (ALPHAP .LT. 0.0) ALPHAP = ALPHAP + 360.0
      DELTAP = DASIND (ZP / RP)
*
      HANGLE = DMOD (THETA0 - ALPHAP, 360.0D0)
      IF (HANGLE .LT. 0.0) HANGLE = HANGLE + 360.0
      ALTITUDE = DASIND (DSIND (DELTAP) * DSIND (PHIP) 
     &          + DCOSD (HANGLE) * DCOSD (DELTAP) * DCOSD (PHIP))
      AZIMUTH = DASIND (- DCOSD (DELTAP) * DSIND (HANGLE) / 
     &                                                 DCOSD (ALTITUDE))
*
      IF ((((0.0 .LE. AZIMUTH) .AND. (AZIMUTH .LT. 90.0)) .AND.	! I quadrant
     &    ((270.0 .LE. HANGLE) .AND. (HANGLE .LT. 360.0))) .OR.	! IV quadrant 
     &  (((-90.0 .LE. AZIMUTH) .AND. (AZIMUTH .LT. 0.0)) .AND.	! IV quadrant
     &      ((0.0 .LE. HANGLE) .AND. (HANGLE .LT. 90.0)))) 	! I quadrant
     &   AZIMUTH = 180.0 - AZIMUTH			! correct the quadrant
      IF (AZIMUTH .LT. 0.0) AZIMUTH = AZIMUTH + 360.0
*
      RETURN
      END
