SUBROUTINE xtereadeph(lue,Jdutc,Frc,Rce,Rcs,Etut,Vce,ierr) c Reads and interpolates ephemeris file (JPL data) c This routine is called by subroutine barycen to bake barycenter corrections. c barycen ignors etut and vce. Rcs is used only in the general relativity c part of subroutine barycen. c File lue should have been opened by the calling routine. C call_var. type I/O description c lue (i) I lu of JPL2000 ephemeris file CP JDUTC I4 I julian day number of ephemeris lookup CP FRC R8 I day fraction, between -0.5 and 0.5 CP RCE(3) R8 O vector from SSBC to Earth (light-s) CP RCS(3) R8 O vector from SSBC to Sun (light-s) CP ETUT R8 O Time difference TDB-UTC (seconds) CP VCE(3) R8 O time derivative of RCE c ierr (i) O error status C C************************ FFORM VERSION 3.0 ********** 31-JAN-90 09:20 C CA author : Lloyd Rawley-STX date: 14-JUL-1989 CA update : GDWirth date: 02-AUG-1989 CA update : TMB date: 27-SEP-1990 16:10 ULTRIX version CU update : TMB date: 01-OCT-1990 15:47 CU update : TMB date: 30-OCT-1990 11:28 no header CU update : eal HSTX date FEB-1994 To read from FITS ephemeris C file CU update : Mike Stark date: 05-FEB-1996 Eliminate UTC->TDT correction. C XTE time is TDT already. c update : Brian K Elza date: 17-APR-1996 Fixed bug left over from c ROSAT version, and notified c proper ROSAT people. C status: Tested C C general description: CG Reads and interpolates ephemeris file (JPL data), returning CG positions of the Earth and Sun with respect to the solar system CG barycenter (SSBC) and the difference between Ephemeris Time (TDB) CG and Universal Time (UTC) at the same epoch. C C C*********************************************************************** C C variables meaning C C NDAYCH I*4 (parameter) number of days of ephemeris C to keep in memory C LUE I*4 (parameter) logical unit for reading C DAYSEC R*8 (parameter) ratio of one second to one day C RCE(3) R*8 vector from SSBC to earth (light-s) C RCS(3) R*8 vector from SSBC to sun (light-s) C VCE(3) R*8 time derivative of RCE C ETUT R*8 time difference TDB-UTC (seconds) C EARTH(3,4,NDAYCH) C R*8 earth position and its first three C derivatives for each day C SUN(3,NDAYCH) R*8 sun positions for each day in memory C TDBTDT(NDAYCH) C R*8 relativistic time correction (TDB-TDT) C TDBDOT(NDAYCH) C R*8 rate of change of TDBTDT C FRC R*8 day fraction (between -0.5 and 0.5) C DT R*8 temporary variable for Taylor series C DT2 R*8 temporary variable for Taylor series C DT3 R*8 temporary variable for Taylor series C NSET I*4 loop index C JDUTC I*4 Julian day number of ephemeris lookup C JD0 I*4 beginning JD of ephemeris file C JD1 I*4 ending JD of ephemeris file C JDCH0 I*4 first JD of ephemeris in memory C JDCH1 I*4 last JD of ephemeris in memory C NREC I*4 pointer to current record in ephemeris file C I I*4 loop index C HEADER C*120 header of ephemeris file C C IMPLICIT NONE INTEGER*4 NDAYCH PARAMETER (NDAYCH=16) REAL*8 DAYSEC PARAMETER (DAYSEC=1.D0/86400.D0) REAL*8 ETATC PARAMETER (ETATC=32.184D0-0.0343817D0) C REAL*8 Rce(3) REAL*8 Rcs(3) REAL*8 Vce(3) REAL*8 earth(3,4,NDAYCH) REAL*8 sun(3,NDAYCH) REAL*8 tdbtdt(NDAYCH) REAL*8 tdbdot(NDAYCH) REAL*8 Etut REAL*8 Frc REAL*8 dt REAL*8 dt2 REAL*8 dt3 INTEGER*4 lue INTEGER*4 nset INTEGER*4 Jdutc INTEGER*4 jd0 INTEGER*4 jd1 INTEGER*4 jdch0 INTEGER*4 jdch1 INTEGER*4 nrec INTEGER*4 i INTEGER*4 ierr character(80) comm,errm double precision nulvd, dbuf data nulvd /-1.2d34/ logical anynul data anynul /.false./ DATA jdch0/0/ SAVE C----------------------------------------------------------------------- C Choose set of coefficents to use; if out of range, read new chunk. C----------------------------------------------------------------------- ierr = 0 nset = 1 + (Jdutc-jdch0) IF ( nset.GT.NDAYCH .OR. nset.LT.1 ) THEN C--------------------------------------------------------- C If this is the first call, read the header record C--------------------------------------------------------- IF ( jdch0.EQ.0 ) THEN CALL ftgkyd(lue,'TSTART',dbuf,comm,ierr) jd0 = int(dbuf) CALL ftgkyd(lue,'TSTOP' ,dbuf,comm,ierr) jd1 = int(dbuf) IF ( jd1-jd0.LT.(NDAYCH-1) ) THEN ierr = 2 errm = 'File is too small for interpolation' GOTO 999 ENDIF ENDIF C------------------------------------------------- C Check if date is in range of ephemeris C------------------------------------------------- IF ( (Jdutc.LT.jd0) .OR. (Jdutc.GT.jd1) ) THEN ierr = 2 errm = 'Out of range of ephemeris file' GOTO 999 ELSE C----------------------------------------------------------------------- C Read earth and sun positions and TDB-TDT from the ephemeris. C Get AT-UT from the leap second table in A1UTCF. C----------------------------------------------------------------------- jdch1 = MIN(Jdutc+(NDAYCH-1),jd1) jdch0 = jdch1 - (NDAYCH-1) DO nset = 1 , NDAYCH nrec = nset + (jdch0-jd0) CALL ftgcvd(lue,1,nrec,1,12,nulvd,earth(1,1,nset) & ,anynul,ierr) CALL ftgcvd(lue,2,nrec,1, 3,nulvd,sun(1,nset) & ,anynul,ierr) CALL ftgcvd(lue,3,nrec,1, 1,nulvd,tdbtdt(nset) & ,anynul,ierr) IF(ierr.ne.0) THEN errm = 'Reading ephemeris file' GOTO 999 ENDIF ENDDO C------------------------------------------------------------------- C Make rough computation of time derivative of TDB-TDT. C------------------------------------------------------------------- DO nset = 1 , NDAYCH - 1 tdbdot(nset) = tdbtdt(nset+1) - tdbtdt(nset) ENDDO tdbdot(NDAYCH) = tdbdot(NDAYCH-1) nset = 1 + (Jdutc-jdch0) ENDIF ENDIF C---------------------------------------------------------------------- C Now use Taylor series to get the coordinates at the desired time. C---------------------------------------------------------------------- dt = Frc dt2 = dt*dt dt3 = dt2*dt Etut = tdbtdt(nset) + dt*tdbdot(nset) DO i = 1 , 3 Rcs(i) = sun(i,nset) Rce(i) = earth(i,1,nset) + dt*earth(i,2,nset) & + 0.5D0*dt2*earth(i,3,nset) + (1.D0/6.D0) & *dt3*earth(i,4,nset) Vce(i) = (earth(i,2,nset) + dt * earth(i,3,nset) + X DT2 / 2.0D0 * Earth(I, 4, NSet)) * DAYSEC ENDDO RETURN 999 continue C CALL xaerror(errm,1) return END