C ***** JPL Planetary and Lunar Ephemerides ***** C C Version : July 8, 1997 C C Program TESTEPH C --------------- C C TESTEPH tests the JPL ephemeris reading and interpolating routine using C examples computed from the original ephemeris. C C TESTEPH contains the reading and interpolating subroutines that are of C eventual interest to the user. Once TESTEPH is working correctly, the C user can extract those subroutines and the installation process is complete. C C You must supply "testpo.XXX" to TESTEPH, via standard input. "testpo.XXX C is the specially formatted text file that contains the test cases for the C ephmeris, DEXXX. C C After the initial identifying text which is concluded by an "EOT" in C columns 1-3, the test file contains the following quantities: C C JPL Ephemeris Number C calendar date C Julian Ephemeris Date C target number (1-Mercury, ...,3-Earth, ,,,9-Pluto, 10-Moon, 11-Sun, C 12-Solar System Barycenter, 13-Earth-Moon Barycenter C 14-Nutations, 15-Librations) C center number (same codes as target number) C coordinate number (1-x, 2-y, ... 6-zdot) C coordinate [au, au/day] C C For each test case input, TESTEPH C C - computes the corresponding state from data contained C in DExxx, C C - compares the two sets, C C - writes an error message if the difference between C any of the state components is greater than 10**(-13). C C - writes state and difference information for every 10th C test case processed. C C C This program is written in standard Fortran-77. C C HOWEVER, there are two parts which are compiler dependent; both have C to do with opening and reading a direct-access file. They are dealt C with in the subroutine FSIZERi, i=1,3. (There are three versions of C this subroutine. C C 1) The parameter RECL in the OPEN statement is the number of units per C record. For some compilers, it is given in bytes; in some, it is given C in single precision words. In the subroutine FSIZER of TESTEPH, the C parameter NRECL must be set to 4 if RECL is given in bytes; NRECL must C be set to 1 if RECL is given in words. (If in doubt, use 4 for UNIX; C 1 for VAX and PC) C C 2) Also for the OPEN statement, the program needs to know the exact value C of RECL (number of single precision words times NRECL). Since this C varies from one JPL ephemeris to another, RECL must be determined somehow C and given to the OPEN statement. There are three methods, depending C upon the compiler. We have included three versions of the subroutine C FSIZER, one for each method. C C a) Use the INQUIRE statement to find the length of the records C automatically before opening the file. This works for VAX's; C not in UNIX. C C b) Open the file with an arbitrary value of RECL, read the first record, C and use the information on that record to determine the exact value C of RECL. Then, close the file and re-open it with the exact value. C This seems to work for UNIX compilers as long as the initial value of C RECL is less than the exact value but large enough to get the required C information from the first file. (For other compilers, this doesn't C work since you can open a file only with the exact value of RECL.) C C c) Hardwire the value of RECL. This number is NRECL*1652 for DE200, C NRECL*2036 for DE405, and NRECL*1456 for DE406. C C C CHARACTER*6 NAMS(400) CHARACTER*3 ALF3 DOUBLE PRECISION DEL DOUBLE PRECISION ET DOUBLE PRECISION R(6) DOUBLE PRECISION SS(3) DOUBLE PRECISION VALS(400) DOUBLE PRECISION XI INTEGER I INTEGER LINE/0/ INTEGER NVS,NTARG,NCTR,NCOORD INTEGER NPT/100/ C Write a fingerprint to the screen. WRITE(*,*) ' JPL TEST-EPHEMERIS program. ' // . ' Last modified July 1997.' C Print the ephemeris constants. CALL Const (NAMS, VALS, SS, NVS) WRITE (*,'(/3F14.2)') SS C WRITE (*,'(200(/2(A8,D24.16)))') (NAMS(I),VALS(I),I=1,NVS) C Skip the file header comments. 1 READ(*,'(a3)')ALF3 IF(ALF3 .NE. 'EOT') GO TO 1 WRITE(*,*) . ' line -- jed -- t# c# x# --- jpl value ---'// . ' --- user value -- -- difference --' WRITE(*,*) C Read a value from the test case; skip if not within the time-range C of the present version of the ephemeris 2 READ(*,'(15X,D10.1,3I3,D22.13)',END=9) . ET,NTARG,NCTR,NCOORD,XI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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 version, we use only the four necessary C arguments and do the testing outside of the subroutine. IF(ET .LT. SS(1)) GO TO 2 IF(ET .GT. SS(2)) GO TO 2 CALL Pleph ( ET, NTARG, NCTR, R ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DEL=DABS(R(NCOORD)-XI) if(ntarg .eq. 15 .and. ncoord .eq. 3) * del=del/(0.23d0*(et-2451545.d0)) LINE=LINE+1 IF ( MOD(LINE,npt) .EQ. 0 ) WRITE(*,200) . LINE,ET,NTARG,NCTR,NCOORD,XI,R(NCOORD),DEL 200 FORMAT(I6,F10.1,3I5,3F20.13) C Print out WARNING if difference greater than tolerance. C IF (DEL .GE. 1.D-13) WRITE(*,201) . LINE,ET,NTARG,NCTR,NCOORD,XI,R(NCOORD),DEL 201 FORMAT(/ '***** WARNING : next difference >= 1.D-13 *****'// . I6,F10.1,3I3,3F20.13/' ') GO TO 2 9 STOP END C++++++++++++++++++++++++ C SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL) C C++++++++++++++++++++++++ C C Version 1.0 uses the INQUIRE statement to find out the the record length C of the direct access file before opening it. This procedure is non-standard, C but seems to work for VAX machines. C C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL. C ***************************************************************** C ***************************************************************** C C THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER C C ***************************************************************** C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE CHARACTER*80 NAMFIL NAMFIL='JPLEPH' C ***************************************************************** C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES C (for a VAX, it is probably 1) C NRECL= 4 C ***************************************************************** C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE NRFILE=12 C ***************************************************************** C FIND THE RECORD SIZE USING THE INQUIRE STATEMENT IRECSZ=0 INQUIRE(FILE=NAMFIL,RECL=IRECSZ) C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0 IF(IRECSZ .LE. 0) write(*,*) . ' INQUIRE STATEMENT PROBABLY DID NOT WORK' KSIZE=IRECSZ/NRECL RETURN END C++++++++++++++++++++++++ C SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL) C C++++++++++++++++++++++++ C THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS C THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE C PRECISION WORDS IN A RECORD. C C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL. IMPLICIT DOUBLE PRECISION(A-H,O-Z) SAVE CHARACTER*6 TTL(14,3),CNAM(400) CHARACTER*80 NAMFIL DIMENSION SS(3) INTEGER IPT(3,13) C ***************************************************************** C ***************************************************************** C C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER C C ***************************************************************** C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES C (for UNIX, it is probably 4) C NRECL=4 C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE NRFILE=12 C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE NAMFIL='JPLEPH' C ***************************************************************** C ***************************************************************** C ** OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO C ** DETERMINE THE SIZE OF THE EPHEMERIS RECORD MRECL=NRECL*1000 OPEN(NRFILE, * FILE=NAMFIL, * ACCESS='DIRECT', * FORM='UNFORMATTED', * RECL=MRECL, * STATUS='OLD') 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) CLOSE(NRFILE) C FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS KMX = 0 KHI = 0 DO I = 1,13 IF (IPT(1,I) .GT. KMX) THEN KMX = IPT(1,I) KHI = I ENDIF ENDDO ND = 3 IF (KHI .EQ. 12) ND=2 KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1) write(*,'(A6,I6,A8,I6,A15,3I8,/,A5,I6,A10,I6)') . 'KHI = ',KHI,', KMX = ',KMX,', IPT(*,KHI) = ', . (IPT(K,KHI),K=1,3),'ND = ',ND,', KSIZE = ',KSIZE RETURN END C++++++++++++++++++++++++ C SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL) C C++++++++++++++++++++++++ C C THE SUBROUTINE SETS THE VALUES OF NRECL, KSIZE, NRFILE, AND NAMFIL. SAVE CHARACTER*80 NAMFIL C ***************************************************************** C ***************************************************************** C C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER C ***************************************************************** C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES C NRECL= C ***************************************************************** C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12) NRFILE=12 C ***************************************************************** C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE NAMFIL='JPLEPH' C ***************************************************************** C KSIZE must be set by the user according to the ephemeris to be read C For de200, set KSIZE to 1652 C For de405, set KSIZE to 2036 C For de406, set KSIZE to 1456 KSIZE = 1652 C ******************************************************************* RETURN END C++++++++++++++++++++++++++ C SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD ) 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 necessary C arguments and do the testing outside of the subroutine. C C C C THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS C AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG' C WITH RESPECT TO 'NCENT'. C C CALLING SEQUENCE PARAMETERS: C C ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION C IS WANTED. C C ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME ** C THE REASON FOR THIS OPTION IS DISCUSSED IN THE C SUBROUTINE STATE C C NTARG = INTEGER NUMBER OF 'TARGET' POINT. C C NCENT = INTEGER NUMBER OF CENTER POINT. C C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS: C 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 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 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13) DIMENSION SS(3),CVAL(400),PVSUN(6) LOGICAL BSAVE,KM,BARY LOGICAL FIRST DATA FIRST/.TRUE./ INTEGER LIST(12),IPT(39),DENUM COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT COMMON/STCOMX/KM,BARY,PVSUN C INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT C ET2(1)=ET ET2(2)=0.D0 GO TO 11 C ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT C (SEE THE DISCUSSION IN THE SUBROUTINE STATE) ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD) ET2(1)=ET2Z(1) ET2(2)=ET2Z(2) 11 DO I=1,6 RRD(I)=0.D0 ENDDO IF(FIRST) CALL STATE(0.D0,0,0D0,0D0) FIRST=.FALSE. write(*,'(A12,/,A6,2F26.12,/,A8,I2,A10,I2)') . 'DPLEPH ENTRY','ET2 = ',ET2(1),ET2(2), . 'NTARG = ',NTARG,', NCENT = ',NCENT 96 IF(NTARG .EQ. NCENT) RETURN DO I=1,12 LIST(I)=0 ENDDO 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 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 FORCE BARYCENTRIC OUTPUT BY 'STATE' 98 BSAVE=BARY BARY=.TRUE. 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 MAKE CALL TO STATE do i=1,13 do j=1,6 pv(j,i)=0 end do end do write(*,'(A13,/,A6,2F26.12,/,A7,12I2)') . 'CALLING STATE','ET2 = ',ET2(1),ET2(2), . 'LIST = ',(LIST(I),I=1,12) CALL STATE(ET2,LIST,PV,RRD) IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN DO I=1,6 PV(I,11)=PVSUN(I) ENDDO ENDIF IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN DO I=1,6 PV(I,12)=0.D0 ENDDO ENDIF IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN DO I=1,6 PV(I,13)=PV(I,3) ENDDO ENDIF 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 IF(LIST(3) .EQ. 2) THEN DO I=1,6 PV(I,3)=PV(I,3)-PV(I,10)/(1.D0+EMRAT) ENDDO ENDIF IF(LIST(10) .EQ. 2) THEN DO I=1,6 PV(I,10)=PV(I,3)+PV(I,10) ENDDO ENDIF 99 DO I=1,6 RRD(I)=PV(I,NTARG)-PV(I,NCENT) ENDDO BARY=BSAVE write(*,'(A18,2(/,3F26.12))') . 'DPLEPH EXIT, RRD =',RRD(1),RRD(2),RRD(3),RRD(4),RRD(5),RRD(6) RETURN END C+++++++++++++++++++++++++++++++++ C SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV) C C+++++++++++++++++++++++++++++++++ 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 C BUF 1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION C 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 C NCF # OF COEFFICIENTS PER COMPONENT C C NCM # OF COMPONENTS PER SET OF COEFFICIENTS C C NA # OF SETS OF COEFFICIENTS IN FULL ARRAY C (I.E., # OF SUB-INTERVALS IN FULL INTERVAL) C C IFL INTEGER FLAG: =1 FOR POSITIONS ONLY C =2 FOR POS AND VEL C C C OUTPUT: C C PV INTERPOLATED QUANTITIES REQUESTED. DIMENSION C EXPECTED IS PV(NCM,IFL), DP. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C SAVE C DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18) 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 write(*,'(A7,/,A4,2F26.12,/,A6,I2,A8,I2,A7,I2,A8,I2)') . 'INTERP:','T = ',T(1),T(2),'NCF = ',NCF, . ', NCM = ',NCM,', NA = ',NA,', IFL = ',IFL DNA=DBLE(NA) DT1=DINT(T(1)) TEMP=DNA*T(1) L=IDINT(TEMP-DT1)+1 write(*,'(A6,F26.12,A8,F26.12,/,A7,F26.12,A6,I6)') . 'DNA = ',DNA,', DT1 = ',DT1,'TEMP = ',TEMP,', L = ',L C TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1) TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0 write(*,'(A5,F26.12)') 'TC = ',TC 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.) IF(TC.NE.PC(2)) THEN write(*,*) 'UPDATING NP/NV/PC(2)/TWOT' 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'. C IF(NP.LT.NCF) THEN write(*,*) 'UPDATING PC' DO 1 I=NP+1,NCF PC(I)=TWOT*PC(I-1)-PC(I-2) write(*,'(A3,I2,A4,F26.12)') 'PC(',I,') = ',PC(I) 1 CONTINUE NP=NCF ENDIF C C INTERPOLATE TO GET POSITION FOR EACH COMPONENT C 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) write(*,'(A10,I2,A6,F26.12,/,A11,I2,A1,I2,A1,I2,A4,F26.12)') . 'INTERP PV(',I,',1) = ',PV(I,1), . ' BUF(',J,',',I,',',L,') = ',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. C VFAC=(DNA+DNA)/T(2) VC(3)=TWOT+TWOT write(*,'(A7,F26.12,A10,F26.12)') . 'VFAC = ',VFAC,', VC(3) = ',VC(3) 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) write(*,'(A10,I2,A4,F26.12)') . 'INTERP VC(',I,') = ',VC(I) write(*,'(A14,F26.12)') . ' TWOT = ',TWOT write(*,'(A10,I2,A4,F26.12)') . ' VC(',I-1,') = ',VC(I-1) write(*,'(A10,I2,A4,F26.12)') . ' PC(',I-1,') = ',PC(I-1) write(*,'(A10,I2,A4,F26.12)') . ' VC(',I-2,') = ',VC(I-2) 4 CONTINUE NV=NCF ENDIF C C INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT C 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) write(*,'(A10,I2,A6,F26.12,/,A11,I2,A1,I2,A1,I2,A4,F26.12)') . 'INTERP PV(',I,',2) = ',PV(I,2), . ' BUF(',J,',',I,',',L,') = ',BUF(J,I,L) 6 CONTINUE PV(I,2)=PV(I,2)*VFAC write(*,'(A10,I2,A6,F26.12)') . 'IFINAL PV(',I,',2) = ',PV(I,2) 5 CONTINUE C write(*,*) 'INTERP EXIT, PV =' write(*,'(3F26.12)') ((PV(K1,K2),K2=1,IFL),K1=1,NCM) RETURN C END C+++++++++++++++++++++++++ C SUBROUTINE SPLIT(TT,FR) C C+++++++++++++++++++++++++ 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 C TT = D.P. INPUT NUMBER C C FR = D.P. 2-WORD OUTPUT ARRAY. C FR(1) CONTAINS INTEGER PART C FR(2) CONTAINS FRACTIONAL PART C 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 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FR(2) C MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS FR(1)=DINT(TT) FR(2)=TT-FR(1) IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN C MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER FR(1)=FR(1)-1.D0 FR(2)=FR(2)+1.D0 RETURN END C++++++++++++++++++++++++++++++++ C SUBROUTINE STATE(ET2,LIST,PV,PNUT) C C++++++++++++++++++++++++++++++++ C C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE C C CALLING SEQUENCE PARAMETERS: C C INPUT: C 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 C A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE C ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0. C 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 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 C LIST 12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION C IS WANTED FOR EACH OF THE BODIES ON THE FILE. C C LIST(I)=0, NO INTERPOLATION FOR BODY I C =1, POSITION ONLY C =2, POSITION AND VELOCITY C C THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS: C 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 C OUTPUT: C 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 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 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 C LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF C LIST(12) IS 1 OR 2. C 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 C D PSI (NUTATION IN LONGITUDE) C D EPSILON (NUTATION IN OBLIQUITY) C D PSI DOT C D EPSILON DOT C C * STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF C RANGE OR I/O ERRORS. C C C COMMON AREA STCOMX: C 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 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 C PVSUN DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND C VELOCITY OF THE SUN. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION ET2(2),PV(6,12),PNUT(4),T(2),PJD(4),BUF(1500), . SS(3),CVAL(400),PVSUN(3,2) INTEGER LIST(12),IPT(3,13) LOGICAL FIRST DATA FIRST/.TRUE./ CHARACTER*6 TTL(14,3),CNAM(400) CHARACTER*80 NAMFIL LOGICAL KM,BARY COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT COMMON/CHRHDR/CNAM,TTL COMMON/STCOMX/KM,BARY,PVSUN C C ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE C IF(FIRST) THEN FIRST=.FALSE. C ************************************************************************ C ************************************************************************ C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1 C ************************************************************************ C CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL) CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL) C CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL) IF(NRECL .EQ. 0) WRITE(*,*)' ***** FSIZER IS NOT WORKING *****' C ************************************************************************ C ************************************************************************ IRECSZ=NRECL*KSIZE NCOEFFS=KSIZE/2 write(*,'(A8,I6,A12,I6,A11,I6)') 'KSIZE = ',KSIZE, . ', NCOEFFS = ',NCOEFFS,', IRECSZ = ',IRECSZ OPEN(NRFILE, * FILE=NAMFIL, * ACCESS='DIRECT', * FORM='UNFORMATTED', * RECL=IRECSZ, * STATUS='OLD') 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) READ(NRFILE,REC=2)CVAL NRL=0 ENDIF C ********** MAIN ENTRY POINT ********** IF(ET2(1) .EQ. 0.D0) RETURN write(*,*) 'STATE ENTERED' write(*,'(A6,2F26.12)') 'ET2 = ',ET2(1),ET2(2) 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) write(*,'(A6,2F26.12,/,A6,2F26.12)') 'PJD = ',PJD(1),PJD(2), . ' ',PJD(3),PJD(4) 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 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) write(*,'(A5,I6,A9,F26.12)') 'NR = ',NR,', T(1) = ',T(1) C READ CORRECT RECORD IF NOT IN CORE IF(NR.NE.NRL) THEN NRL=NR write(*,*) 'READING' READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS) ENDIF IF(KM) THEN T(2)=SS(3)*86400.D0 AUFAC=1.D0 ELSE T(2)=SS(3) AUFAC=1.D0/AU ENDIF write(*,'(A7,F26.12,A10,F26.12)') 'T(2) = ',T(2), . ', AUFAC = ',AUFAC C INTERPOLATE SSBARY SUN write(*,*) 'INTERP FOR SSBARY SUN (11)' CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN) DO I=1,6 PVSUN(I,1)=PVSUN(I,1)*AUFAC ENDDO write(*,'(A8,2(/,3F26.12))') 'PVSUN = ',(PVSUN(K,1),K=1,6) C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED DO 4 I=1,10 IF(LIST(I).EQ.0) GO TO 4 write(*,'(A11,I2)') 'INTERP FOR ',I CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I), & LIST(I),PV(1,I)) 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 write(*,'(A3,I2,A4,2(/,3F26.12))') 'PV(',I,') = ',PV(1,I), . PV(2,I),PV(3,I),PV(4,I),PV(5,I),PV(6,I) 4 CONTINUE C DO NUTATIONS IF REQUESTED (AND IF ON FILE) IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0) * CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12), * LIST(11),PNUT) C GET LIBRATIONS IF REQUESTED (AND IF ON FILE) IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0) * CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13), * LIST(12),PV(1,11)) write(*,*) 'STATE EXIT, PV =' write(*,'(3F26.12)') ((PV(K2,K1),K2=1,6),K1=1,12) RETURN 98 WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2) 198 format(' *** Requested JED,',f12.2, * ' not within ephemeris limits,',2f12.2,' ***') stop 99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE' STOP END C+++++++++++++++++++++++++++++ C SUBROUTINE CONST(NAM,VAL,SSS,N) C C+++++++++++++++++++++++++++++ C C THIS ENTRY OBTAINS THE CONSTANTS FROM THE EPHEMERIS FILE C C CALLING SEQEUNCE PARAMETERS (ALL OUTPUT): C C NAM = CHARACTER*6 ARRAY OF CONSTANT NAMES C C VAL = D.P. ARRAY OF VALUES OF CONSTANTS C C SSS = D.P. JD START, JD STOP, STEP OF EPHEMERIS C C N = INTEGER NUMBER OF ENTRIES IN 'NAM' AND 'VAL' ARRAYS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE CHARACTER*6 NAM(*),TTL(14,3),CNAM(400) DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(400) INTEGER IPT(3,13),DENUM COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT COMMON/CHRHDR/CNAM,TTL C CALL STATE TO INITIALIZE THE EPHEMERIS AND READ IN THE CONSTANTS CALL STATE(0.D0,0,0.D0,0.D0) N=NCON DO I=1,3 SSS(I)=SS(I) ENDDO DO I=1,N NAM(I)=CNAM(I) VAL(I)=CVAL(I) ENDDO RETURN END