C************************************************************************* SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD ) Implicit None C The following code is from JPL with only a few necessary modifications: C************************************************************************* C C NOTE : Over the years, different versions of PLEPH have had a fifth argument: C sometimes, an error return statement number; sometimes, a logical denoting C whether or not the requested date is covered by the ephemeris. We apologize C for this inconsistency; in this present version, we use only the four C necessary arguments and do the testing outside of the subroutine. C C THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS AND GIVES THE C POSITION AND VELOCITY OF THE POINT 'NTARG' WITH RESPECT TO 'NCENT'. C C CALLING SEQUENCE PARAMETERS: C C ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION IS WANTED. C C ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME ** C THE REASON FOR THIS OPTION IS DISCUSSED IN THE SUBROUTINE STATE C C NTARG = INTEGER NUMBER OF 'TARGET' POINT. C NCENT = INTEGER NUMBER OF CENTER POINT. C C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS: C 1 = MERCURY 8 = NEPTUNE C 2 = VENUS 9 = PLUTO C 3 = EARTH 10 = MOON C 4 = MARS 11 = SUN C 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER C 6 = SATURN 13 = EARTH-MOON BARYCENTER C 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ) C 15 = LIBRATIONS, IF ON EPH FILE C (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS, C SET NTARG = 15. SET NCENT=0.) C C RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY C OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND C AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS C PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF C RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF C RADIANS AND RADIANS/DAY. C C The option is available to have the units in km and km/sec. C For this, set km=.true. in the STCOMX common block. C REAL*8 RRD(6),ET2Z(2),ET2(2),PV(6,13), ET, AU, EMRAT REAL*8 SS(3),CVAL(400),PVSUN(6) LOGICAL BSAVE,KM,BARY LOGICAL FIRST DATA FIRST/.TRUE./ INTEGER*4 LIST(12), IPT(39), DENUM, NTARG, NCENT, NCON, I, K COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT COMMON/STCOMX/KM,BARY,PVSUN C INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT ET2(1)=ET ET2(2)=0.D0 C GO TO 11 C C ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT C (SEE THE DISCUSSION IN THE SUBROUTINE STATE) C ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD) C ET2(1)=ET2Z(1) C ET2(2)=ET2Z(2) C 11 DO I=1,6 RRD(I)=0.D0 ENDDO C IF(FIRST) CALL STATE(0.D0,0,0,0) FIRST=.FALSE. C 96 IF(NTARG .EQ. NCENT) RETURN C DO I=1,12 LIST(I)=0 ENDDO C C CHECK FOR NUTATION CALL IF(NTARG.NE.14) GO TO 97 IF(IPT(35).GT.0) THEN LIST(11)=2 CALL STATE(ET2,LIST,PV,RRD) RETURN ELSE WRITE(6,297) 297 FORMAT(' ***** NO NUTATIONS ON THE EPHEMERIS FILE *****') STOP ENDIF C C CHECK FOR LIBRATIONS 97 IF(NTARG.NE.15) GO TO 98 IF(IPT(38).GT.0) THEN LIST(12)=2 CALL STATE(ET2,LIST,PV,RRD) DO I=1,6 RRD(I)=PV(I,11) ENDDO RETURN ELSE WRITE(6,298) 298 FORMAT(' ***** NO LIBRATIONS ON THE EPHEMERIS FILE *****') STOP ENDIF C C FORCE BARYCENTRIC OUTPUT BY 'STATE' 98 BSAVE=BARY BARY=.TRUE. C C SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL DO I=1,2 K=NTARG IF(I .EQ. 2) K=NCENT IF(K .LE. 10) LIST(K)=2 IF(K .EQ. 10) LIST(3)=2 IF(K .EQ. 3) LIST(10)=2 IF(K .EQ. 13) LIST(3)=2 ENDDO C C MAKE CALL TO STATE CALL STATE(ET2,LIST,PV,RRD) C IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN DO I=1,6 PV(I,11)=PVSUN(I) ENDDO ENDIF C IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN DO I=1,6 PV(I,12)=0.D0 ENDDO ENDIF C IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN DO I=1,6 PV(I,13)=PV(I,3) ENDDO ENDIF C IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN DO I=1,6 PV(I,3)=0.D0 ENDDO GO TO 99 ENDIF C IF(LIST(3) .EQ. 2) THEN DO I=1,6 PV(I,3)=PV(I,3)-PV(I,10)/(1.D0+EMRAT) ENDDO ENDIF C IF(LIST(10) .EQ. 2) THEN DO I=1,6 PV(I,10)=PV(I,3)+PV(I,10) ENDDO ENDIF C 99 DO I=1,6 RRD(I)=PV(I,NTARG)-PV(I,NCENT) ENDDO C BARY=BSAVE C RETURN END C C************************************************************************** SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV) Implicit None C C THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A C SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY C C CALLING SEQUENCE PARAMETERS: C C INPUT: C BUF 1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION C T T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY C COEFFICIENTS AT WHICH INTERPOLATION IS WANTED C (0 .LE. T(1) .LE. 1). T(2) IS DP LENGTH OF WHOLE C INTERVAL IN INPUT TIME UNITS. C NCF # OF COEFFICIENTS PER COMPONENT C NCM # OF COMPONENTS PER SET OF COEFFICIENTS C NA # OF SETS OF COEFFICIENTS IN FULL ARRAY C (I.E., # OF SUB-INTERVALS IN FULL INTERVAL) C IFL INTEGER FLAG: =1 FOR POSITIONS ONLY C =2 FOR POS AND VEL C C OUTPUT: C PV INTERPOLATED QUANTITIES REQUESTED. DIMENSION C EXPECTED IS PV(NCM,IFL), DP. C SAVE C WEW - move following line from after to before REAL*8 declaration INTEGER*4 NCF, NCM, NA, IFL, L, NP, NV, I, J REAL*8 BUF(NCF,NCM,*), T(2), PV(NCM,*), PC(18), VC(18), * DNA, DT1, TEMP, TC, TWOT, VFAC C INTEGER*4 NCF, NCM, NA, IFL, L, NP, NV, I, J C DATA NP/2/ DATA NV/3/ DATA TWOT/0.D0/ DATA PC(1),PC(2)/1.D0,0.D0/ DATA VC(2)/1.D0/ C C ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET C OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME C WITHIN THAT SUBINTERVAL. C DNA=DBLE(NA) DT1=DINT(T(1)) TEMP=DNA*T(1) L=IDINT(TEMP-DT1)+1 C TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1) TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0 C C CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED, C AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS. C (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE C CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.) C IF(TC.NE.PC(2)) THEN NP=2 NV=3 PC(2)=TC TWOT=TC+TC ENDIF C C BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED C AND ARE STORED IN THE ARRAY 'PC'. IF(NP.LT.NCF) THEN DO 1 I=NP+1,NCF PC(I)=TWOT*PC(I-1)-PC(I-2) 1 CONTINUE NP=NCF ENDIF C C INTERPOLATE TO GET POSITION FOR EACH COMPONENT DO 2 I=1,NCM PV(I,1)=0.D0 DO 3 J=NCF,1,-1 PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L) 3 CONTINUE 2 CONTINUE IF(IFL.LE.1) RETURN C C IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH C DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED. VFAC=(DNA+DNA)/T(2) VC(3)=TWOT+TWOT IF(NV.LT.NCF) THEN DO 4 I=NV+1,NCF VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2) 4 CONTINUE NV=NCF ENDIF C C INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT DO 5 I=1,NCM PV(I,2)=0.D0 DO 6 J=NCF,2,-1 PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L) 6 CONTINUE PV(I,2)=PV(I,2)*VFAC 5 CONTINUE C RETURN END C C*********************************************************************** SUBROUTINE SPLIT(TT,FR) Implicit None C C THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER C AND A D.P. FRACTIONAL PART. C C CALLING SEQUENCE PARAMETERS: C TT = D.P. INPUT NUMBER C FR = D.P. 2-WORD OUTPUT ARRAY. C FR(1) CONTAINS INTEGER PART. C FR(2) CONTAINS FRACTIONAL PART. C FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT C MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION. C C CALLING SEQUENCE DECLARATIONS REAL*8 TT, FR(2) C C MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS FR(1)=DINT(TT) FR(2)=TT-FR(1) C IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN C C MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER FR(1)=FR(1)-1.D0 FR(2)=FR(2)+1.D0 C RETURN END C C******************************************************************** SUBROUTINE STATE(ET2,LIST,PV,PNUT) Implicit None C C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE C C CALLING SEQUENCE PARAMETERS: C C INPUT: C ET2 DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION C IS WANTED. ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS C WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH. C A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE C ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0. C B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) = C THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION C EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY C ELAPSED BETWEEN ET2(1) AND EPOCH. C C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET C ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION, C AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH. C LIST 12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION C IS WANTED FOR EACH OF THE BODIES ON THE FILE. C LIST(I)=0, NO INTERPOLATION FOR BODY I C =1, POSITION ONLY C =2, POSITION AND VELOCITY C THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS: C I = 1: MERCURY C = 2: VENUS C = 3: EARTH-MOON BARYCENTER C = 4: MARS C = 5: JUPITER C = 6: SATURN C = 7: URANUS C = 8: NEPTUNE C = 9: PLUTO C =10: GEOCENTRIC MOON C =11: NUTATIONS IN LONGITUDE AND OBLIQUITY C =12: LUNAR LIBRATIONS (IF ON FILE) C C OUTPUT: C PV DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED C QUANTITIES. THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS C STATE IN THE ARRAY STARTING AT PV(1,I). (ON ANY GIVEN C CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE AFFECTED BY THE C FIRST 10 'LIST' ENTRIES (AND BY LIST(12) IF LIBRATIONS ARE C ON THE FILE) ARE SET. THE REST OF THE 'PV' ARRAY C IS UNTOUCHED.) THE ORDER OF COMPONENTS STARTING IN C PV(1,I) IS: X,Y,Z,DX,DY,DZ. C ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN C EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR C GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200. C THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES C ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC, C DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW). C LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF C LIST(12) IS 1 OR 2. C NUT DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES, C DEPENDING ON THE SETTING OF LIST(11). THE ORDER OF C QUANTITIES IN NUT IS: C D PSI (NUTATION IN LONGITUDE) C D EPSILON (NUTATION IN OBLIQUITY) C D PSI DOT C D EPSILON DOT C * STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF C RANGE OR I/O ERRORS. C C COMMON AREA STCOMX: C KM LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT C STATES. KM = .TRUE., KM AND KM/SEC C = .FALSE., AU AND AU/DAY C DEFAULT VALUE = .FALSE. (KM DETERMINES TIME UNIT C FOR NUTATIONS AND LIBRATIONS. ANGLE UNIT IS ALWAYS RADIANS.) C BARY LOGICAL FLAG DEFINING OUTPUT CENTER. C ONLY THE 9 PLANETS ARE AFFECTED. C BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER C = .FALSE. =\ CENTER IS SUN C DEFAULT VALUE = .FALSE. C PVSUN DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND C VELOCITY OF THE SUN. C SAVE REAL*8 ET2(2),PV(6,12),PNUT(4),T(2),PJD(4),BUF(1500), . SS(3),CVAL(400),PVSUN(6,1), AU, EMRAT, S, AUFAC C . SS(3),CVAL(400),PVSUN(3,2), AU, EMRAT, S, AUFAC C INTEGER*4 LIST(12),IPT(3,13), NUMDE, NCON, NRECL, KSIZE, NRFILE, INTEGER*4 LIST(12),IPT(3,13), NUMDE, NCON, NRECL, NRFILE, * IRECSZ, NCOEFFS, I, J, K, NRL, NR INTEGER*4 I2, I3 DATA I2 /2/ DATA I3 /3/ LOGICAL FIRST DATA FIRST/.TRUE./ CHARACTER*6 TTL(14,3),CNAM(400) C CHARACTER*80 NAMFIL LOGICAL KM,BARY C INCLUDE 'param.i' C Variables from: C 1. JPL_eph - Character string giving the complete path name of C the JPL DE403 ephemeris. C COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT COMMON/CHRHDR/CNAM,TTL COMMON/STCOMX/KM,BARY,PVSUN DATA KM/.TRUE./ C C ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE IF(FIRST) THEN FIRST=.FALSE. C NRECL=4 NRFILE=12 C NAMFIL= 'JPL.DE403' C KSIZE = 1796 C C IRECSZ=NRECL*KSIZE C NCOEFFS=KSIZE/2 NCOEFFS = IRECL/NRECL/2 C C OPEN(NRFILE, FILE=NAMFIL, ACCESS='DIRECT', FORM='UNFORMATTED', C * RECL=IRECSZ, STATUS='OLD') OPEN(NRFILE, FILE=JPL_eph, ACCESS='DIRECT', FORM='UNFORMATTED', * RECL=IRECL, STATUS='OLD') C READ(NRFILE,REC=1)TTL,CNAM,SS,NCON,AU,EMRAT, . ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3) C READ(NRFILE,REC=2)CVAL c wew: Add write(6,'(/" Opened JPL DE200 Ephemeris - MJD range: ", . 2f12.1/)' ) SS(1), SS(2) write(6,'("TTL: ", 3(14(a6)/))') TTL write(6,'("CNAM: ",40(10(a6,1x)/))') CNAM write(6,'("SS: ",3d25.16)') SS write(6,'("NCON, NUMDE: ",2i10)') NCON, NUMDE write(6,'("AU, EMRAT: ",2d25.16)') AU, EMRAT write(6,'("IPT: ",3(13i8,/))') IPT NRL=0 C ENDIF C C ********** MAIN ENTRY POINT ********** IF(ET2(1) .EQ. 0.D0) RETURN C S=ET2(1)-.5D0 CALL SPLIT(S,PJD(1)) CALL SPLIT(ET2(2),PJD(3)) PJD(1)=PJD(1)+PJD(3)+.5D0 PJD(2)=PJD(2)+PJD(4) CALL SPLIT(PJD(2),PJD(3)) PJD(1)=PJD(1)+PJD(3) C C ERROR RETURN FOR EPOCH OUT OF RANGE IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98 C C CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL NR=IDINT((PJD(1)-SS(1))/SS(3))+3 IF(PJD(1).EQ.SS(2)) NR=NR-1 T(1)=((PJD(1)-(DBLE(NR-3)*SS(3)+SS(1)))+PJD(4))/SS(3) C C READ CORRECT RECORD IF NOT IN CORE IF(NR.NE.NRL) THEN NRL=NR C print *, '!!!!!!!!!! STATE: NR = ', NR READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS) ENDIF C IF(KM) THEN T(2)=SS(3)*86400.D0 AUFAC=1.D0 ELSE T(2)=SS(3) AUFAC=1.D0/AU ENDIF C C INTERPOLATE SSBARY SUN C CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN) CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),I3,IPT(3,11),I2,PVSUN) DO I=1,6 PVSUN(I,1)=PVSUN(I,1)*AUFAC ENDDO C C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED C DO 4 I=1,10 IF(LIST(I).EQ.0) GO TO 4 C CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I), C & LIST(I),PV(1,I)) CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),I3,IPT(3,I), & LIST(I),PV(1,I)) C DO J=1,6 IF(I.LE.9 .AND. .NOT.BARY) THEN PV(J,I)=PV(J,I)*AUFAC-PVSUN(J,1) ELSE PV(J,I)=PV(J,I)*AUFAC ENDIF ENDDO C 4 CONTINUE C C DO NUTATIONS IF REQUESTED (AND IF ON FILE) IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0) THEN C CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12), C * LIST(11),PNUT) CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),I2,IPT(3,12), * LIST(11),PNUT) ENDIF C C GET LIBRATIONS IF REQUESTED (AND IF ON FILE) IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0) THEN C CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13), C * LIST(12),PV(1,11)) CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),I3,IPT(3,13), * LIST(12),PV(1,11)) ENDIF C RETURN C 98 WRITE(6,198)ET2(1)+ET2(2),SS(1),SS(2) 198 format(' *** Requested JED,',f12.2, * ' not within ephemeris limits,',2f12.2,' ***') C STOP C 99 WRITE(6,'(2F12.2," ERROR RETURN IN STATE")') ET2 STOP END