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
      CHARACTER IHDR(45)*4,DSN*50,TYPE*4
      INTEGER*2 NTIME(6),INAVTIME(6)
      INTEGER*4 INAV(126,2)
      REAL*4 RNAV(126,2),NAV(6),PV(3,3),RH(3),VH(3),
     &       MTB(3,3),MHG(3,3),MTB5(3,3),ANG(2),SPV(6)
      REAL*8 REALTIME,NAVTIME
C
      DATA ICALL/0/,ARAD/1.495985E8/,PV/9*0.0/
C
      EQUIVALENCE ( RNAV(1,1), INAV(1,1) ), ( IHDR(1), INAV(1,1) )
C
      WRITE(6,*)
      WRITE(6,*) 'ENTER INPUT DSN'
      READ(5,'(A)') DSN
      OPEN(41,FILE=DSN,STATUS='OLD',FORM='UNFORMATTED',
     &     RECORDTYPE='VARIABLE',RECL=8191,READONLY)
C
      WRITE(6,*)
      WRITE(6,*) 'ENTER RECORD STEP COUNT >= 1'
      READ(5,*) ISTEP
      IF ( ISTEP.LE.0 ) ISTEP = 1
C
C GET NAVIGATION FILE HEADER RECORD.  CONTAINS CHARACTER AND INTEGER DATA.
C
      READ(41) IHDR
      WRITE(6,*)
      WRITE(6,800) IHDR(1),IHDR(2),IHDR(4),IHDR(5),INAV(3,1),IHDR(14)
      WRITE(6,*)
      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.
C
      IWRITE = 0
      NCNT = 0
   10 CONTINUE
       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
       NCNT = NCNT + 1
C
C GET NAVIGATION RECORD TIME TAG.  USE 2 DIGIT YEAR.
C
       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 EME'50 SUN-CENTERED X,Y,Z COORDINATE S/C POSITION
C
        NAV(1) = RNAV(13,1)/ARAD
        NAV(2) = RNAV(14,1)/ARAD
        NAV(3) = RNAV(15,1)/ARAD
C
C EME'50 SUN-CENTERED X,Y,Z COMPONENT S/C VELOCITY
C
        NAV(4) = RNAV(16,1)
        NAV(5) = RNAV(17,1)
        NAV(6) = RNAV(18,1)
C
C RANGE  SUN - S/C
C
        IF ( TYPE.EQ.'LAUN' .OR.
     &       TYPE.EQ.'CRUI' .OR.
     &       TYPE.EQ.'XCRU' ) THEN
         RANGE = RNAV(51,1)/ARAD
        ELSE IF ( TYPE.EQ.'JUPI' ) THEN
         RANGE = RNAV(47,2)/ARAD
        ELSE IF ( TYPE.EQ.'SATU' ) THEN
         RANGE = RNAV(82,1)/ARAD
        ELSE IF ( TYPE.EQ.'URAN' ) THEN
         RANGE = RNAV(24,2)/ARAD
        ELSE IF ( TYPE.EQ.'NEPT' ) THEN
         RANGE = RNAV(24,2)/ARAD
        ELSE
         WRITE(6,*)
         WRITE(6,*) 'UNRECOGNIZED NAV DATA BLOCK ID:  ',TYPE
         RANGE = SQRT(NAV(1)**2+NAV(2)**2+NAV(3)**2)
        END IF
C
C       CALL SEDRPROC(NAV,PV,RH,VH,THETA,BETA,MTB,MHG,MTB5)
        CALL SEDRCRU(NTIME,RNAV,PV,SPV,ANG,MTB,MHG,MTB5,RANGE)
        THETA = ANG(2)
        BETA = ANG(1)
        RH(1) = SPV(1)
        RH(2) = SPV(2)
        RH(3) = SPV(3)
C
       IF (MOD(NCNT-1,ISTEP).EQ.0) THEN
        IF ( IWRITE.EQ.0 ) THEN
         WRITE(6,*)
         WRITE(6,*)
         WRITE(6,808)
         WRITE(6,809)
         WRITE(6,*)
        END IF
        BETA = BETA * 180.0/3.14159
        THETA = THETA * 180.0/3.14159
        WRITE(6,810) NTIME,RANGE,BETA,THETA,(RH(I),I=1,3)
c        WRITE(6,'(2X,3(3(1X,E10.4),3X))') (MTB(I,1),I=1,3),
c     &   (MTB5(1,I),I=1,3),(MHG(I,1),I=1,3)
c        WRITE(6,'(2X,3(3(1X,E10.4),3X))') (MTB(I,2),I=1,3),
c     &   (MTB5(2,I),I=1,3),(MHG(I,2),I=1,3)
c        WRITE(6,'(2X,3(3(1X,E10.4),3X)/)') (MTB(I,3),I=1,3),
c     &   (MTB5(3,I),I=1,3),(MHG(I,3),I=1,3)
        IWRITE = IWRITE + 1
       END IF
C
      GOTO 10
  100 CONTINUE
      WRITE(6,*)
      WRITE(6,*) '***WARNING!  END OF NAVIGATION FILE REACHED***'
  900 CONTINUE
C
  800 FORMAT(1X,A4,3X,A4,3X,2A4,3X,'FLIGHT ',I2,3X,A4)
  808 FORMAT(8X,'TIME',11X,'RANGE',3X,'BETA',4X,'THETA',
     &       6X,'X',9X,'Y',9X,'Z')
  809 FORMAT(24X,'(AU)',3X,'(deg)',3X,'(deg)',5X,'(AU)',6X,'(AU)',
     &       6X,'(AU)')
  810 FORMAT(1X,I2.2,'-',I3.3,'/',I2.2,2(':',I2.2),'.',I3.3,1X,
     &       F8.4,2(1X,F7.3),3(1X,F9.4))
      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 IN TO NEXT
C CENTURY.
C
      IF (TIME(1).LT.77) REALTIME = REALTIME + 100.0D0
C
      RETURN
      END
      SUBROUTINE SEDRPROC(NVBLK,MPE,RH,VH,THETA,BETA,MTB,MHG,MTB5)
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 (NVBLK) AND POINTING VECTOR DATA (MPE)
C
C MATRICES ARE REPRESENTED IN ROW MAJOR ORDER, IE;
C    11 12 13
C    21 22 23
C    31 32 33
C
      REAL*4 NVBLK(6),RE(3),RH(3),VE(3),VH(3)
      REAL*4 MHG(3,3),MTB(3,3),MSUN(3,3),M5(3,3),MPE(3,3),MTB5(3,3)
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
      DATA M5 / 0.2576, -0.9585,  0.1219,
     &          0.9662,  0.2555, -0.0325,
     &          0.0   ,  0.1262,  0.9920 /
C 
C LOAD S/C POSITION VECTOR INTO ARRAY RE.
C
      RE(1) = NVBLK(1)
      RE(2) = NVBLK(2)
      RE(3) = NVBLK(3)
C
C LOAD S/C VELOCITY VECTOR INTO ARRAY VE.
C
      VE(1) = NVBLK(4)
      VE(2) = NVBLK(5)
      VE(3) = NVBLK(6)
C
C GET IHG SPACECRAFT LATITUDE (THETA) AND LONGITUDE (BETA), AND 
C IHG TO HG TRANSFORMATION MATRIX (MTB).
C 
      CALL POS(RE,VE,M5,RH,VH,THETA,BETA,MTB)
C
C GET PAYLOAD TO IHG TRANSFORMATION MATRIX (MSUN)
C
      CALL IHG(MSUN,M5,MPE)
C
C GET PAYLOAD TO HG TRANSFORMATION MATRIX (MHG)
C
      CALL HG(MHG,MTB,MSUN)
C
C GET EME'50 TO HG TRANSFORMATION MATRIX (MTB5)
C
      CALL MPRD33(MTB5,MTB,M5)
C
      RETURN
      END
      SUBROUTINE IHG(MSUN,M5,MPE)
C
C COMPUTE THE PAYLOAD TO IHG TRANSFORMATION MATRIX (MSUN).
C
      REAL*4 MSUN(3,3),M5(3,3),MPE(3,3)
C
      MSUN(1,1) = M5(1,1)*MPE(1,1) + M5(1,2)*MPE(2,1) + M5(1,3)*MPE(3,1)
      MSUN(1,2) = M5(1,1)*MPE(1,2) + M5(1,2)*MPE(2,2) + M5(1,3)*MPE(3,2)
      MSUN(1,3) = M5(1,1)*MPE(1,3) + M5(1,2)*MPE(2,3) + M5(1,3)*MPE(3,3)
C
      MSUN(2,1) = M5(2,1)*MPE(1,1) + M5(2,2)*MPE(2,1) + M5(2,3)*MPE(3,1)
      MSUN(2,2) = M5(2,1)*MPE(1,2) + M5(2,2)*MPE(2,2) + M5(2,3)*MPE(3,2)
      MSUN(2,3) = M5(2,1)*MPE(1,3) + M5(2,2)*MPE(2,3) + M5(2,3)*MPE(3,3)
C
      MSUN(3,1) = M5(3,1)*MPE(1,1) + M5(3,2)*MPE(2,1) + M5(3,3)*MPE(3,1)
      MSUN(3,2) = M5(3,1)*MPE(1,2) + M5(3,2)*MPE(2,2) + M5(3,3)*MPE(3,2)
      MSUN(3,3) = M5(3,1)*MPE(1,3) + M5(3,2)*MPE(2,3) + M5(3,3)*MPE(3,3)
C
      RETURN
      END 
      SUBROUTINE POS(RE,VE,M5,RH,VH,THETA,BETA,MTB)
C
C ROTATE SPACECRAFT POSITION AND VELOCITY FROM EME'50 COORDINATES 
C TO IHG COORDINATES.
C COMPUTE SPACECRAFT LATITUDE AND LONGITUDE IN IHG COORDINATES.
C COMPUTE THE IHG TO HG ROTATION MATRIX (MTB). 
C 
      REAL*4 RE(3),VE(3),M5(3,3),RH(3),VH(3),MTB(3,3)
C
      DATA TPIE/6.2831853E0/
C
C ROTATE S/C POSITION VECTOR INTO IHG COORDINATES.
C
      RH(1) = M5(1,1)*RE(1) + M5(1,2)*RE(2) + M5(1,3)*RE(3)
      RH(2) = M5(2,1)*RE(1) + M5(2,2)*RE(2) + M5(2,3)*RE(3)
      RH(3) = M5(3,1)*RE(1) + M5(3,2)*RE(2) + M5(3,3)*RE(3)
C
C ROTATE S/C VELOCITY VECTOR INTO IHG COORDINATES.
C
      VH(1) = M5(1,1)*VE(1) + M5(1,2)*VE(2) + M5(1,3)*VE(3)
      VH(2) = M5(2,1)*VE(1) + M5(2,2)*VE(2) + M5(2,3)*VE(3)
      VH(3) = M5(3,1)*VE(1) + M5(3,2)*VE(2) + M5(3,3)*VE(3)
C
      X = RH(1)
      Y = RH(2)
      Z = RH(3)
C
C SPACECRAFT RANGE
C
      R = SQRT(X**2 + Y**2 + Z**2)
C
C SPACECRAFT LATITUDINAL POSITION ANGLE
C
      THETA = ASIN(Z/R)
C
C SPACECRAFT LONGITUDINAL POSITION ANGLE
C
      BETA = ATAN2(Y,X)
      IF (BETA.LT.0.0) BETA = BETA + TPIE
C
      CSB = X/SQRT(X**2 + Y**2)
      SNB = Y/SQRT(X**2 + Y**2)
      SNT = Z/R
      CST = SQRT(1.0 - SNT**2)
C
      MTB(1,1) =  CST*CSB
      MTB(1,2) =  CST*SNB
      MTB(1,3) =  SNT
C
      MTB(2,1) = -SNB
      MTB(2,2) =  CSB
      MTB(2,3) =  0.0
C
      MTB(3,1) = -SNT*CSB
      MTB(3,2) = -SNT*SNB
      MTB(3,3) =  CST
C
      RETURN
      END
      SUBROUTINE HG(MHG,MTB,MSUN)
C
C COMPUTE PAYLOAD TO HG ROTATION MATRIX (MHG).
C
      REAL*4 MHG(3,3),MTB(3,3),MSUN(3,3)
C
      MHG(1,1) = MTB(1,1)*MSUN(1,1) + MTB(1,2)*MSUN(2,1) + 
     &           MTB(1,3)*MSUN(3,1)
      MHG(1,2) = MTB(1,1)*MSUN(1,2) + MTB(1,2)*MSUN(2,2) + 
     &           MTB(1,3)*MSUN(3,2)
      MHG(1,3) = MTB(1,1)*MSUN(1,3) + MTB(1,2)*MSUN(2,3) + 
     &           MTB(1,3)*MSUN(3,3)
C
      MHG(2,1) = MTB(2,1)*MSUN(1,1) + MTB(2,2)*MSUN(2,1) + 
     &           MTB(2,3)*MSUN(3,1)
      MHG(2,2) = MTB(2,1)*MSUN(1,2) + MTB(2,2)*MSUN(2,2) + 
     &           MTB(2,3)*MSUN(3,2)
      MHG(2,3) = MTB(2,1)*MSUN(1,3) + MTB(2,2)*MSUN(2,3) + 
     &           MTB(2,3)*MSUN(3,3)
C
      MHG(3,1) = MTB(3,1)*MSUN(1,1) + MTB(3,2)*MSUN(2,1) + 
     &           MTB(3,3)*MSUN(3,1)
      MHG(3,2) = MTB(3,1)*MSUN(1,2) + MTB(3,2)*MSUN(2,2) + 
     &           MTB(3,3)*MSUN(3,2)
      MHG(3,3) = MTB(3,1)*MSUN(1,3) + MTB(3,2)*MSUN(2,3) + 
     &           MTB(3,3)*MSUN(3,3)
C
      RETURN
      END
      SUBROUTINE MPRD31(R,A,B)
C
C PERFORM MATRIX MULTIPLICATION (3X3 BY 3X1)
C
      REAL*4 R(3,1),A(3,3),B(3,1)
C
      R(1,1) = A(1,1)*B(1,1) + A(1,2)*B(2,1) + A(1,3)*B(3,1)
C
      R(2,1) = A(2,1)*B(1,1) + A(2,2)*B(2,1) + A(2,3)*B(3,1)
C
      R(3,1) = A(3,1)*B(1,1) + A(3,2)*B(2,1) + A(3,3)*B(3,1)
C
      RETURN
      END
      SUBROUTINE MPRD33(R,A,B)
C
C PERFORM MATRIX MULTIPLICATION (3X3 BY 3X3)
C
      REAL*4 R(3,3),A(3,3),B(3,3)
C
      R(1,1) = A(1,1)*B(1,1) + A(1,2)*B(2,1) + A(1,3)*B(3,1)
      R(1,2) = A(1,1)*B(1,2) + A(1,2)*B(2,2) + A(1,3)*B(3,2)
      R(1,3) = A(1,1)*B(1,3) + A(1,2)*B(2,3) + A(1,3)*B(3,3)
C
      R(2,1) = A(2,1)*B(1,1) + A(2,2)*B(2,1) + A(2,3)*B(3,1)
      R(2,2) = A(2,1)*B(1,2) + A(2,2)*B(2,2) + A(2,3)*B(3,2)
      R(2,3) = A(2,1)*B(1,3) + A(2,2)*B(2,3) + A(2,3)*B(3,3)
C
      R(3,1) = A(3,1)*B(1,1) + A(3,2)*B(2,1) + A(3,3)*B(3,1)
      R(3,2) = A(3,1)*B(1,2) + A(3,2)*B(2,2) + A(3,3)*B(3,2)
      R(3,3) = A(3,1)*B(1,3) + A(3,2)*B(2,3) + A(3,3)*B(3,3)
C
      RETURN
      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
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