w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-08 10:06:11 -05:00
Code Issues Packages Projects Releases Wiki Activity
5a8fa2f379
WSJT-X/lib/jplsubs.f

900 lines
24 KiB
FortranFixed
Raw Normal View History

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 12:43:40 -04:00
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 Copy Permalink