real*8 function angdst( ra1, dec1, ra2, dec2 )
!
!     Return the angular distance (radians) between two points on
!     the sphere.
!
!     B. Knapp, 1996-01-29, 2001-01-03
!
!     Reference: J. Meeus, Astronomical Algorithms, Willman-Bell, 1991.
!
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      implicit none
!
!     Input:
      real*8 ra1, dec1, ra2, dec2
!
!     Constants:
      real*8 pi
      parameter (pi=3.14159265358979324d0)
!
!     Local:
      real*8 c1, c2, d, ra, dec
!
!     Functions:
      real*8 dhav, dahav
      external dhav, dahav
!
!
      c1 = cos( dec1 )
      c2 = cos( dec2 )
      d = dahav( dhav(dec1-dec2)+c1*c2*dhav(ra1-ra2) )
      if (d .gt. pi) then
!
!        Move one point to its polar opposite
         dec = -dec2
         ra = ra2+pi
         d = pi - dahav( dhav(dec1-dec)+c1*c2*dhav(ra1-ra) )
      endif
      angdst = d
      return
      end