C$Procedure MKCOEF (call Chebyshev coefficients calculation
C                   and write CK4 data to segment of result DAF file)

      SUBROUTINE MKCOEF ( ICOUNT, TOL, ITYPE, U, F, IVEL,
     .                                             HANDL1, NUMINT )

C$ Abstract
C
C      Subroutine MKCOEF calls subroutine Chebsh
C     to calculate Chebyshev polynomials coefficients.
C      Then it puts sets of coefficients and reference times
C     to segment of result DAF file
C
C$ Declarations


      IMPLICIT                NONE
      
      INCLUDE                 'ckparam.inc'
      INCLUDE                 'makeck4.inc'

      INTEGER                 ICOUNT
      DOUBLE PRECISION        TOL( * )
      INTEGER                 ITYPE
      DOUBLE PRECISION        U( NUMPNT )
      DOUBLE PRECISION        F( QAVSIZ, NUMPNT )
      INTEGER                 IVEL
      INTEGER                 HANDL1
      INTEGER                 NUMINT

C$ Copyright
C
C      Copyright (1997-1998), Jet Propulsion Laboratory
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     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     HANDL1     I   Handle for output DAF file
C     NUMINT    I/O  Number of subintervals of approximation
C
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     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
C     HANDL1     Handle of output DAF-file to write data for
C                for generic segments
C
C     NUMINT     Integrated value of subintervals number
C
C
C


      INTEGER                 NDPVAL
      PARAMETER             ( NDPVAL = CK4RSZ )

      INTEGER                 NPKTS
      PARAMETER             ( NPKTS = 1 )

      INTEGER                 COUNT
      INTEGER                 I
      INTEGER                 ICOUN1
      INTEGER                 IFLAG
      INTEGER                 IL
      INTEGER                 IK
      INTEGER                 IQ
      INTEGER                 IS
      INTEGER                 LASTAD
      INTEGER                 NCI( NUMSUB, QAVSIZ )
      INTEGER                 NCOEF
      INTEGER                 NE( NUMSUB )
      INTEGER                 NSTEP
      INTEGER                 NSUB0
      INTEGER                 PKTSIZ ( NPKTS )
      INTEGER                 PSIZ
      INTEGER                 JJ

      DOUBLE PRECISION        B( QAVSIZ, NUMSUB, CK4MXD + 1 )
      DOUBLE PRECISION        PKTDAT ( NDPVAL )
      DOUBLE PRECISION        SCLKDP ( NDPVAL )
      DOUBLE PRECISION        VALUE

      LOGICAL                 RETURN


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

      PSIZ = QAVSIZ

      LASTAD = 9

      IFLAG = 5

      DO WHILE ( IFLAG .EQ. 5 )

C
C        Initial lenght of subinterval
C
         NSUB0 = ICOUNT / 2 + 1

C
C        Initial step for INCR/DECR subinterval lenght
C
         NSTEP = NSUB0 /  3 + 1

C
C        Approximation is impossible if subinterval
C        contains ZERO points
C
         IF ( ICOUNT .EQ. 0 ) THEN

            CALL TOSTDO ( 'Done' )
            CALL CHKOUT ( 'MKCOEF' )
            RETURN

         END IF

         IF ( ICOUNT .EQ. 1 ) THEN

            ICOUN1 = 2

            U( 2 ) = U( 1 ) + 1.D-1
            U( 1 ) = U( 1 ) - 1.D-1

            DO IK = 1, QAVSIZ

               F(IK, 2) = F(IK, 1)

            END DO

         ELSE

            ICOUN1 = ICOUNT

         END IF

C         WRITE(*,*) 'DATA:', U(1), U(ICOUNT)

         CALL CHEBYS ( U, F, ICOUN1, IVEL, NSUB0, NSTEP, TOL,
     .                   ITYPE, NE, IS, B, NCI, IFLAG )

         DO IL = 1, IS

C
C           Assign mid time of subinterval to PKTDAT
C
            VALUE = U( NE( IL ) ) + ( U( NE( IL + 1 ) )
     .                                       - U( NE( IL ) ) ) / 2.d0
            PKTDAT (1) = VALUE

C
C           Assign radii of subinterval to PKTDAT
C

            VALUE = ( U( NE( IL + 1 ) ) - U( NE( IL ) ) ) / 2.d0

            PKTDAT (2) = VALUE


C
C           Pack degrees into record. Compute total
C           number of element in the record.
C
            PKTSIZ (1) = LASTAD

            DO IK = 1, QAVSIZ

               PKTDAT ( IK + 2 ) = NCI( IL, IK )
               PKTSIZ (1) = PKTSIZ (1) + NCI( IL, IK )

            END DO

C
C           Output sets of coefficients for Q0, Q1, Q2, Q3
C           and  for AV1, AV2, AV3 to
C           output DAF file
C
            COUNT = LASTAD

            DO IQ= 1,PSIZ

               DO NCOEF= 1, NCI( IL, IQ )

                  COUNT = COUNT + 1

                  PKTDAT( COUNT ) = B( IQ, IL, NCOEF )

               END DO

C
C              End of cycle for output data
C
            END DO

C
C           Compute start of the interval
C
            SCLKDP (1) = PKTDAT( 1 ) - PKTDAT( 2 )

C
C           Add data packets to the segment of DAF file
C
C            WRITE(*,*) PKTDAT(1)-PKTDAT(2), PKTDAT(1)+PKTDAT(2)
            CALL CKW04A ( HANDL1, NPKTS, PKTSIZ, PKTDAT, SCLKDP )

         END DO

         NUMINT = NUMINT + IS

         IF ( IFLAG .EQ. 5 ) THEN

            DO I = NE( NUMSUB ), ICOUNT

               U( I - NE( NUMSUB ) + 1 ) = U( I )
               
               DO JJ = 1, QAVSIZ
               
                  F( JJ, I - NE( NUMSUB ) + 1 ) =  F( JJ, I )
                 
               END DO

            END DO

            ICOUNT = ICOUNT - NE( NUMSUB ) + 1

         END IF

      END DO

      CALL TOSTDO ( 'Done' )
      CALL CHKOUT ( 'MKCOEF' )
      RETURN

      END

