w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-23 20:58:55 -05:00
Code Issues Packages Projects Releases Wiki Activity
4f9cfafec9
WSJT-X/lib/jplsubs.f
Joe Taylor 43856b5eda EME Doppler calculations now use JPL Planetary Ephemeris. 
Also new, simplified routines for Doppler spread.
Beware! Not yet fully tested ...



git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@5496 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2015-06-02 16:43:40 +00:00

900 lines
24 KiB
Fortran
Raw Blame History

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
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
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*80 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*80 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
common/jplcom/jpleph_file_name
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.
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,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*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)
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
Reference in New Issue View Git Blame Copy Permalink