mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2024-08-17 20:16:35 -04:00
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6122 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
900 lines
24 KiB
Fortran
900 lines
24 KiB
Fortran
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*256 NAMFIL
|
|
|
|
c 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
|
|
|
|
c NRFILE=12
|
|
|
|
C *****************************************************************
|
|
|
|
C FIND THE RECORD SIZE USING THE INQUIRE STATEMENT
|
|
|
|
|
|
c 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
|
|
if(nrfile.eq.-99) stop !silence compiler warning
|
|
|
|
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
|
|
|
|
INTEGER OLDMAX
|
|
PARAMETER (OLDMAX = 400)
|
|
INTEGER NMAX
|
|
PARAMETER (NMAX = 1000)
|
|
CHARACTER*6 TTL(14,3),CNAM(NMAX)
|
|
CHARACTER*256 NAMFIL,jpleph_file_name
|
|
DIMENSION SS(3)
|
|
INTEGER IPT(3,13)
|
|
common/jplcom/jpleph_file_name
|
|
|
|
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'
|
|
NAMFIL=jpleph_file_name
|
|
|
|
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(K),K=1,OLDMAX),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)
|
|
|
|
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*256 NAMFIL,jpleph_file_name
|
|
common/jplcom/jpleph_file_name
|
|
|
|
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
|
|
|
|
NRECL=4
|
|
|
|
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'
|
|
NAMFIL=jpleph_file_name
|
|
|
|
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
|
|
C For de414, set KSIZE to 2036
|
|
C For de418, set KSIZE to 2036
|
|
C For de421, set KSIZE to 2036
|
|
C For de422, set KSIZE to 2036
|
|
C For de423, set KSIZE to 2036
|
|
C For de424, set KSIZE to 2036
|
|
C For de430, set KSIZE to 2036
|
|
|
|
KSIZE = 2036
|
|
|
|
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 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)
|
|
|
|
INTEGER NMAX
|
|
PARAMETER (NMAX = 1000)
|
|
|
|
DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13)
|
|
DIMENSION PVST(6,11),PNUT(4)
|
|
DIMENSION SS(3),CVAL(NMAX),PVSUN(6),ZIPS(2)
|
|
DATA ZIPS/2*0.d0/
|
|
|
|
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 IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
|
|
FIRST=.FALSE.
|
|
|
|
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,PVST,PNUT)
|
|
DO I=1,4
|
|
RRD(I)=PNUT(I)
|
|
ENDDO
|
|
RRD(5) = 0.d0
|
|
RRD(6) = 0.d0
|
|
RETURN
|
|
ELSE
|
|
DO I=1,4
|
|
RRD(I)=0.d0
|
|
ENDDO
|
|
WRITE(6,297)
|
|
297 FORMAT(' ***** NO NUTATIONS ON THE EPHEMERIS FILE *****')
|
|
STOP
|
|
ENDIF
|
|
|
|
C CHECK FOR LIBRATIONS
|
|
|
|
97 CONTINUE
|
|
DO I=1,6
|
|
RRD(I)=0.d0
|
|
ENDDO
|
|
|
|
IF(NTARG.NE.15) GO TO 98
|
|
IF(IPT(38).GT.0) THEN
|
|
LIST(12)=2
|
|
CALL STATE(ET2,LIST,PVST,PNUT)
|
|
DO I=1,6
|
|
RRD(I)=PVST(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
|
|
|
|
CALL STATE(ET2,LIST,PVST,PNUT)
|
|
|
|
DO I=1,10
|
|
DO J = 1,6
|
|
PV(J,I) = PVST(J,I)
|
|
ENDDO
|
|
ENDDO
|
|
|
|
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) = PVST(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)=PVST(I,3)-PVST(I,10)/(1.D0+EMRAT)
|
|
ENDDO
|
|
ENDIF
|
|
|
|
IF(LIST(10) .EQ. 2) THEN
|
|
DO I=1,6
|
|
PV(I,10) = PV(I,3)+PVST(I,10)
|
|
ENDDO
|
|
ENDIF
|
|
|
|
99 DO I=1,6
|
|
RRD(I)=PV(I,NTARG)-PV(I,NCENT)
|
|
ENDDO
|
|
|
|
BARY=BSAVE
|
|
|
|
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
|
|
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 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
|
|
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
|
|
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
|
|
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)
|
|
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
|
|
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
|
|
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)
|
|
6 CONTINUE
|
|
PV(I,2)=PV(I,2)*VFAC
|
|
5 CONTINUE
|
|
C
|
|
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 OUTPUT:
|
|
C
|
|
C PV DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED
|
|
C QUANTITIES (OTHER THAN NUTATION, STOERD IN PNUT).
|
|
C THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS
|
|
C STATE IN THE ARRAY STARTING AT PV(1,I).
|
|
C (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE
|
|
C AFFECTED BY THE FIRST 10 'LIST' ENTRIES, AND BY LIST(12)
|
|
C IF LIBRATIONS ARE ON THE FILE, ARE SET.
|
|
C THE REST OF THE 'PV' ARRAYIS UNTOUCHED.)
|
|
C THE ORDER OF COMPONENTS STARTING IN 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 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
|
|
|
|
INTEGER OLDMAX
|
|
PARAMETER ( OLDMAX = 400)
|
|
INTEGER NMAX
|
|
PARAMETER ( NMAX = 1000)
|
|
|
|
DIMENSION ET2(2),PV(6,11),PNUT(4),T(2),PJD(4),BUF(1500),
|
|
. SS(3),CVAL(NMAX),PVSUN(6)
|
|
|
|
INTEGER LIST(12),IPT(3,13)
|
|
|
|
LOGICAL FIRST
|
|
DATA FIRST/.TRUE./
|
|
|
|
CHARACTER*6 TTL(14,3),CNAM(NMAX)
|
|
CHARACTER*256 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)
|
|
C CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
|
|
CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
|
|
|
|
IF(NRECL .EQ. 0) WRITE(*,*)' ***** FSIZER IS NOT WORKING *****'
|
|
|
|
C ************************************************************************
|
|
C ************************************************************************
|
|
|
|
IRECSZ=NRECL*KSIZE
|
|
NCOEFFS=KSIZE/2
|
|
|
|
OPEN(NRFILE,
|
|
* FILE=NAMFIL,
|
|
* ACCESS='DIRECT',
|
|
* FORM='UNFORMATTED',
|
|
* RECL=IRECSZ,
|
|
* STATUS='OLD')
|
|
|
|
READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
|
|
& ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
|
|
& ,(CNAM(L),L=OLDMAX+1,NCON)
|
|
|
|
IF(NCON .LE. OLDMAX)THEN
|
|
READ(NRFILE,REC=2)(CVAL(I),I=1,OLDMAX)
|
|
ELSE
|
|
READ(NRFILE,REC=2)(CVAL(I),I=1,NCON)
|
|
ENDIF
|
|
|
|
NRL=0
|
|
|
|
ENDIF
|
|
|
|
C ********** MAIN ENTRY POINT **********
|
|
|
|
IF(ET2(1) .EQ. 0.D0) RETURN
|
|
|
|
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 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
|
|
|
|
tmp1 = DBLE(NR-3)*SS(3) + SS(1)
|
|
tmp2 = PJD(1) - tmp1
|
|
T(1) = (tmp2 + PJD(4))/SS(3)
|
|
|
|
C READ CORRECT RECORD IF NOT IN CORE
|
|
|
|
IF(NR.NE.NRL) THEN
|
|
NRL=NR
|
|
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
|
|
|
|
C INTERPOLATE SSBARY SUN
|
|
|
|
CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)
|
|
|
|
DO I=1,6
|
|
PVSUN(I)=PVSUN(I)*AUFAC
|
|
ENDDO
|
|
|
|
C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED
|
|
|
|
DO 4 I=1,10
|
|
IF(LIST(I).EQ.0) GO TO 4
|
|
|
|
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)
|
|
ELSE
|
|
PV(J,I)=PV(J,I)*AUFAC
|
|
ENDIF
|
|
ENDDO
|
|
|
|
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))
|
|
|
|
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
|
|
|
|
INTEGER NMAX
|
|
PARAMETER (NMAX = 1000)
|
|
|
|
CHARACTER*6 NAM(*),TTL(14,3),CNAM(NMAX)
|
|
|
|
DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(NMAX),ZIPS(2)
|
|
DOUBLE PRECISION PVST(6,11),PNUT(4)
|
|
DATA ZIPS/2*0.d0/
|
|
|
|
INTEGER IPT(3,13),DENUM,LIST(12)
|
|
logical first
|
|
data first/.true./
|
|
|
|
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
|
|
|
|
IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
|
|
first=.false.
|
|
|
|
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
|