C$Procedure      STDIFF ( State differences )
 
      SUBROUTINE STDIFF ( STA, STB, NITR, TIMES, UNIT, WORST, NONZER )
 
C$ Abstract
C
C     Computes the differences between the state vectors in STA
C     and those in STB and writes the summary to the file
C     connected to UNIT.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     EPHEMERIS
C
C$ Declarations
 
      DOUBLE PRECISION      STA   ( 6, * )
      DOUBLE PRECISION      STB   ( 6, * )
      INTEGER               NITR
      DOUBLE PRECISION      TIMES (    * )
      INTEGER               UNIT
      DOUBLE PRECISION      WORST ( 6 )
      LOGICAL               NONZER
 
      INTEGER               MAXN
      PARAMETER           ( MAXN = 1000 )
 
C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     STA        I   Array of state vectors.
C
C     STB        I   Another array of state vectors that represent
C                    the same states as those in STA but generated
C                    in a different way.
C
C     NITR       I   The number of states in STA or STB.
C
C     TIMES      I   Epochs of the states.
C
C     UNIT       I   Logical unit connected to a text file.
C
C     WORST      O   Array containing statistics for the worst case.
C
C     NONZER     O   Flag indicating if all calculations were
C                    performed.
C
C     MAXN       P   Maximum value for N.
C
C$ Detailed_Input
C
C     STA,
C     STB         are arrays of state vectors.  A state vector
C                 is a six element array that describes the position
C                 and velocity of one body with respect to another
C                 at a particular time.  The first three
C                 elements are the X, Y, and Z coordinates
C                 of the position, and the last three elements
C                 are the coordinates of the velocity.
C
C                 The states in STA and STB should be for the same
C                 body, with respect to the same center, in the same
C                 reference frame, and associated with the same epochs.
C                 They differ only in that they come from different
C                 sources.
C
C     NITR        The number of states in STA or STB.
C
C     TIMES       is an array of epochs for which the states in
C                 STA and STB are valid.
C
C     UNIT        is the logical unit connected to a text file.
C                 The states in STA are compared to those in STB
C                 and a summary of the tests is written to the file.
C
C$ Detailed_Output
C
C     See UNIT above.
C
C     WORST       is an array which returns some of the statistics
C                 generated from the pair of states with the largest
C                 relative difference in position.
C
C                 WORST ( 1 ) = relative difference in position
C                 WORST ( 2 ) = magnitute of position difference
C                 WORST ( 3 ) = down track difference in position
C                 WORST ( 4 ) = out of plane difference in position
C                 WORST ( 5 ) = in plane difference in position
C                 WORST ( 6 ) = down track difference divided by speed
C
C     NONZER      flag indicating if view frame components could be
C                 computed.
C
C$ Parameters
C
C     MAXN        the maximum number of states is this value plus 1.
C
C$ Exceptions
C
C     None.
C
C$ Files
C
C     See UNIT above.
C
C$ Particulars
C
C     STA and STB are arrays of NITR states of some body with
C     respect to some center at epochs given by the array TIMES.
C     These vectors should all be in the same inertial reference
C     frame.  This routine compares the states given by the two arrays
C     and then writes to a text file the statistics generated by the
C     comparison tests.
C
C     The following comparison tests are performed on the states from
C     the two arrays:
C
C     1) The average and maximum relative differences in the position
C        and velocity vectors.
C
C     2) The average and maximum magnitudes of the position and
C        velocity difference vectors.
C
C     3) The components of the position difference vectors are
C        calculated in "view frame" coordinates.  The output listing
C        includes:
C
C           a) The components in view frame coordinates.
C
C           b) The absolute value of these components.
C
C           c) The components for the difference vector associated
C              with the states that had the largest relative difference
C              in position.
C
C$ Examples
C
C     See the programs TEST_NIOSPK and TEST_SPKSPK
C
C$ Restrictions
C
C     The states in STA and STB should be for the same body, with
C     respect to the same center, in the same reference frame, and
C     associated with the same epochs.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     J.M. Lynch     (JPL)
C     J.E. McLean    (JPL)
C
C$ Version
C
C-    Beta Version 3.0.3, 08-OCTOBER-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitely given.  New
C        environments are WIN-NT
C
C-    Beta Version 3.0.2, 28-JULY-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitely given.  New
C        environments are PC-DIGITAL, SGI-O32 and SGI-N32.
C
C-    Beta Version 3.0.1, 18-MARCH-1999 (WLT)
C
C        The environment lines were expanded so that the supported
C        environments are now explicitely given.  Previously,
C        environments such as SUN-SUNOS and SUN-SOLARIS were implied
C        by the environment label SUN.
C
C-    Beta Version 3.0.0, 01-MAY-1992 (JML)
C
C        Added tests in view frame coordinates.
C
C        Deleted tests on individual components of states.
C
C        Added output arguments which return the worst case results.
C
C        Changed the indexing of the arrays from 0:NITR to 1:NITR.
C
C-    Beta Version 2.0.0, 21-MAY-1991 (JEM)
C
C        Added computation of average absolute deviation.
C
C-    Beta Version 1.0.0, 01-MAY-1991 (JEM)
C
C-&
 
C
C     SPICELIB functions
C
      LOGICAL               RETURN
 
      DOUBLE PRECISION      VNORM
      DOUBLE PRECISION      VREL
      DOUBLE PRECISION      VDOT
 
C
C     Local variables
C
      INTEGER               J
      INTEGER               I
      INTEGER               NLBMAX
      INTEGER               NBMAX
 
      CHARACTER*(20)        OUTCH    ( 4 )
      CHARACTER*(70)        POSSTR
      CHARACTER*(70)        VELSTR
 
      DOUBLE PRECISION      DIFF     ( 6 )
      DOUBLE PRECISION      POSDIF   ( 3 )
      DOUBLE PRECISION      VELDIF   ( 3 )
      DOUBLE PRECISION      POSCMP   ( 3 )
      DOUBLE PRECISION      VELCMP   ( 3 )
      DOUBLE PRECISION      SUMP     ( 3 )
      DOUBLE PRECISION      SUMV     ( 3 )
      DOUBLE PRECISION      SUMPA    ( 3 )
      DOUBLE PRECISION      SUMVA    ( 3 )
      DOUBLE PRECISION      MAXPOS   ( 3 )
      DOUBLE PRECISION      MAXVEL   ( 3 )
      DOUBLE PRECISION      AVPO     ( 3 )
      DOUBLE PRECISION      AVVE     ( 3 )
      DOUBLE PRECISION      AVPOA    ( 3 )
      DOUBLE PRECISION      AVVEA    ( 3 )
      DOUBLE PRECISION      AXIS     ( 3, 3 )
      DOUBLE PRECISION      DELTA
      DOUBLE PRECISION      POSMAG
      DOUBLE PRECISION      VELMAG
      DOUBLE PRECISION      POSSUM
      DOUBLE PRECISION      VELSUM
      DOUBLE PRECISION      SUMDT
      DOUBLE PRECISION      SUMADT
      DOUBLE PRECISION      AVPOS
      DOUBLE PRECISION      AVVEL
      DOUBLE PRECISION      AVADT
      DOUBLE PRECISION      AVDT
      DOUBLE PRECISION      POSREL
      DOUBLE PRECISION      VELREL
      DOUBLE PRECISION      RELPSM
      DOUBLE PRECISION      RELVSM
      DOUBLE PRECISION      MAXPMG
      DOUBLE PRECISION      MAXVMG
      DOUBLE PRECISION      MAXPRL
      DOUBLE PRECISION      MAXVRL
      DOUBLE PRECISION      MAXADT
      DOUBLE PRECISION      AVRELP
      DOUBLE PRECISION      AVRELV
 
      LOGICAL               MAXTRU
 
C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'STDIFF' )
      END IF
 
C
C     Initial values.
C
      DO I = 1, 3
         SUMP   (I)   =   0.D0
         SUMV   (I)   =   0.D0
         SUMPA  (I)   =   0.D0
         SUMVA  (I)   =   0.D0
         MAXPOS (I)   =   0.D0
         MAXVEL (I)   =   0.D0
      END DO
 
      DO I = 1, 6
         WORST (I) = 0.D0
      END DO
 
      SUMDT   =   0.D0
      SUMADT  =   0.D0
      POSSUM  =   0.D0
      VELSUM  =   0.D0
      RELPSM  =   0.D0
      RELVSM  =   0.D0
      MAXPMG  =   0.D0
      MAXVMG  =   0.D0
      MAXPRL  =   0.D0
      MAXVRL  =   0.D0
      MAXADT  =   0.D0
 
      NLBMAX  =   0
      NBMAX   =   0
 
      NONZER  =  .TRUE.
 
C
C     After subracting the state vectors for the Jth
C     time, first compare the individual components of the difference
C     vector and then perform the tests on the vector itself.
C
      DO J = 1, NITR
 
C
C        Find the difference vector between the two states (STA-STB).
C
         CALL VSUBG ( STA(1,J), STB(1,J), 6, DIFF )
 
C
C        Now perform tests on the actual difference vectors.
C
         CALL VEQU ( DIFF(1), POSDIF )
         CALL VEQU ( DIFF(4), VELDIF )
 
C
C        Find the magnitudes and relative error of the position and
C        velocity differecne vectors.
C
         POSMAG = VNORM ( POSDIF )
         VELMAG = VNORM ( VELDIF )
 
         POSREL = VREL ( STA(1,J), STB(1,J) )
 
         VELREL = VREL ( STA(4,J), STB(4,J) )
 
         POSSUM = POSSUM + POSMAG
         VELSUM = VELSUM + VELMAG
 
         RELPSM = RELPSM + POSREL
         RELVSM = RELVSM + VELREL
 
         MAXTRU = .FALSE.
 
         IF ( POSMAG .GT. MAXPMG ) THEN
            MAXPMG    = POSMAG
         END IF
 
         IF ( VELMAG .GT. MAXVMG ) THEN
            MAXVMG = VELMAG
         END IF
 
         IF ( POSREL .GT. MAXPRL ) THEN
 
C
C           We are going to return information on the case with the
C           largest relative difference in position.
C
            MAXTRU    = .TRUE.
            MAXPRL    = POSREL
            WORST(1)  = POSREL
            WORST(2)  = POSMAG
 
         END IF
 
         IF ( VELREL .GT. MAXVRL ) THEN
            MAXVRL = VELREL
         END IF
 
C
C        Compute the components of the position and velocity difference
C        vectors in view frame coordinates.
C
C           AXIS(1): is parallel to the direction of motion given by
C                    STA
C
C           AXIS(2): is normal to the orbit plane defined by STA
C
C           AXIS(3): is AXIS(1) X AXIS(2)
C
C
         CALL VEQU ( STA(4,J), AXIS(1,1) )
 
         CALL VCRSS ( STA (1,J), AXIS(1,1), AXIS(1,2) )
 
         CALL VCRSS ( AXIS(1,1), AXIS(1,2), AXIS(1,3) )
 
C
C        Can't do the tests if any of the axis are zero.
C
         IF ( ( VNORM ( AXIS(1,1) ) .EQ. 0.D0 )   .OR.
     .        ( VNORM ( AXIS(1,2) ) .EQ. 0.D0 )   .OR.
     .        ( VNORM ( AXIS(1,3) ) .EQ. 0.D0 ) ) THEN
 
                NONZER = .FALSE.
         END IF
 
         IF ( NONZER ) THEN
 
C
C           Find the components of the difference vector w/r to view
C           frame axes.
C
            DO I = 1, 3
 
               POSCMP(I) = VDOT ( POSDIF,AXIS(1,I) ) / VNORM(AXIS(1,I))
 
               VELCMP(I) = VDOT ( VELDIF,AXIS(1,I) ) / VNORM(AXIS(1,I))
 
            END DO
C
C           Divide the downtrack difference in the position by the
C           speed to give the error in time along the flight path.
C
            DELTA  = POSCMP(1) / VNORM ( AXIS(1,1) )
 
C
C           Keep track of the differences.
C
            SUMDT  = SUMDT     + DELTA
 
            SUMADT = SUMADT    + ABS ( DELTA )
 
            DO I = 1, 3
 
               SUMP  (I) = SUMP(I) + POSCMP(I)
 
               SUMV  (I) = SUMV(I) + VELCMP(I)
 
               SUMPA (I) = SUMPA(I) + ABS ( POSCMP(I) )
 
               SUMVA (I) = SUMVA(I) + ABS ( VELCMP(I) )
 
            END DO
 
C
C           For the worst relative difference in the position, record
C           the worst components in view frame coordinates.
C
            IF ( MAXTRU ) THEN
 
               DO I = 1, 3
                  MAXPOS(I)   = POSCMP(I)
                  MAXVEL(I)   = VELCMP(I)
                  WORST(I+2)  = POSCMP(I)
               END DO
 
               MAXADT    = DELTA
               WORST(6)  = DELTA
 
C
C              SIGBIT calculates the number of bits in time the delta
C              time represents.
C
               CALL SIGBIT_1  ( MAXADT, TIMES(J), NLBMAX, NBMAX )
 
            END IF
 
         END IF
 
      END DO
 
 
C
C     Find the average values of all the statistics computed.
C
      AVPOS = POSSUM / DBLE ( NITR )
 
      AVVEL = VELSUM / DBLE ( NITR )
 
      AVRELP = RELPSM / DBLE ( NITR )
 
      AVRELV = RELVSM / DBLE ( NITR )
 
      IF ( NONZER ) THEN
 
         DO I = 1, 3
            AVPO  (I) = SUMP  (I) / DBLE ( NITR )
            AVVE  (I) = SUMV  (I) / DBLE ( NITR )
            AVPOA (I) = SUMPA (I) / DBLE ( NITR )
            AVVEA (I) = SUMVA (I) / DBLE ( NITR )
         END DO
 
         AVDT   = SUMDT  / DBLE ( NITR )
         AVADT  = SUMADT / DBLE ( NITR )
 
      END IF
 
C
C     Write out the test results.
C
C     For the case where we are printing out two numbers on a line
C     convert the dp's to strings so that they will look nice.
C
      CALL DPSTR ( MAXPRL, 14, OUTCH(1) )
      CALL DPSTR ( AVRELP, 14, OUTCH(2) )
      CALL DPSTR ( MAXVRL, 14, OUTCH(3) )
      CALL DPSTR ( AVRELV, 14, OUTCH(4) )
 
      POSSTR = '  Position:             #      # '
 
      CALL REPMC ( POSSTR, '#', OUTCH(1), POSSTR )
      CALL REPMC ( POSSTR, '#', OUTCH(2), POSSTR )
 
      VELSTR = '  Velocity:             #      # '
 
      CALL REPMC ( VELSTR, '#', OUTCH(3), VELSTR )
      CALL REPMC ( VELSTR, '#', OUTCH(4), VELSTR )
 
      WRITE (UNIT,*)
      WRITE (UNIT,*)
      WRITE (UNIT,*) 'Relative differences in state vectors: '
      WRITE (UNIT,*)
      WRITE (UNIT,*) '                              maximum       ' //
     .               '          average'
      WRITE (UNIT,*)
      WRITE (UNIT,*) POSSTR
      WRITE (UNIT,*) VELSTR
      WRITE (UNIT,*)
      WRITE (UNIT,*)
 
      CALL DPSTR ( MAXPMG, 14, OUTCH(1) )
      CALL DPSTR ( AVPOS,  14, OUTCH(2) )
      CALL DPSTR ( MAXVMG, 14, OUTCH(3) )
      CALL DPSTR ( AVVEL,  14, OUTCH(4) )
 
      POSSTR = '  Position (km):        #      # '
 
      CALL REPMC ( POSSTR, '#', OUTCH(1), POSSTR )
      CALL REPMC ( POSSTR, '#', OUTCH(2), POSSTR )
 
      VELSTR = '  Velocity (km/s):      #      # '
 
      CALL REPMC ( VELSTR, '#', OUTCH(3), VELSTR )
      CALL REPMC ( VELSTR, '#', OUTCH(4), VELSTR )
 
      WRITE (UNIT,*) 'Magnitude of state difference vectors:'
      WRITE (UNIT,*)
      WRITE (UNIT,*) '                              maximum       ' //
     .               '          average'
      WRITE (UNIT,*)
      WRITE (UNIT,*) POSSTR
      WRITE (UNIT,*) VELSTR
 
C      WRITE (UNIT,*) 'Relative differences in state vectors: '
C      WRITE (UNIT,*)
C      WRITE (UNIT,*) '                               maximum       '//
C     .               '          average'
C      WRITE (UNIT,*)
C
C      IF ( ( MAXPRL .EQ. 0.D0 ) .OR. ( AVRELP .EQ. 0.D0 ) ) THEN
C         WRITE (UNIT,*) '  Position:             0.00000000000000' //
C     .                  'E+00    0.00000000000000E+00'
C      ELSE
C         WRITE (UNIT,*) '  Position:           ', MAXPRL, AVRELP
C      END IF
C
C      IF ( ( MAXVRL .EQ. 0.D0 ) .OR. ( AVRELV .EQ. 0.D0 ) ) THEN
C         WRITE (UNIT,*) '  Velocity:             0.00000000000000' //
C     .                  'E+00    0.00000000000000E+00'
C      ELSE
C         WRITE (UNIT,*) '  Velocity:           ', MAXVRL, AVRELV
C      END IF
 
 
 
 
 
C      WRITE (UNIT,*) 'Magnitude of state difference vectors:'
C      WRITE (UNIT,*)
C      WRITE (UNIT,*) '                               maximum       '//
C     .               '          average'
C      WRITE (UNIT,*)
C
C      IF ( ( MAXPMG .EQ. 0.D0 ) .OR. ( AVPOS .EQ. 0.D0 ) ) THEN
C         WRITE (UNIT,*) '  Position (km):        0.00000000000000' //
C     .                  'E+00    0.00000000000000E+00'
C      ELSE
C         WRITE (UNIT,*) '  Position (km):      ', MAXPMG, AVPOS
C      END IF
C
C      IF ( ( MAXVMG .EQ. 0.D0 ) .OR. ( AVVEL .EQ. 0.D0 ) ) THEN
C         WRITE (UNIT,*) '  Velocity (km/s):      0.00000000000000' //
C     .                  'E+00    0.00000000000000E+00'
C      ELSE
C         WRITE (UNIT,*) '  Velocity (km/s):    ', MAXVMG, AVVEL
C      END IF
 
      WRITE (UNIT,*)
      WRITE (UNIT,*)
      WRITE (UNIT,*)
 
 
      IF ( NONZER ) THEN
 
         WRITE (UNIT,*) '1) Average components of position difference '
     .               //   'vectors in '
         WRITE (UNIT,*) '   view frame coordinates:'
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  1a) Down track (km):                       ',
     .                     AVPO(1)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  1b) In orbit plane (km):                   ',
     .                     AVPO(3)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  1c) Normal to orbit plane (km):            ',
     .                     AVPO(2)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  1d) Average delta time down track (sec):   ',
     .                     AVDT
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '2) Average |components| of position difference'
     .               //  ' vectors in'
         WRITE (UNIT,*) '   view frame coordinates:'
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  2a) Down track (km):                       ',
     .                                    AVPOA(1)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  2b) In orbit plane (km):                   ',
     .                                    AVPOA(3)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  2c) Normal to orbit plane (km):            ',
     .                                    AVPOA(2)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  2d) Average |delta time| down track (sec): ',
     .                     AVADT
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '3) Components of the position difference ' //
     .                  'vector (in view'
         WRITE (UNIT,*) '   frame coordinates) for the states with ' //
     .                  'the MAXIMUM RELATIVE '
         WRITE (UNIT,*) '   difference in position: '
 
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  3a) Down track (km):                       ',
     .                     MAXPOS(1)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  3b) In orbit plane (km):                   ',
     .                     MAXPOS(3)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  3c) Normal to orbit plane (km):            ',
     .                     MAXPOS(2)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  3d) Delta time down track (sec):           ',
     .                     MAXADT
         WRITE (UNIT,*)
         WRITE (UNIT,*) '  3e) Bits of precision w/r ET time:          '
     .               // '           ', NBMAX
         WRITE (UNIT,*)
         WRITE (UNIT,*)
         WRITE (UNIT,*)
 
      END IF
 
      CALL CHKOUT ( 'STDIFF' )
      RETURN
      END