REAL *8 FUNCTION JD ( IDAY, IMON, IYEAR )
C
C           Returns the julian day for 0 hours ut.
      IMPLICIT NONE           
      SAVE
      INTEGER *4 IDAY, IMON, IYEAR, YEAR0, YEAR1, DOY
      PARAMETER (YEAR0=1986)
      PARAMETER (YEAR1=1995)
      REAL *8 JD0(YEAR0:YEAR1)
      EXTERNAL DOY
      DATA JD0 / 2446430.5D0,  2446795.5D0, 2447160.5, 2447526.5,
     $           2447891.5D0,  2448256.5D0, 2448621.5, 2448987.5,
     $           2449352.5D0,  2449717.5D0 /
      IF (( IYEAR .LT. YEAR0 ).OR.( IYEAR .GT. YEAR1 ) ) THEN
          WRITE(0,'(A,I5)') ' INVALID YEAR FOR COMPUTING JD:', IYEAR
      ELSE
          JD = JD0(IYEAR) + DOY( IDAY, IMON, IYEAR )
      END IF
      RETURN
      END
      INTEGER *4 FUNCTION DOY( IDAY, IMON, IYEAR )
C
C     Returns Day of Year for date in arguements ( 1 Jan => doy=1 )
      IMPLICIT NONE           
      SAVE
      INTEGER *4 IDAY, IMON, IYEAR, LEAP, I, NDAY(12,0:1)
      DATA NDAY / 0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30,
     $            0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30 /
C        LEAP=0 for leap lear
      LEAP = MIN( 1, MOD( IYEAR, 4 ) )
      DOY = IDAY
      DO 100 I = 1, IMON
          DOY = DOY + NDAY(I, LEAP)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE NEXTDAY ( IDAY, IMON, IYEAR )
C
C   Increments the date by one
C
      IMPLICIT NONE           
      INTEGER *4 IDAY, IMON, IYEAR, DAYS(12,0:1), LEAP
      DATA DAYS / 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31,
     $            31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /
      LEAP = MIN ( 1, MOD( IYEAR, 4 ) )
      IDAY = IDAY + 1
      IF ( IDAY .GT. DAYS(IMON,LEAP) ) THEN
          IMON = IMON + 1
          IDAY = 1
          IF ( IMON .EQ. 13 ) THEN
              IMON = 1
              IYEAR = IYEAR + 1
          END IF
      END IF
      RETURN
      END
      SUBROUTINE DIURN( TJD, GAST, RAA, DCA, RAT, DCT )
C
C           Calculates diurnal abberation
C    TJD   Julian date ( in Terestial dynamical time )
C    GAST  Greenwich aparent sideral time
C    RAA   Apparent geocentric right ascen.
C    DCA   Apparent geocentric decl.
C    RAT   Topocentric right ascen. ( output )
C    DCT   Topocentric decl. ( output )
C
      IMPLICIT NONE           
      SAVE
      REAL *8 PI, TWOPI, RADCON, C, ERAD, AU, VROT, DLAT, DLON, HT
      REAL *8 SINLAT, COSLAT, RHO, ROTVEL, ABRCON
      REAL *8 TJD, RAA, DCA, RAT, DCT, GAST, HA
      REAL *8 SINHA, COSHA, SINDEC, COSDEC, ALT, SINALT
      PARAMETER ( PI = 3.141592653589793D0 )
      PARAMETER ( TWOPI = 2.D0*PI )
      PARAMETER ( RADCON = PI / 180.D0 )
      PARAMETER ( C = 2.99792538D8 )
      PARAMETER ( ERAD = 6378140.D0 )
      PARAMETER ( AU = 1.49597870D11 )
      PARAMETER ( VROT = 7.2921151467D-5 )
C           Mt Wilson geocentric coordinates.
      PARAMETER ( DLAT =  34.2166944D0 )
      PARAMETER ( DLON = 118.0591670D0 )
      PARAMETER ( HT   = 1742.0D0 )
      LOGICAL FIRST
      DATA FIRST / .TRUE. /
C
      IF ( FIRST ) THEN 
          FIRST  = .FALSE.
          SINLAT = DSIN ( DLAT * RADCON )
          COSLAT = DCOS ( DLAT * RADCON )
          RHO = ( ERAD + HT  - 21400.D0 * SINLAT*SINLAT ) / ERAD
          ROTVEL = VROT * RHO * ERAD
          ABRCON = ROTVEL * COSLAT / C
      END IF
      HA  = ( (GAST - RAA ) * 15.D0 - DLON ) * RADCON
      SINHA  = DSIN( HA )
      COSHA  = DCOS( HA )
      SINDEC = DSIN( DCA * RADCON )
      COSDEC = DCOS( DCA * RADCON )
C
C         Test for object above the hroizon
C
      SINALT = SINLAT * SINDEC + COSLAT * COSDEC * COSHA
      ALT    = DASIN ( SINALT ) / RADCON
      IF ( ALT .LE. -1.D0 ) THEN
          RAT = 99.9D0
          DCT = 99.9D0
      ELSE
C           Apply correction for dirunal abberation
          RAT = RAA + ( ABRCON * COSHA / COSDEC ) /(15.D0*RADCON)
          DCT = DCA + ( ABRCON * SINHA * SINDEC ) / RADCON 
      END IF
      RETURN
      END