SUBROUTINE CHEBYS ( U,F,ICOUNT,IVEL,NSUB0,NSTEP0,TOL,ITYPE,
     .                                            NE,IS,B,NCI,IFLAG)
C$ Abstract
C
C     Objective: to divide given interval into several subintervals
C     taking into account the tolerances for given criterion of
C     accuracy of approximation by Chebyshev's polynomials and
C     determine a corresponding set of factors
C
C$ Declarations
C
      IMPLICIT              NONE
      
      INCLUDE               'ckparam.inc'
      INCLUDE               'makeck4.inc'

      DOUBLE PRECISION      U     ( * )
      DOUBLE PRECISION      F     ( QAVSIZ, * )
      INTEGER               ICOUNT
      INTEGER               IVEL
      INTEGER               NSUB0
      INTEGER               NSTEP0
      DOUBLE PRECISION      TOL   ( * )
      INTEGER               ITYPE ( * )
      INTEGER               NE    ( * )
      INTEGER               IS
      DOUBLE PRECISION      B     ( QAVSIZ, NUMSUB, * )
      INTEGER               NCI   ( NUMSUB, * )
      INTEGER               IFLAG
C
C$ Copyright
C
C     Copyright (1997-1998), Jet Propulsion Laboratory
C
C$ Required_Reading
C
C     CK.REQ
C
C  Keywords
C
C     SUBINTERVAL
C     CHEBYSHEV
C    POLYNOMIAL APPROXIMATION
C

C
C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     U          I   array with given nodes (arguments).
C     F          I   array with given function values, corresponding
C                    to the given nodes.
C     ICOUNT     I   number of instances.
C     IVEL       I   angular velocities flag:
C                    IVEL = 1 if input data contains angular
C                    velocities,
C                    IVEL = 0 when input data doesn't contain angular
C                    velocities.
C     NSUB0      I   initial length of subinterv.(number_of_nodes-1).
C     NSTEP0     I   increment for increasing/decreasing subinterval
C                    length.
C     TOL        I   array of required tolerances.
C     ITYPE      I   indication of a tolerance criteria:
C                    ITYPE = 1 for mean square root criterion,
C                    ITYPE = 2 for module of maximum deviation.
C     IS         O   number of subintervals.
C     NE         O   array with numbers of edge nodes of
C                    subintervals.
C     B          O   array with Chebyshev's polynomials factors.
C     NCI        O   array with recommended numbers of Chebyshev's
C                    polynomials factors for each subinterval and
C                    each parameter.
C     IFLAG      O
C$ Detailed_Input
C
C     U          is array, containing the given nodes (arguments, or
C                time of measurements), extracted from input
C                file. Time format do not matter (tics, seconds,
C                etc.). The only requirement is that values must be
C                ordered, i.e. U (i + 1) > U (i).
C
C     F          is array with given function values, corresponding
C                to the given nodes U(i). This array has two entries:
C                F j,i), where j is a number of parameters underlying
C                to approximation (j = 4 for Q0,Q1,Q2,Q3 or j=7 when
C                data for angular velocities are included: Q0,Q1,Q2,
C                Q3,A1,A2,A3); i is a number of node (number of
C                instance from 1 to ICOUNT).
C
C     ICOUNT     total number of nodes, or number of measurements.
C
C     IVEL       indication of presence or absence of data about
C                angular velocities in the corresponding source file:
C                IVEL = 1, when data contains angular velocities;
C                one instance contain seven parameters (Q0,Q1,Q2,Q3,
C                A1,A2,A3);
C
C                IVEL = 0, when data doesn't contain angular
C                velocities; one instance contains four parameters
C                (Q0,Q1,Q2,Q3).
C
C     NSUB0      initial length of subinterval (number of nodes -1).
C                Subroutine CHEBYSHEV starts analysis for required
C                length of each immediate subinterval taking i0 as a
C                beginning of the subinterval and (i0 + NSUB0) as the
C                end of the subinterval. So, amount (NSUB + 1) can be
C                considered as expected typical reasonable length of
C                the subinterval.
C
C     NSTEP0     increment for increasing/decreasing subinterval length
C                If you have more or less modern computer we reccommend
C                to take NSTEP = 1. This will provide the best data

C                compression.
C
C                At the beginning of analysis of the immediate
C                subinterval it started at i0 and ended at (i0 + NSUB).
C                If the accuracy is not corresponding to a given
C                tolerance, the end of the subinterval must be shifted
C                to decrease the length of the immediate subinterval.
C                The end of the subinterval will be changed to
C                (i0 + NSUB0 - NSTEP). If the accuracy still will be
C                more then given tolerance, the right end of the
C                subinterval will be shifted to (i0 + NSUB0 -2*NSTEP)
C                and so on until the accuracy will fit to the tolerance
C
C                If the accuracy at the beginning of analysis of the
C                immediate subinterval will be better than requested by
C                the tolerance, than the end of the subinterval will be
C                shifted to (i0 + NSUB0 + NSTEP). If the accuracy still
C                will be better, than requested, the right end of the
C                subinterval will changed to (i0 + NSUB0 + 2*NSTEP) and
C                so on.
C
C     TOL        array of acceptable tolerances for each of quaternion
C                components. This array must contain four or seven
C                values correspondingly for IVEL = 0 or IVEL = 1.
C
C     ITYPE      indication of a tolerance nature. We considered two
C                reasonable criteria: mean square root deviation and
C                maximum module of deviation for all nodes.
C
C                For mean square root criterion ITYPE = 1, and
C
C                              --
C                              \
C                error =     ( /  (Fi - Gi)**2  )/N,
C                              --
C                            i=1,2,...,N
C
C                where Fi - given function value in the i-node,
C                      Gi - approach to the Fi,
C                      N  - number of nodes in the subinterval.
C
C
C                In case of maximum module of deviation critaria
C                ITYPE = 2:
C
C                 error =    MAX |Fi - Gi|.
C                          i=1,2,...,N
C
C$ Detailed_Output
C
C     NE         array with numbers of selected edge points for
C                subintervals. The first point is begin of the interval
C                and begin of the first subinterval. The second point
C                is the end of the first subinterval and the beginning
C                of the second subinterval. The third point is the end
C                of the second subinterval and the beginning of the
C                third subinterval and so on. The last point in the
C                array is the end of the last subinterval corresponding
C                to the end of the interval.
C
C     B          outcome array with factors for Chebyshev's polynomials
C                for approximation of the input functions for each
C                subinterval and each quaternion component.
C                This array has three entries (i,j,k): number of the
C                quaternion component, number of the subinterval and
C                number of the factor for corresponding Chebyshev's
C                polynomial. Array B contains complete set of factors,
C                corresponding to maximum available polynomial's
C                degree.
C
C     NCI        array with recommended numbers of factors to
C                Chebyshev's polynomials. This array has two entries:
C                number of subinterval and number of the quaternion's
C                component.
C                NCI is intended for extraction from the array B only
C                reasonable number of coefficients for adequate
C                representation of the quaternions components (i.e.
C                with accordance to the given tolerance).
C
C     IFLAG
C
C$ Parameters
C
C     None.
C
C
C$ Exceptions
C
C     None.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine allows to construct from the given interval
C     an optimal subintervals that approximates the functions (i.e.
C     quaternion components) by a set of the factors to Chebyshev's
C     polynomials; minimum number of factors provided required accuracy
C     will be determined. This allows to minimize required memory for
C     an adequate representation of data about pointing.
C
C$ Examples
C
C     Suppose you wished biuld approximation to pionting data,
C     represented by quaternions and angular velocities in CK-file
C     (IVEL = 1).
C     First of all you have to fill out the array U with values of time
C     and array F with quaternion components and angular velocity
C     values, corresponding to the nodes in the array U. Be sure that
C     the values in the array U are ordered and all instances are
C     presented, from 1 to ICOUNT. Suppose you estimated the size of
C     the subintervals as equal to 51 nodes (ie. NSUB0 = 51 - 1 = 50),
C     suggested that the step for adjusting length of a subinterval is
C     equal 1. Suppose that you want consider as a criterion for
C     accuracy a maximum module of deviation (ITYPE = 2) and you only
C     want control the accuracy for Q0,Q1,Q2,Q3, not for A1,A2 and A3.
C     For instance, you selected for the tolerance 0.01 (the same for
C     all quaternion components); for the rest of parameters (A1,A2,A3)
C     you have to select just big amount, for instance, 1000.
C     So, you can write
C        DO I = 1,4
C        TOL(I) = 0.01
C        ENDDO
C
C        DO I = 5,7
C        TOL(I) = 1000.
C        ENDDO
C
C        CALL  CHEBYSHEV(U,F,ICOUNT,1,50,1,TOL,1, NE,IS,B,NCI)
C
C$ Restrictions
C
C     None
C
C$ Literature_References
C
C     None
C
C$ Author_and_Institution
C
C
C$ Version
C
C-    SPICELIB Version ?.?.?, (to be determined)
C
C-&

C& Index_Entries
C
C     subintervals with Chebushev's polynomial approximation
C
C-&

C$ Revisions
C
C-    Beta Version
C
C         Examples section completed.
C-&

C
C     SPICELIB  functions
C

C
C     Parameters from 'ckparam.inc'
C
C     CK4MXD         Max degree of Chebyshev polynomials.
C
C     QSIZ           count of the quaternion components.
C
C     QAVSIZ         count of the quaternion and angular rate 
C                    components together.
C

C
C     Local parameters
C
C     NCOEFF         Max number of coefficients in each polynomial.
C
C     SETLIM         Max available number of coefficients sets for a
C                    one call of subroutine CHEBYSHEV
C
C
      INTEGER               NCOEFF
      PARAMETER           ( NCOEFF = CK4MXD + 1 )

      INTEGER               SETLIM
      PARAMETER           ( SETLIM = NUMSUB )

C
C     Local variables
C
      DOUBLE PRECISION      D     ( QAVSIZ, NCOEFF )
      DOUBLE PRECISION      DM    ( QAVSIZ, NCOEFF )
      DOUBLE PRECISION      C     ( QAVSIZ, NCOEFF )
      INTEGER               NC    ( QAVSIZ )
      INTEGER               I0
      INTEGER               IEND
      INTEGER               I
      INTEGER               NSUB
      INTEGER               NSTEP
      INTEGER               J
      INTEGER               NPAR

      LOGICAL               ACCURA
      LOGICAL               RETURN


      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'CHEBYS' )
      END IF

C
C     Set the number of parameters underlying to approximation
C
      IF ( IVEL .EQ. 0 ) NPAR = QSIZ
      IF ( IVEL .EQ. 1 ) NPAR = QAVSIZ

C
C     Set initial begin of the subinterval I0
C
      I0 = 1

C
C     Set initial value for subintervals counter IS
C
      IS = 1
      IEND = 0

C
C     Set initial value of the flag, characterizing the state of
C     immediate subinterval determination state (IFLAG)
C
      IFLAG = 0

C
C     Set initial value of NSUB0
C
      NSUB = NSUB0

C
C     Set initial increment for increasing/decreasing subinterval
C     length
C
      NSTEP = NSTEP0

C
C     Counter for number of subinterval
C

C
C     Check if NSUB is too much big
C
      IF ( ( I0 + NSUB ) .GT. ICOUNT ) THEN

         NSUB = ICOUNT - I0

      END IF

      DO WHILE ( IEND .NE. 1 )

C
C        Call the routine for determination of the given subinterval
C        set of parameters
C
         CALL SUBINT ( I0, NSUB, U, F, IVEL, D, DM, C)

C
C        Recording values of factors to Chebyshev's polynomials
C

C
C        Cycle for number of pointing parameters NPAR (J)
C
         DO J = 1, NPAR

C
C              Cycle for number of factor to Chebyshev's polynomial (I)
C
               DO I = 1, NCOEFF

                  B( J, IS, I ) = C( J, I )

               END DO

         END DO

C
C        Check if the current subinterval is done (IFLAG = 3);
C        If Yes, to do the next one!
C
         IF ( ( IFLAG .EQ. 3 ) .OR. ( IFLAG .EQ. 5 ) ) THEN

            NE( IS ) = I0

C
C           Determination of the reasonable numbers of factors:
C
            CALL FILTER ( ITYPE, IVEL, TOL, D, DM, NC )

            DO J = 1, NPAR
               NCI( IS, J ) = NC( J )
            END DO

C
C           Reserved numbers of subintervals is full! Terminate
C           calculation and return to main program
C
            IF ( IFLAG .EQ. 5 ) THEN

               NE( IS + 1 )= I0 + NSUB
               IEND = 1

            ELSE

C
C              Shift initial number of the next subinterval begin
C
               I0 = I0 + NSUB

C
C              Set NSUB to initial value NSUB0 for processing the
C              next subinterval
C
               NSUB = NSUB0

C
C              Check if determined begin of the next interval not
C              exceed ICOUNT. If yes, then correct it
C
               IF ( ( I0 + NSUB ) .GT. ICOUNT ) THEN

                  NSUB = ICOUNT - I0

               END IF

C
C              Set initial value IFLAG for next interval
C
               IFLAG = 0

C
C              Go to computing parameters of the next subinterval
C
               IS = IS + 1

            END IF

         ELSE

C
C           Check, if the accuracy for maximum polynom. degree is OK?
C
            IF ( ACCURA ( ITYPE, IVEL, TOL, D, DM ) ) THEN

C
C              Is the interval processing complete (IFLAG = 4)?
C
               IF ( IFLAG .EQ. 4 ) THEN

C
C                 work can be completed
C

C
C                 Record the begin and the end nodes for the interval
C
                  NE( IS ) = I0
                  NE( IS + 1 ) = ICOUNT

C
C                 Determination of the reasonable numbers of factors
C                 for last subinterval:
C
                  CALL FILTER ( ITYPE, IVEL, TOL, D, DM, NC)

C
C                 Recording reccommended numbers of factors
C
                  DO J = 1, NPAR

                     NCI( IS, J ) = NC( J )

                  END DO

C
C                 Recording values of factors to Cheby polynomials
C
                  DO  J = 1, NPAR

                     DO I = 1, NCOEFF

                        B( J, IS, I ) = C( J, I )

                     END DO

                  END DO

                  IEND = 1
 
               ELSE

C
C                 Increasing the subinterval length
C
                  IF ( IFLAG .EQ. 2 ) THEN

                     NSTEP = NSTEP / 3 + 1

                  END IF

                  NSUB  = NSUB + NSTEP

C
C                 Check if I0+NSUB exceeds  ICOUNT
C
                  IF ( ( I0 + NSUB ) .GE. ICOUNT ) THEN

C
C                    correction of the NSUB
C
                     NSUB = ICOUNT - I0

C
C                    set IFLAG=4, i.e. "processing of subinterval must
C                    be completed"
C
                     IFLAG = 4

                  ELSE

C
C                    Change of a flag, characterized current state of
C                    the subinterval processing to IFLAG = 1, i.e.
C                    "we are expanding the immediate interval"
C
                     IFLAG = 1

                  END IF

               END IF

            ELSE

C
C              Accuracy is insufficient, so the immediate subinterval
C              must be less in size; correction NSUB:
C
               IF ( IFLAG .EQ. 1 .OR. IFLAG .EQ. 4 ) THEN

                  IF ( NSTEP .EQ. 1 ) THEN

                     NSUB = NSUB - NSTEP
                     IFLAG = 3    

C
C                    End of value reserved in
C                    array description
C
                     IF ( IS .GE. SETLIM - 1 ) THEN

                        IFLAG= 5

                     END IF

                     NSTEP = NSTEP0

                  ELSE

                     NSTEP = NSTEP / 3 + 1

                     IF ( NSUB .LE. NSTEP ) THEN

                        NSTEP = NSUB / 3 + 1

                     END IF

                     NSUB = NSUB - NSTEP

                     IFLAG = 2

                  END IF

               ELSE

                  IF ( NSUB .LE. NSTEP ) THEN

                     NSTEP = NSUB / 3 + 1

                  END IF

                  NSUB = NSUB - NSTEP
                  IFLAG = 2

               END IF

            END IF

         END IF

C
C        Compute again al subinterval parameters for changed
C        subinterval length.
C
      END DO

      CALL CHKOUT ( 'CHEBYS' )
      RETURN
      END