CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C objectdist : program to compute the angular distance between [RA,DEC]
C              coordinates.
C
C By Don Jennings, GSFC/HSTX, 04/14/95
C
C USAGE: objectdist MAXDIST RA1 DEC1 RA2 DEC2 .... RAn DECn
C
C        MAXDIST, RA*, DEC* in decimal degrees
C
C OUTPUT: angular distance between each pair of [RA,DEC] values in decimal
C         degrees and a warning if any distance value is greater than MAXDIST
C
C The formula for calculating the angular difference between ra,dec positions
C was taken from "Computational Spherical Astronomy" by Laurence G. Taff,
C page 16.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      program objectdist
C
      implicit     none
      integer      MAX
      parameter    (MAX=20)
      real         VER
      parameter    (VER=1.0)
C
      real*8       ra(MAX),dec(MAX),factor1,factor2,factor3,
     *             distance,maxdist
      integer      iargc,num,rnum,diff,i,j,k
      character*20 arg
C
      k = 0
C 
      num  = iargc()
      rnum = num/2
      diff = num - rnum*2
      if(num .lt. 5 .or. diff .eq. 0)then
         print *,'USAGE: objdist MAXDIST RA1 DEC1 ... RAn DECn'
      else
         write(6,50)VER
         call getarg(1,arg)
         read(arg,*) maxdist
         write(6,40) maxdist
         do i = 2,num,2
            k = k + 1
            call getarg(i,arg)
            read(arg,*) ra(k)
            call getarg(i+1,arg)
            read(arg,*) dec(k)
            write(6,30),k,ra(k),dec(k)
         enddo
C
         do i = 1,k-1
            do j = i+1,k
               factor1 = dec(i)
               factor2 = dec(j)
               factor3 = ra(i) - ra(j)
               distance = dacosd(dsind(factor1)*dsind(factor2) + 
     *                  dcosd(factor1)*dcosd(factor2)*dcosd(factor3))
               write(6,10) i,j,distance
               if(distance .gt. maxdist)then
                  write(6,20) i,j,maxdist
               endif
            enddo
         enddo
C
      endif
C
 10   format('objdist: distance from ',I3,' to ',I3,' = ',F10.6,
     *' degrees')
 20   format('objdist: WARNING: distance from ',I3, ' to ',I3,
     *' more than ',F10.6,' degrees')
 30   format('objdist: point ',I3,': RA = ',F10.6,' deg, DEC = ',
     *       F10.6,' deg')
 40   format('objdist: maximum angular distance between points = ',
     *       F10.6,' degrees')
 50   format('OBJDIST VERSION ',F4.2)
C
      end