      SUBROUTINE RISE_SET(JD, P, LAT, LON, RISE, TRANSIT, SET, STATUS)
C
C     GIVEN A JULIAN DATE (DAY NUMBER AND FRACTION), SOLAR SYSTEM
C     BODY NUMBER P, AND GEOGRAPHIC LATITUDE AND LONGITUDE (DEG),
C     THIS ROUTINE RETURNS THE JD OF THE TRANSIT NEAREST IN TIME TO
C     THE GIVEN JULIAN DATE, AND THE JD OF THE RISE AND SET WHICH
C     IMMEDIATELY PRECEDE AND FOLLOW THIS TRANSIT.  (ALGORITHM OF
C     J. MEEUS, ASTRONOMICAL ALGORITHMS, WILLMANN-BELL, 1991.)
C
C     B. KNAPP, 1998-06-25, 2001-04-24 (i*4),
C               2004-03-25 (EXTERNAL COSD, SIND, ACOSD, ASIND)
C
C     BODY NUMBERS: 0=SUN, 1=MERCURY, 2=VENUS, 4=MARS, 5=JUPITER,
C                   6=SATURN, 7=URANUS, 8=NEPTUNE, 9=PLUTO, 10=MOON
C
C     THE OUTPUT VARIABLE STATUS INDICATES SUCCESS/FAILURE, AND ALSO
C     WHETHER OR NOT THE REQUESTED OBJECT IS CIRCUM-POLAR (I.E., DOES
C     NOT SET) AT THE REQUESTED TIME AND GEOGRAPHIC POSITION:
C
C        STATUS = 0 ==> SUCCESS, OBJECT NOT CIRCUM-POLAR
C        STATUS = 1 ==> SUCCESS, OBJECT CIRCUM-POLAR
C        STATUS = 2 ==> FAILURE COMPUTING OBJECT'S RA AND DEC
C        STATUS = 3 ==> FAILURE (INVALID BODY NUMBER)
C        STATUS = 4 ==> FAILURE (NO TRANSIT IN JD +/- ND/2 DAYS)
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 JD, LAT, LON
      INTEGER*4 P
C
C     OUTPUT
      REAL*8 RISE, TRANSIT, SET
      INTEGER*4 STATUS
C
C     LOCAL
      INTEGER*4 ND
      PARAMETER (ND=6)
      REAL*8 DOFFSET
      PARAMETER (DOFFSET=ND/2.D0+0.5D0)
      REAL*8 ST0(ND), RA(ND), DEC(ND), DT(ND), DIST(ND), ID, DJ, JDE,
     &   YR8, DM, MINSEP, ST, M, N, A_RA, B_RA, C_RA, A_DEC, B_DEC,
     &   C_DEC, R, D, MT, MR, MS, COSLAT, SINLAT, COSH0, H0, H, ALT,
     &   HZPX, ALTP, PV(6)
      INTEGER*4 J, K
C
C     CONSTANTS
      REAL*8 SDAY, RE, AU
      PARAMETER (SDAY = 360.985647D0, RE=6378.135D0, AU=149597870.66D0)
      REAL*8 STDALT(0:10)/-0.83333333D0, -0.56666667D0, -0.56666667D0,
     &      0.00000000D0, -0.56666667D0, -0.56666667D0, -0.56666667D0,
     &     -0.56666667D0, -0.56666667D0, -0.56666667D0,  0.12500000D0/
      INTEGER*4 APPARENT, EQUATORIAL, SPHERICAL
      PARAMETER (APPARENT=2, EQUATORIAL=1, SPHERICAL=1)
C
C     EXTERNAL FUNCTIONS, SUBROUTINES
      REAL*8 COSD, SIND, ACOSD, ASIND, DELTAT, SIDTIM, TD2UT
      EXTERNAL COSD, SIND, ACOSD, ASIND, DELTAT, EPHEM, SIDTIM, TD2UT
C
C     CHECK BODY NUMBER
      IF (P .LT. 0 .OR. P .GT. 10 .OR. P .EQ. 3) THEN
         STATUS = 3
         RETURN
      ENDIF
C
C     FIND THE TRANSIT NEAREST IN TIME TO JD
      ID = DNINT( JD )
D     WRITE(*,10) JD, ID
D10   FORMAT(/' JD =',F11.1,'  ID =',F11.1/)
D     WRITE(*,20)
D    &   'J','K','DJ','RA(J)','DEC(J)','ST0(J)','M','DJ+M','MINSEP'
D20   FORMAT(2A3,3X,A2,10X,A5,7X,A6,5X,A6,7X,A1,9X,A4,12X,A6)
      MINSEP = 8.D0
      DO J=1,ND
         DJ = ID+J-DOFFSET
         ST0(J) = SIDTIM( DJ )
         JDE = TD2UT( DJ )
         YR8 = 2000.0D0+(DJ-2451545.0D0)/365.25D0
         DT(J) = DELTAT( YR8 )
         CALL EPHEM( JDE, P, APPARENT, EQUATORIAL, SPHERICAL, 
     &      PV, STATUS )
         IF ( STATUS .NE. 0 ) THEN
            STATUS = 2
            RETURN
         ENDIF
         RA(J) = PV(1)
         DEC(J) = PV(2)
         DIST(J) = PV(3)
         M = (RA(J)+LON-ST0(J))/360.
         IF      ( M .LT. 0.D0 ) THEN
            M = M+1.D0
         ELSE IF ( M .GT. 1.D0 ) THEN
            M = M-1.D0
         ENDIF
         IF ( ABS( DJ+M-JD ) .LT. MINSEP ) THEN
            MINSEP = ABS( DJ+M-JD )
            MT = M
            K = J
         ENDIF
D        WRITE(*,'(2I3,F12.1,3F12.6,F11.6,F16.6,F10.6)')
D    &      J, K, DJ, RA(J), DEC(J), ST0(J), M, DJ+M, MINSEP
      ENDDO
C
      IF ( K .EQ. 1 .OR. K .EQ. ND ) THEN
         STATUS = 4
         RETURN
      ENDIF
      A_RA = RA(K)-RA(K-1)
      IF      ( A_RA .LT. -180.D0 ) THEN
         A_RA = A_RA+360.D0
      ELSE IF ( A_RA .GT.  180.D0 ) THEN
         A_RA = A_RA-360.D0
      ENDIF
      B_RA = RA(K+1)-RA(K)
      IF      ( B_RA .LT. -180.D0 ) THEN
         B_RA = B_RA+360.D0
      ELSE IF ( B_RA .GT.  180.D0 ) THEN
         B_RA = B_RA-360.D0
      ENDIF
      C_RA = B_RA-A_RA
      A_DEC = DEC(K)-DEC(K-1)
      B_DEC = DEC(K+1)-DEC(K)
      C_DEC = B_DEC-A_DEC
C
C     TRANSIT
      DJ = ID+K-DOFFSET
      ST = ST0(K)+SDAY*MT
      N = MT+DT(K)/86400.d0
      R = RA(K) + (N/2.D0)*(A_RA+B_RA+N*C_RA)
      H = MOD( ST-LON-R, 360.D0 )
      IF ( H .LT. -180.D0 ) THEN
         H = H+360.D0
      ELSE IF ( H .GT. 180.D0 ) THEN
         H = H-360.D0
      ENDIF
      DM = -H/360.D0
      TRANSIT = DJ+MT+DM
C
C     STANDARD ALTITUDE OF OBJECT
      IF ( P .NE. 10 ) THEN
         ALTP = STDALT(P)
      ELSE
         HZPX = ASIND( RE/(DIST(K)*AU) )
         ALTP = 0.7275*HZPX-0.56666667D0
      ENDIF
C     WRITE(*,*) K, DIST(K), HZPX, ALTP
      COSLAT = COSD( LAT )
      SINLAT = SIND( LAT )
      COSH0 = (SIND( ALTP )-SINLAT*SIND( DEC(K) ))/
     &        (COSLAT*COSD( DEC(K) ))
      IF ( ABS( COSH0 ) .GT. 1.D0 ) THEN
C
C        OBJECT IS CIRCUM-POLAR
         RISE = TRANSIT-0.5D0
         SET = TRANSIT+0.5D0
         STATUS = 1
         RETURN
      ENDIF
      H0 = ACOSD( COSH0 )
D     WRITE(*,30) COSH0, H0
D30   FORMAT( /' COS(H0) =',F12.8,'  H0 =',F12.6/ )
C
D     WRITE(*,40) 'K','DJ','ST','M','N','R','D','H','DM'
D40   FORMAT( A6,2X,A2,9X,A2,10X,A1,9X,A1,10X,A1,11X,A1,10X,A1,11X,A2)
D     WRITE(*,50) 'T', K, DJ, ST, MT, N, R, 0., H, DM
D50   FORMAT( A3, I3, F11.1, F12.6, 2F10.6, 3F12.6, F10.6 )
C
C     RISE
      MR = MT-H0/360.D0
      IF ( MR .LT. 0.D0 ) MR = MR+1.D0
      ST = ST0(K) + SDAY*MR
      N = MR+DT(K)/86400.d0
      R = RA(K) + (N/2.D0)*(A_RA+B_RA+N*C_RA)
      H = MOD( ST-LON-R, 360.D0 )
      IF ( H .LT. -180.D0 ) THEN
         H = H+360.D0
      ELSE IF ( H .GT. 180.D0 ) THEN
         H = H-360.D0
      ENDIF
      D = DEC(K) + (N/2.D0)*(A_DEC+B_DEC+N*C_DEC)
      ALT = ASIND( SINLAT*SIND( D ) + COSLAT*COSD( D )*COSD( H ) )
      DM =(ALT-ALTP)/(COSD( D )*COSLAT*SIND( H )*360.D0)
      RISE = DJ + MR + DM
D     WRITE(*,50) 'R', K, DJ, ST, MR, N, R, D, H, DM
C
C     SET
      MS = MT+H0/360.D0
      IF ( MS .GT. 1.D0 ) MS = MS-1.D0
      ST = ST0(K) + SDAY*MS
      N = MS+DT(K)/86400.d0
      R = RA(K) + (N/2.D0)*(A_RA+B_RA+N*C_RA)
      H = MOD( ST-LON-R, 360.D0 )
      IF ( H .LT. -180.D0 ) THEN
         H = H+360.D0
      ELSE IF ( H .GT. 180.D0 ) THEN
         H = H-360.D0
      ENDIF
      D = DEC(K) + (N/2.D0)*(A_DEC+B_DEC+N*C_DEC)
      ALT = ASIND( SINLAT*SIND( D ) + COSLAT*COSD( D )*COSD( H ) )
      DM = (ALT-ALTP)/(COSD( D )*COSLAT*SIND( H )*360.D0)
      SET = DJ + MS + DM
D     WRITE(*,50) 'S', K, DJ, ST, MS, N, R, D, H, DM
C
      STATUS = 0
      RETURN
      END
