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