      SUBROUTINE GETSEDR(EDRTIME,TD,TN,TP,ESPV,ERANGE,EANG,
     &                   MAT1,MAT2,MAT3)
C
C EXTRACT POINTING VECTOR AND NAVIGATION WORDS NECESSARY FOR THE COMPUTATION
C OF MATRICES TO ROTATE FIELD DATA FROM PAYLOAD TO IHG AND HG COORDINATES.
C
C ORIGINAL CODE BY SANDY KRAMER, HUGHES STX, CODE 692 NASA GSFC, 10/25/93
C HANDLE ENCOUNTER MODES, 05/01/95 - SBK
C INTERPOLATE POINTING VECTOR AND NAVIGATION DATA, 07/31/95 - SBK
C DETECT HEADER RECORDS IN CONCATENATED NAV FILES, 08/05/96 - SBK
C DETECT HEADER RECORDS IN CONCATENATED PTG FILES, 08/06/96 - SBK
C DETECT BAD NAV REC TIME TAGS, 08/22/96 - SBK
C RESTORE DETECTION OF FILE COMMENTS, 09/11/96 - SBK
C
      CHARACTER IHDR(45)*4,TYPE*4
      CHARACTER STAR*12,PVRD*77
      INTEGER*2 EDRTIME(6),PTIME(6),NTIME(6),LAUNCH(6)
      INTEGER*4 INAV(126,2)
      REAL*4 RNAV(126,2),NAV(252),PV(3,3),OPV(3,3),
     &       RH(3),VH(3),SPV(6),RANGE,ANG(2),OANG(2),EANG(2),
     &       OSPV(6),ESPV(6),MTB(3,3),MTB5(3,3),MHG(3,3),
     &       OMTB(3,3),OMTB5(3,3),OMHG(3,3),
     &       MAT1(3,3),MAT2(3,3),MAT3(3,3),PMAT(3,3)
      REAL*8 TD,TP,TN,REALTIME,DATATIME,
     &       PVTIME,NAVTIME,OTN,OTP
C
      DATA ICALL/0/,ARAD/1.495985E8/
C
      EQUIVALENCE ( IHDR(1), INAV(1,1) ),
     &            ( NAV(1), RNAV(1,1) )
C
      SAVE PVRD,PVTIME,PV,NAVTIME,INAV,NAV,ICALL
C
C      INCLUDE 'UNPACK.INC'
C
C DIAGNOSTIC OUTPUT UNIT
C
      IUNIT = 68
C
C COUNT CALLS FOR SEDR DATA
C
      ICALL = ICALL + 1
C
C LAUNCH TIME FOR CALCULATION OF NAV AND PNTG EPIC TIMES
C
      LAUNCH(1) = 77
      LAUNCH(2) = 232
      LAUNCH(3) = 0
      LAUNCH(4) = 0
      LAUNCH(5) = 0
      LAUNCH(6) = 0
C
C DETERMINE DESIRED COORDINATE PROCESSING
C
c      IF ( SYS2(17) ) THEN
c       ICOORD = 4  ! JUPITER
c      ELSE IF ( SYS2(18) ) THEN
c       ICOORD = 3  ! SATURN
c      ELSE IF ( SYS2(19) ) THEN
c       ICOORD = 5  ! URANUS
c      ELSE IF ( SYS2(20) ) THEN
c       ICOORD = 6  ! NEPTUNE
c      ELSE
c       ICOORD = 2  ! HG (DEFAULT)
c      END IF
c
      icoord =2
C
C CONVERT CALENDAR DATA TIME TO DECIMAL YEAR.
C
      DATATIME = REALTIME(EDRTIME)
C
C COMPUTE ELAPSED EDR TIME
C
      CALL ELPSTIME(LAUNCH,EDRTIME,TD)
C
      IF ( ICALL.GT.1 .AND. DATATIME.GT.PVTIME ) GOTO 7  ! GET NEW PTG VECTOR
      IF ( ICALL.GT.1 .AND. DATATIME.LE.PVTIME ) GOTO 12 ! CHECK NAV DATA
C
C READ POINTING VECTOR FILE HEADER.  EXPECT $$VGR AT HEADER START AND
C $$EOH AT HEADER END.
C
      READ(40,'(A77)') PVRD
    4 CONTINUE
      IF ( PVRD(1:5).NE.'$$VGR' ) THEN 
       WRITE(IUNIT,804)
       STOP
      END IF
      WRITE(IUNIT,*)
      WRITE(IUNIT,*) 'POINTING VECTOR FILE INFO'
      DO WHILE ( PVRD(1:5).NE.'$$EOH' )
       READ(40,'(A77)') PVRD
       IF ( PVRD(1:1).NE.'*' ) WRITE(IUNIT,*) PVRD
      END DO
C
C READ POINTING VECTOR DATA
C
    5 CONTINUE
      write(6,*) 'reading pointing vector file...'
      READ(40,'(A77)',END=200) PVRD
      IF (PVRD(1:5).EQ.'$$VGR') GOTO 4  ! NEW HEADER BLOCK
      IF (PVRD(1:1).EQ.'*') THEN
       WRITE(IUNIT,*) 
       WRITE(IUNIT,*) 'POINTING VECTOR FILE INFO'
       WRITE(IUNIT,'(A77)') PVRD 
       GOTO 5
      END IF
C
C GET POINTING VECTORS' TIME TAG
C
      CALL EXTRACTDATE(PTIME,PVRD)
      PVTIME = REALTIME(PTIME)
C
C COMPUTE ELAPSED POINTING VECTOR TIME
C
      CALL ELPSTIME(LAUNCH,PTIME,TP)
C
      IF (PVTIME.LT.DATATIME) THEN
C
C LOAD 3X3 SENSOR POINTING MATRIX.  INSTRUMENT AXIS VECTORS ARE COLUMN
C ORIENTED.
C
C X AXIS VECTOR TO COLUMN 1
C
       READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,1),I=1,3)
C
C Y AXIS VECTOR TO COLUMN 2
C
       READ(40,'(A77)',END=200) PVRD
       READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,2),I=1,3)
C
C Z AXIS VECTOR TO COLUMN 3
C
       READ(40,'(A77)',END=200) PVRD
       READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,3),I=1,3)
C
      END IF
    7 CONTINUE
      IF (PVTIME.LT.DATATIME) THEN
C
C SAVE POINTING VECTOR DATA BEFORE GETTING NEXT RECORD
C
       DO J = 1,3
        DO I = 1,3
         OPV(I,J) = PV(I,J)
        END DO
       END DO
       OTP = TP
       GOTO 5
      END IF
C
C GET POINTING VECTOR WITH TIME TAG FOLLOWING DATA TIME TAG
C
C LOAD 3X3 SENSOR POINTING MATRIX.  INSTRUMENT AXIS VECTORS ARE COLUMN
C ORIENTED.
C
C X AXIS VECTOR TO COLUMN 1
C
      READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,1),I=1,3)
C
C Y AXIS VECTOR TO COLUMN 2
C
      READ(40,'(A77)',END=200) PVRD
      READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,2),I=1,3)
C
C Z AXIS VECTOR TO COLUMN 3
C
      READ(40,'(A77)',END=200) PVRD
      READ(PVRD,'(34X,2(E13.12,2X),E13.12)') (PV(I,3),I=1,3)
C
C INTERPOLATE POINTING VECTOR DATA
C
      CALL SEDRIP(OTP,OPV,TP,PV,TD,PMAT,9,0)
C
      IF (ICALL.GT.1) GOTO 12
C
C GET NAVIGATION FILE HEADER RECORD.  CONTAINS CHARACTER AND INTEGER DATA.
C
      READ(41) IHDR
      IF ( IHDR(2).NE.'SEDR' ) THEN
       WRITE(IUNIT,805)
       STOP
      END IF
    9 CONTINUE
      WRITE(IUNIT,*)
      WRITE(IUNIT,*) 'NAVIGATION FILE INFO'
      WRITE(IUNIT,'(1X,''PROJECT ID:  '',A4)') IHDR(1) 
      WRITE(IUNIT,'(1X,''FILE TYPE:   '',A4)') IHDR(2)
      WRITE(IUNIT,'(1X,''FLIGHT:      '',I2.2)') INAV(3,1)
      WRITE(IUNIT,'(1X,''FILE ID:     '',2A4)') IHDR(4),IHDR(5)
      WRITE(IUNIT,'(1X,''FILE FORMAT: '',2A4)') IHDR(14),IHDR(15)
      TYPE = IHDR(14)
C
C GET NAVIGATION DATA RECORDS.  CONTAINS INTEGER AND FLOATING POINT DATA.
C FLOATING POINT VALUES ARE CONVERTED FROM VAXG TO IEEE REPRESENTATION BY
C DEC FORTRAN CONVERT PARAMETER IN OPEN STATEMENT WHEN SOURCE CODE IS 
C COMPILED USING THE FLOAT=IEEE FLAG.
C
   10 CONTINUE
       write(6,*) 'reading navigation data file...'
       READ(41,END=100,ERR=11) (INAV(I,1),I=1,6),(RNAV(I,1),I=7,126),
     &                         (RNAV(I,2),I=1,126)
   11  CONTINUE
C
C GET NAVIGATION RECORD TIME TAG.  USE 2 DIGIT YEAR.
C
       IF ( IHDR(2).EQ.'SEDR' ) GOTO 9  ! DETECT HDR REC IN CONCAT SEDRS
       NTIME(1) = INAV(1,1) - 1900
       IF (NTIME(1).GT.99) NTIME(1) = NTIME(1) - 100
       NTIME(2) = INAV(2,1)
       NTIME(3) = INAV(3,1)
       NTIME(4) = INAV(4,1)
       NTIME(5) = INAV(5,1)
       NTIME(6) = INAV(6,1)
       NAVTIME = REALTIME(NTIME)
C
C COMPUTE ELAPSED NAVIGATION TIME
C
       CALL ELPSTIME(LAUNCH,NTIME,TN)
C
C DETECT BAD NAV REC BASED ON VALID TIME TAG
C
       IF ( TN.EQ.-1.0D0 ) GOTO 10
C
C STORE CONTENTS OF NAVIGATION RECORD WHILE NAVIGATION TIME IS LESS THAN
C EDR TIME
C
       IF (NAVTIME.LT.DATATIME) THEN
C
C CALL SEDR PROCESSING ROUTINE 
C
        IF ( ICOORD.EQ.2 ) THEN
         CALL SEDRCRU(NTIME,NAV,PMAT,SPV,ANG,MTB,MHG,MTB5,RANGE)
        ELSE IF ( ICOORD.GE.3 ) THEN
         CALL SEDRENC(NTIME,NAV,PMAT,SPV,ANG,MTB,MHG,MTB5,RANGE)
        END IF
C
       END IF
   12  CONTINUE
C
C SAVE PROCESSED SEDR VALUES BEFORE GETTING NEW NAV DATA
C
       IF (NAVTIME.LT.DATATIME) THEN
        DO J = 1,3
         DO I = 1,3
          OMTB(I,J) = MTB(I,J)
          OMTB5(I,J) = MTB5(I,J)
          OMHG(I,J) = MHG(I,J)
         END DO
        END DO
        OTN = TN
        OSPV(1) = SPV(1)
        OSPV(2) = SPV(2)
        OSPV(3) = SPV(3)
        OSPV(4) = SPV(4)
        OSPV(5) = SPV(5)
        OSPV(6) = SPV(6)
        ORANGE = RANGE
        OANG(1) = ANG(1)
        OANG(2) = ANG(2)
        GOTO 10  ! CHECK NEXT NAV RECORD
       END IF
C
C CALL SEDR PROCESSING ROUTINE 
C
       IF ( ICOORD.EQ.2 ) THEN
        CALL SEDRCRU(NTIME,NAV,PMAT,SPV,ANG,MTB,MHG,MTB5,RANGE)
       ELSE IF ( ICOORD.GE.3 ) THEN
        CALL SEDRENC(NTIME,NAV,PMAT,SPV,ANG,MTB,MHG,MTB5,RANGE)
       END IF
C
C INTERPOLATE ROTATION MATRICES, ANGLES, POSITION AND VELOCITY
C
       CALL SEDRIP(OTN,OMTB,TN,MTB,TD,MAT1,9,0) ! ROTATION MATRIX 1
       CALL SEDRIP(OTN,OMTB5,TN,MTB5,TD,MAT2,9,0) ! ROTATION MATRIX 2
       CALL SEDRIP(OTN,OMHG,TN,MHG,TD,MAT3,9,0) ! ROTATION MATRIX 3
       CALL SEDRIP(OTN,ORANGE,TN,RANGE,TD,ERANGE,1,0) ! RANGE
       CALL SEDRIP(OTN,OSPV,TN,SPV,TD,ESPV,6,0) ! S/C POSITION/VELOCITY
       CALL SEDRIP(OTN,OANG(1),TN,ANG(1),TD,EANG(1),1,3) ! LONGITUDE
       CALL SEDRIP(OTN,OANG(2),TN,ANG(2),TD,EANG(2),1,4) ! LATITUDE
       GOTO 900
  100 CONTINUE
      WRITE(IUNIT,800)
      WRITE(IUNIT,801) NTIME
      GOTO 900
  200 CONTINUE
      WRITE(IUNIT,802)
      WRITE(IUNIT,803) PTIME
  900 CONTINUE
C
      RETURN
  800 FORMAT(1X,'*GETSEDR*  WARNING!  ',
     &       'END OF NAVIGATION FILE REACHED.')
  801 FORMAT(1X,'*GETSEDR*  LAST NAVIGATION TIME TAG AT',
     &       1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3) 
  802 FORMAT(1X,'*GETSEDR*  WARNING!  '
     &       'END OF POINTING VECTOR FILE REACHED.')
  803 FORMAT(1X,'*GETSEDR*  LAST POINTING VECTOR TIME TAG AT',
     &       1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3) 
  804 FORMAT(1X,'*GETSEDR*  ERROR!  ',
     &       'INVALID POINTING VECTOR FILE HEADER.')
  805 FORMAT(1X,'*GETSEDR*  ERROR!  ',
     &       'INVALID NAVIGATION FILE HEADER.')
      END
      SUBROUTINE SEDRCRU(GMT,NAV,MPE,SPV,ANG,MTB,MHG,MTB5,RANGE)
C
C COMPUTE THE TRANFORMATION FROM THE EARTH MEAN ECLIPTIC AND EQUINOX
C OF 1950.0 (EME'50) COORDINATES TO HELIOGRAPHIC COORDINATES.
C
C INPUT NAVIGATION RECORD DATA (NAV) AND POINTING VECTOR DATA (MPE)
C
      INTEGER*2 GMT(6)
C
      REAL*4 NAV(252),RE(3),RH(3),VE(3),VH(3),SPV(6),ANG(2),      
     &       MHG(3,3),MTB(3,3),M5(3,3),MPE(3,3),MTB5(3,3)
C
      DATA TPIE / 6.2831853E0 /, ARAD / 1.495985E8 /
C
C M5 IS A PSEUDO-CONSTANT ORTHOGONAL MATRIX THAT MAPS VECTORS FROM THE EARTH-
C MEAN-ECLIPTIC, EQUINOX OF 1950 (EME'50) INERTIAL SYSTEM TO AN INERTIAL 
C HELIOGRAPHIC SYSTEM (IHG).
C
C MTB5 IS THE EME'50 TO HG TRANSFORMATION MATRIX (MTB5)
C
C MHG IS THE PAYLOAD TO HG TRANSFORMATION MATRIX (MHG)
C
      DATA M5 / 0.2576, 0.9662, 0.0   ,  
     &         -0.9585, 0.2555, 0.1262,  
     &          0.1219,-0.0325, 0.9920 /
C
      INCLUDE 'UNPACK.INC'
C 
C LOAD S/C POSITION VECTOR INTO ARRAY RE.
C
      RE(1) = NAV(13)
      RE(2) = NAV(14)
      RE(3) = NAV(15)
C
C LOAD S/C VELOCITY VECTOR INTO ARRAY VE.
C
      VE(1) = NAV(16)
      VE(2) = NAV(17)
      VE(3) = NAV(18)
C
      DO I = 1,3
       RH(I) = ( M5(1,I)*RE(1) + M5(2,I)*RE(2) + M5(3,I)*RE(3) ) / ARAD
       VH(I) =   M5(1,I)*VE(1) + M5(2,I)*VE(2) + M5(3,I)*VE(3)
      END DO
C
C GET IHG SPACECRAFT LATITUDE (THETA) AND LONGITUDE (BETA), AND 
C IHG TO HG TRANSFORMATION MATRIX (MTB).
C 
      X = RH(1)
      Y = RH(2)
      Z = RH(3)
      RANGE = SQRT(RH(1)**2+RH(2)**2+RH(3)**2)
C
C SPACECRAFT LATITUDINAL POSITION ANGLE
C
      THETA = ASIN(Z/RANGE)
C
C SPACECRAFT LONGITUDINAL POSITION ANGLE
C
      BETA = ATAN2(Y,X)
      IF (BETA.LT.0.0) BETA = BETA + TPIE
C
      DEN = SQRT(X**2 + Y**2)
      CSB = X/DEN
      SNB = Y/DEN
      SNT = Z/RANGE
      CST = SQRT(1.0 - SNT**2)
C
      MTB(1,1) =  CST*CSB
      MTB(1,2) = -SNB
      MTB(1,3) = -SNT*CSB
C
      MTB(2,1) =  CST*SNB
      MTB(2,2) =  CSB
      MTB(2,3) = -SNT*SNB
C
      MTB(3,1) =  SNT
      MTB(3,2) =  0.0
      MTB(3,3) =  CST
C
      SPV(1) = RH(1)
      SPV(2) = RH(2)
      SPV(3) = RH(3)
      SPV(4) = VH(1)
      SPV(5) = VH(2)
      SPV(6) = VH(3)
      ANG(1) = BETA
      ANG(2) = THETA
C
C GET EME'50 TO HG TRANSFORMATION MATRIX (MTB5)
C
      DO I = 1,3
       DO J = 1,3
        MTB5(I,J) = MTB(1,I)*M5(J,1) + 
     &              MTB(2,I)*M5(J,2) + 
     &              MTB(3,I)*M5(J,3)
       END DO
      END DO
C
C GET PAYLOAD TO HG TRANSFORMATION MATRIX (MHG)
C
      DO I = 1,3
       DO J = 1,3
        MHG(J,I) = MTB5(I,1)*MPE(1,J) + 
     &             MTB5(I,2)*MPE(2,J) + 
     &             MTB5(I,3)*MPE(3,J)
       END DO
      END DO
C
      RETURN
      END
C
C COMPUTE THE MATRICES FOR TRANSFORMATION FROM THE EARTH MEAN ECLIPTIC AND 
C EQUINOX OF 1950 (EME'50) COORDINATES TO PLANETARY SPHERICAL COORDINATES.
C
C MODIFIED TO SWAP VALUES OF LATITUDE AND LONGITUDE BEFORE RETURN 
C 07/31/95 SBK
C
      SUBROUTINE SEDRENC(GMT,NAV,PV,RJVE,ANG,S,ROT,Q,RANGE)
C
C NAV - Navigation record of either 126 or 252 word length.
C
C PV - Pointing vector data from words 12-20 of original MVS pointing
C      vector block.  
C         PV(1,1) = PITCH X COMPONENT
C         PV(2,1) = PITCH Y COMPONENT
C         PV(3,1) = PITCH Z COMPONENT
C         PV(1,2) = YAW X COMPONENT
C         PV(2,2) = YAW Y COMPONENT
C         PV(3,2) = YAW Z COMPONENT
C         PV(1,3) = ROLL X COMPONENT
C         PV(2,3) = ROLL Y COMPONENT
C         PV(3,3) = ROLL Z COMPONENT
C
C RJVE - Spacecraft range and velocity in planetary coordinates
C
      INTEGER*2 GMT(6)
C
      REAL*4 NAV(252),RJVE(6),ANG(2),PV(3,3),Q(3,3),ROT(3,3),S(3,3)
C
      REAL*8 DAYS,LONG,DRANGE,RA
C
C P-MATRICES TRANSFORM FROM CATERSIAN COMPONENTS REFERRED TO BY THE
C EARTH MEAN ECLIPTIC AND EQUINOX OF 1950 (EME'50) SYSTEM TO A 
C PLANETOCENTRIC CARTESIAN VERNAL EQUINOX SYSTEM (INERTIAL).
C
      REAL*4 PJ(3,3) / -7.214966E-1, +6.922823E-1, +1.370326E-2,
     &                 -6.922676E-1, -7.207863E-1, -3.511083E-2,
     &                 -1.442949E-2, -3.481867E-2, +9.992895E-1 /
C
      REAL*4 PS(3,3) / +9.915140E-1, -1.244780E-1, -3.749200E-2,
     &                 +9.252900E-2, +8.783090E-1, -4.690540E-1,
     &                 +9.131600E-2, +4.616040E-1, +8.823730E-1 /
C
      REAL*4 PU(3,3) / +9.737952E-1, -2.270260E-1, -1.349334E-2,
     &                 -4.308389E-2, -1.258954E-1, -9.911075E-1,
     &                 +2.233085E-1, +9.657171E-1, -1.323775E-1 /
C
      REAL*4 PN(3,3) / -7.093737E-1, -7.041542E-1, +3.091500E-2,
     &                 +6.092710E-1, -6.346587E-1, -4.753915E-1,
     &                 +3.543694E-1, -3.183947E-1, +8.792310E-1 / 
C
      REAL*4 ARADJ/71372.0E0/, ARADS/60330.0E0/, ARADU/25600.E0/,
     &       ARADN/24765.0E0/, DTR/0.17453293E0/
C
      INCLUDE 'UNPACK.INC'
C 
      IF ( ICOORD.EQ.4 ) THEN
C
C JUPITER
C
       DO I=1,3
        DO J=1,3
         Q(J,I) = PJ(1,I)*PV(1,J) + PJ(2,I)*PV(2,J) + PJ(3,I)*PV(3,J)
        END DO
        RJVE(I) = 
     &   (PJ(1,I)*NAV(19) + PJ(2,I)*NAV(20) + PJ(3,I)*NAV(21))/ARADJ
        RJVE(I+3) = 
     &   (PJ(1,I)*NAV(22) + PJ(2,I)*NAV(23) + PJ(3,I)*NAV(24))/ARADJ
       END DO
       ANG(1) = NAV(157)
       ANG(2) = NAV(159)
       RANGE = NAV(176)/ARADJ
C
      ELSE IF ( ICOORD.EQ.3 ) THEN
C
C SATURN
C
       DO I=1,3
        DO J=1,3
         Q(J,I) = PS(1,I)*PV(1,J) + PS(2,I)*PV(2,J) + PS(3,I)*PV(3,J)
        END DO
        RJVE(I) = 
     &   (PS(1,I)*NAV(19) + PS(2,I)*NAV(20) + PS(3,I)*NAV(21))/ARADS
        RJVE(I+3) = 
     &   (PS(1,I)*NAV(22) + PS(2,I)*NAV(23) + PS(3,I)*NAV(24))/ARADS
       END DO
       ANG(1) = NAV(76)
       ANG(2) = NAV(77) ! NEED TO SEE ORIGINAL CODE
       DAYS = GMT(2) - 1.0D0 + (GMT(3) + (GMT(4) + (GMT(5) +
     &        GMT(6)/1000.0D0)/60.0D0)/60.0D0)/24.0D0
       IF ( GMT(1).EQ.1981 ) DAYS = DAYS + 366.0D0
       DRANGE = DBLE(NAV(83))
       RA = NAV(103)
       LONG = 810.76D0*(DAYS - 3.860695545D-11*DRANGE) - RA
       ANG(2) = DMOD(LONG,360.D0)
       RANGE = NAV(83)/ARADS
C
      ELSE IF ( ICOORD.EQ.5 ) THEN
C
C URANUS
C
       DO I=1,3
        DO J=1,3
         Q(J,I) = PU(1,I)*PV(1,J) + PU(2,I)*PV(2,J) + PU(3,I)*PV(3,J)
        END DO
        RJVE(I) = 
     &   (PU(1,I)*NAV(19) + PU(2,I)*NAV(20) + PU(3,I)*NAV(21))/ARADU
        RJVE(I+3) = 
     &   (PU(1,I)*NAV(22) + PU(2,I)*NAV(23) + PU(3,I)*NAV(24))/ARADU
       END DO
       ANG(1) = NAV(139)
       ANG(2) = NAV(140)
       RANGE = NAV(153)/ARADU
C
      ELSE IF ( ICOORD.EQ.6 ) THEN
C
C NEPTUNE
C
       DO I=1,3
        DO J=1,3
         Q(J,I) = PN(1,I)*PV(1,J) + PN(2,I)*PV(2,J) + PN(3,I)*PV(3,J)
        END DO
        RJVE(I) = 
     &   (PN(1,I)*NAV(19) + PN(2,I)*NAV(20) + PN(3,I)*NAV(21))/ARADN
        RJVE(I+3) = 
     &   (PN(1,I)*NAV(22) + PN(2,I)*NAV(23) + PN(3,I)*NAV(24))/ARADN
       END DO
       ANG(1) = NAV(139)
       ANG(2) = NAV(140)
       RANGE = NAV(153)/ARADN
C
      END IF
C 
C COMPUTE PSI AND OMEGA
C
      SQ = SQRT(RJVE(1)**2+RJVE(2)**2+RJVE(3)**2)
      IF ( SQ.NE.0.0 ) THEN
       PSI = ACOS(RJVE(3)/SQ)
      ELSE
       PSI = 0.0
      END IF
      OMG = ATAN2(RJVE(2),RJVE(1))
C
C COMPUTES SINES AND COSINES OF ANGLES PSI AND OMEGA
C
      CPSI = COS(PSI)
      SPSI = SIN(PSI)
      COMG = COS(OMG)
      SOMG = SIN(OMG)
C
C COMPUTE S  MATRIX
C          1
C
      S(1,1) = SPSI*COMG
      S(1,2) = CPSI*COMG
      S(1,3) = -SOMG
      S(2,1) = SPSI*SOMG
      S(2,2) = CPSI*SOMG
      S(2,3) = COMG
      S(3,1) = CPSI
      S(3,2) = -SPSI
      S(3,3) = 0.0
C
C                             T
C COMPUTE ROTATION MATRIX S PA
C                          1
C
      DO I = 1,3
       DO J = 1,3
        ROT(J,I) = S(1,I)*Q(J,1) + S(2,I)*Q(J,2) + S(3,I)*Q(J,3)
       END DO
      END DO
C
C COMPUTE S  MATRIX
C          2
C
      COMG = COS(DTR*(360.0-ANG(2)))
      SOMG = SIN(DTR*(360.0-ANG(2)))
      S(1,1) = SPSI*COMG
      S(1,2) = CPSI*COMG
      S(1,3) = -SOMG
      S(2,1) = SPSI*SOMG
      S(2,2) = CPSI*SOMG
      S(2,3) = COMG
C
C          T    T
C COMPUTE S S PA
C          2 1
C
      DO I = 1,3
       DO J = 1,3
        Q(J,I) = S(I,1)*ROT(J,1) + S(I,2)*ROT(J,2) + S(I,3)*ROT(J,3)
       END DO
      END DO
C
C SWAP LATITUDE AND LONGITUDE FOR CONFORMANCE TO CRUISE MODE ROUTINES.
C
      TEMP = ANG(1)
      ANG(1) = ANG(2)
      ANG(2) = TEMP
C
      RETURN
      END
      SUBROUTINE SEDRIP (T1,B1,T2,B2,T,B,N,MD)
C
C****  =U2ADS.SEDRIP       05-OCT-79     ALLAN SILVER
C
      REAL*8 T1,T2,T
      REAL*4 B1(1),B2(1),B(1)
C --> REAL*4 RL(3)/0.,-3.1415927E0,0./       <-- REPLACED 86-APR-19,
C --> REAL*4 RU(3)/+3.1415927E0,+3.1415927E0,+360./    JAJONES, SAR
      REAL*4 RL(4)/0.,-3.1415927E0,0.,-90./
      REAL*4 RU(4)/+3.1415927E0,+3.1415927E0,+360.,+90./
C***********************************************************************
C******** LINEAR INTERPOLATION WITHIN LIMITS & NORMALIZATION ***********
C***********************************************************************
C****  INDEPENDENT VARIABLE     FUNCTION VALUE(S)                      *
C****    T1 (INPUT)              B1(1...N) (INPUT)                     *
C****    T2 (INPUT)              B2(1...N) (INPUT)                     *
C****     T (INPUT)               B(1...N) (OUTPUT; LINEAR INTERP.)    *
C***********************************************************************
C****         MD         LIMITS OF FUNCTION                            *
C****          0          UNLIMITED                                    *
C****          1          0 TO PI                                      *
C****          2          - PI TO + PI                                 *
C****          3          0 TO 360                                     *
C****          4          - 90 TO + 90  ( <-- ADDED BY JAJ 86-APR19 )  *
C***********************************************************************
C****  T1<0  ---> DO NOT USE B1  ****  T2<0  ---> DO NOT USE B2        *
C****                    T1=T2  ---> USE B1 ONLY                       *
C***********************************************************************
C**** FOR N=9 & MD=0, TREAT B AS A 3X3 MATRIX & NORMALIZE BY ROWS      *
C***********************************************************************
C
C**** COMPUTE RANGE AND HALF RANGE
C
      IF (MD.NE.0)  THEN
       RG = RU(MD) - RL(MD)
       RGH = .5E0*RG
      END IF
C
C**** COMPUTE LINEAR INTERPOLATION COEFFICIENTS
C
      IF (T1.GE.0.D0) THEN
       IF ( T2.GE.0.D0 .AND. T1.NE.T2 ) THEN 
        ALF = (T2 - T)/(T2 - T1)
        BET = (T - T1)/(T2 - T1)
       ELSE
        ALF = 1.E0
        BET = 0.E0
       END IF
      ELSE IF (T2.LT.0.D0)  THEN
       WRITE (6,2001) T1,T2  ! ERROR
       RETURN
      ELSE
       ALF = 0.E0
       BET = 1.E0
      END IF
C
C**** DO THE INTERPOLATION
C
      DO I = 1,N
       B(I) = ALF*B1(I) + BET*B2(I)
       IF ( MD.NE.0 ) THEN
C**** PUT B1 AND B2 WITHIN HALF-RANGE OF EACH OTHER
        K = INT((B2(I) - B1(I))/RGH)
        B(I) = B(I) + RGH*ALF*(K + MOD(K,2))
C**** PUT B WITHIN LIMITS
        B(I) = B(I) - RG*INT((B(I) - RL(MD))/RG)
        IF ( B(I).LT.RL(MD) ) B(I) = B(I) + RG
       END IF
      END DO
      IF (N.NE.9.OR.MD.NE.0) RETURN  ! EXIT ROUTINE IF NOT 3X3 MATRIX
C
C**** NORMALIZE THE MATRIX
C
      DO J = 1,7,3
       FCT = SQRT(B(J)**2 + B(J+1)**2 + B(J+2)**2)
       IF ( FCT.NE.0.E0 ) THEN
        B(J) = B(J)/FCT
        B(J+1) = B(J+1)/FCT
        B(J+2) = B(J+2)/FCT
       END IF
      END DO
C
      RETURN
 2001 FORMAT (' SEDRIP ERROR: ',1P2D15.7)
      END
      SUBROUTINE EXTRACTDATE(TIME,PVRD)                  
C
C EXTRACT POINTING VECTOR RECORD TIME
C
      INTEGER*2 TIME(6)
      CHARACTER PVRD*77
C                                                 
      READ(UNIT=PVRD,FMT=89) TIME
C
      RETURN                                                            
  89  FORMAT(I2,1X,I3,3(1X,I2),1X,I3)
      END
      SUBROUTINE ELPSTIME(START,TIME,DAYS)
C
C THIS ROUTINE COMPUTES A POSITIVE ELAPSED TIME BETWEEN TWO INTEGER
C ARRAYED TIME VALUES WHERE:
C
C   TIME ELEMENT 1 = YEAR
C   TIME ELEMENT 2 = DAY
C   TIME ELEMENT 3 = HOUR
C   TIME ELEMENT 4 = MINUTE
C   TIME ELEMENT 5 = SECOND
C   TIME ELEMENT 6 = MILLISECOND
C
      INTEGER*2 START(6),TIME(6)
      REAL*8 DAYS,REALTIME
C
C TEST FOR PROPER INPUT.  2 DIGIT YEAR EXPECTED (LESS THAN 77 ASSUMED 21ST
C CENTURY)
C
      IF ((MOD(START(1),4).EQ.0.AND.
     &     START(2).GT.366).OR.
     &    (MOD(START(1),4).NE.0.AND.
     &     START(2).GT.365).OR.
     &    START(1).GT.99.OR.
     &    START(3).GT.23.OR.
     &    START(4).GT.59.OR. 
     &    START(5).GT.59.OR.
     &    START(6).GT.999.OR.
     &    START(1).LT.0.OR.
     &    START(2).LT.1.OR.
     &    START(3).LT.0.OR.
     &    START(4).LT.0.OR.
     &    START(5).LT.0.OR.
     &    START(6).LT.0) THEN
       WRITE(6,'(1X,''*ELPSTIME*  INVALID START TIME: '',
     &       I2,''-'',I3.3,''/'',I2.2,'':'',I2.2,'':'',
     &       I2.2,''.'',I3.3)') START
       DAYS = -1.0D0
       RETURN
      END IF
C
      IF ((MOD(TIME(1),4).EQ.0.AND.
     &     TIME(2).GT.366).OR.
     &    (MOD(TIME(1),4).NE.0.AND.
     &     TIME(2).GT.365).OR.
     &    TIME(1).GT.99.OR.
     &    TIME(3).GT.23.OR.
     &    TIME(4).GT.59.OR. 
     &    TIME(5).GT.59.OR.
     &    TIME(6).GT.999.OR.
     &    TIME(1).LT.0.OR.
     &    TIME(2).LT.1.OR.
     &    TIME(3).LT.0.OR.
     &    TIME(4).LT.0.OR.
     &    TIME(5).LT.0.OR.
     &    TIME(6).LT.0) THEN
       WRITE(6,'(1X,''*ELPSTIME*  INVALID CURRENT TIME: '',
     &       I2,''-'',I3.3,''/'',I2.2,'':'',I2.2,'':'',
     &       I2.2,''.'',I3.3)') TIME
       DAYS = -1.0D0
       RETURN
      END IF
C
      DAYS = DBLE(START(2)) + 
     &       DBLE(START(3))/24.0D0 +
     &       DBLE(START(4))/24.0D0/60.0D0 +
     &       DBLE(START(5))/24.0D0/60.0D0/60.0D0 +
     &       DBLE(START(6))/24.0D0/60.0D0/60.0D0/1000.0D0
      DAYS = - DAYS
C
C CONVERT YEARS DIFFERENCE INTO DECIMAL DAYS
C
      IF (TIME(1).LT.77) THEN
       ICENT = 100
      ELSE
       ICENT = 0
      END IF
C
      DO I = START(1), TIME(1) - 1 + ICENT
C
       IF (MOD(I,4).EQ.0) THEN
        DAYS = DAYS + 366.0D0
       ELSE 
        DAYS = DAYS + 365.0D0
       END IF
C
      END DO
C
C ADD DECIMAL DAYS OF CURRENT YEAR TO SUM
C
      DAYS = DAYS +
     &       DBLE(TIME(2)) + 
     &       DBLE(TIME(3))/24.0D0 +
     &       DBLE(TIME(4))/24.0D0/60.0D0 +
     &       DBLE(TIME(5))/24.0D0/60.0D0/60.0D0 +
     &       DBLE(TIME(6))/24.0D0/60.0D0/60.0D0/1000.0D0
C
      RETURN
      END
      REAL*8 FUNCTION REALTIME(TIME)
C
C CONVERT INTEGER CALENDAR TIME INTO DECIMAL YEAR REAL TIME.
C
      INTEGER*2 TIME(6)
      REAL*8 DAYS
C
      DAYS = 365.0D0
      IF (MOD(TIME(1),4).EQ.0) DAYS = 366.0D0
      REALTIME = DBLE(TIME(1)) + 
     &           DBLE(TIME(2)-1)/DAYS +
     &           DBLE(TIME(3))/24.0D0/DAYS +
     &           DBLE(TIME(4))/60.0D0/24.0D0/DAYS +
     &           DBLE(TIME(5))/60.0D0/60.0D0/24.0D0/DAYS +
     &           DBLE(TIME(6))/1000.0D0/60.0D0/60.0D0/24.0D0/DAYS
C
C ASSUME 2 DIGIT YEAR.  ANY YEAR BEFORE VOYAGER LAUNCH (77) IS INTO NEXT
C CENTURY.
C
      IF (TIME(1).LT.77) REALTIME = REALTIME + 100.0D0
C
      RETURN                                                            
      END
