C$Procedure CHEBPL ( Calculate Chebyshev coefficients
C                     and errors of approximation )
      SUBROUTINE  CHEBPL(U,F,INST,D,DM,A)
C
C$ Abstract
C
C     Objective: to build approximation for a vector function F
C     given as a table with vector of nodes U using Chebyshev's
C     polynomials and to determine a set of factors and achived
C     accuracies
C
C$ Declarations

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

      DOUBLE PRECISION    U     ( * )
      DOUBLE PRECISION    F     ( * )
      INTEGER             INST
      DOUBLE PRECISION    D     ( * )
      DOUBLE PRECISION    DM    ( * )
      DOUBLE PRECISION    A     ( * )
      
C$ Copyright
C
C     Copyright (1997-1998), Jet Propulsion Laboratory
C
C$ Required_Reading
C
C
C  Keywords
C
C     CHEBYSHEV
C     POLYNOMIAL APPROXIMATION
C
C
C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     U          I   array of input arguments (time of telemetric
C                    measurements).
C     F          I   array with function values, corresponding to the
C                    input arguments.
C     INST       I   number of instances.
C     D          O   array, containing mean square root errors for
C                    different orders of Chebyshev polynomials.
C     DM         O   array, containing maximum modules of deviations
C                    for different orders of Chebyshev polynomials.
C     A          O   array with factors to Chebyshev polynomials
C$ Detailed_Input
C
C     U          is array, containing the given nodes (arguments, or
C                time of measurements). Values in the array must be
C                ordered, i.e. U (i + 1) > U (i).
C
C     F          is array with given function values, corresponding to
C                the given nodes U (i). This array has one entry:
C                F (i), where i is a number of the node (number of
C                instance from 1 to INST).
C
C     INST       total number of instances, or number of measurements.
C
C
C$ Detailed_Output
C
C     D          array of mean square root deviations for
C                approximation by different orders of Chebyshev
C                polynomials. D(i), where (i) is an order of Ti.
C
C     DM         array of maximum modules of deviations of
C                approximated values from the given values F(i) for
C                different orders of Chebyshev's polynomials. The array
C                has one entry  DM(i), where (i) is the order of Ti.
C
C     A          outcome array with factors for Chebyshev's polynomials
C                (Ti) for approximation of the input function for each
C                order of Ti.
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 you determine factors for Chebyshev's
C     polynomials for approximation of the input function F by
C     different orders of polynomials and corresponding accuracies.
C
C$ Examples
C
C     Suppose you wish to approximate a function given by a
C     tables of arguments (array U) and function values (array F).
C     You know that there are INST argument and function values.
C
C     You can write
C
C        CALL  CHEBPL(U,F,INST,D,DM,A)
C
C$ Restrictions
C
C     None
C
C$ Literature_References
C
C     A.Korn, M.Korn
C     Mathematical Handbook
C     for scientists and engineers. Definitions, Theorems and Formulas
C     For Reference and Review. Second, enlargend and revised edition
C     McGraw-Hill Book Company, 1968
C
C$ Author_and_Institution
C
C
C
C$ Version
C
C-    SPICELIB Version ?.?.?, (to be determined)
C
C-&

C& Index_Entries
C
C     determining a set of factors for Chebushev's polynomial
C     approximation
C
C-&


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


C
C     SPICELIB  functions
C
      DOUBLE PRECISION VDOTG
      
C
C     Parameters from 'ckparam.inc'
C
C     CK4MXD         Max degree of Chebyshev polynomials.
C

C
C     Local parameters
C
C     M - number of maximum Chebyshev's polynomial order + 1
C
C     N - number of "Chebyshev's" nodes
C
C     U1 - begin of subinterval
C
C     U2 - end of subinterval

      INTEGER             M
      PARAMETER         ( M = CK4MXD + 1 )

      INTEGER             N
      PARAMETER         ( N = CK4MXD + 1 )

      DOUBLE PRECISION    XX  ( NUMPNT )
      DOUBLE PRECISION    DF  ( NUMPNT, N )
      DOUBLE PRECISION    R   ( NUMPNT )
      DOUBLE PRECISION    X   ( N )
      DOUBLE PRECISION    TJ  ( N )
      DOUBLE PRECISION    Y   ( N )
      DOUBLE PRECISION    V   ( N )
      DOUBLE PRECISION    W   ( N )
      DOUBLE PRECISION    FC  ( NUMPNT )
      DOUBLE PRECISION    U1
      DOUBLE PRECISION    U2
      DOUBLE PRECISION    WI
      DOUBLE PRECISION    ARG
      DOUBLE PRECISION    DOT1
      DOUBLE PRECISION    DOT2
      DOUBLE PRECISION    A1
      DOUBLE PRECISION    A2
      DOUBLE PRECISION    F1
      DOUBLE PRECISION    F2
      DOUBLE PRECISION    W1
      DOUBLE PRECISION    W2
      DOUBLE PRECISION    DFRES
      DOUBLE PRECISION    DMAX

      INTEGER             I
      INTEGER             IMAX
      INTEGER             J
      INTEGER             J1
      
      DOUBLE PRECISION    PI
      PARAMETER         ( PI = 3.141592653589793D0 )

      LOGICAL                 RETURN


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

      U1 = U(1)
      U2 = U(INST)

C
C     Computing roots of the polynomial with order (M + 1).
C     They represent "the best" nodes for computing the factors
C
      DO J = 1, N
      
         I = N + 1 - J

C
C        Non-dimensioned arguments:
C
         X( J ) = DCOS ( ( 2 * ( I - 1 ) + 1.D0 ) / ( 2 * N ) * PI )

C
C        Dimensioned arguments:
C
         W( J ) = ( X( J ) * ( U2 - U1 ) + ( U2 + U1 ) ) / 2
         
      END DO

C
C     Computing the function values in the "Chebyshev's" nodes
C     using LINEAR interpolation from given input data
C     Set an initial value of the counter J
C
      J = 2
      
      DO I = 1, N

         WI = W( I )
         
         DO WHILE ( U( J ) .LT. WI )
         
            J = J + 1
            
         END DO

C
C     Because of the nature of "Chebyshev's" nodes, all of them are
C     situated inside of the segment [U(1),U(INST)], so we have
C     U( 1 ) < W( 1 ) and (J - 1) > or = 1.
C

C
C     Computing cefficients for linear interpolation
C
         W1 = U( J - 1 )
         W2 = U( J )
         F1 = F( J - 1 )
         F2 = F( J )

         A1 = ( F2 - F1 ) / ( W2 - W1 )
         A2 = F1 - A1 * W1

C
C        Linear interpolation: V = A1*W + A2
C
         V( I ) = A1 * WI + A2

      END DO

C
C     Computing Chebyshev's factors
C
C     Cycle for number of factor (Aj), J=1,2,3,...,M
C
      DO J = 1, M

C
C        Cycle for nodes
C
         DO I = 1, N

            ARG = X( I )

C
C           Computing Chebyshev polynomial values for argument ARG and
C           orders of the polynomials from 0 to M-1:
C
            CALL CNP ( Y, ARG, M-1 )

C
C           TJ(I) is the Chebyshev's polynomial of order J for node I
C
            TJ( I ) = Y( J )

         END DO

C
C        SUM Fi(Xi)*TJ(Xi), I=1,2,3,...,N
C
         DOT1 = VDOTG ( V, TJ, N )

C
C        SUM (TJ(Xi))**2, I=1,2,3,...,N
C
         DOT2 = VDOTG ( TJ, TJ, N )

C
C        COEFFICIENT Aj
C
         A( J ) = DOT1 / DOT2

      END DO

C
C     Computing deviations in the input nodes U(I)
C

C
C     Computing scaled input nodes: U(I) ==> XX(I)
C
      DO I = 1, INST
         XX( I ) = ( 2 * U( I ) - U2 - U1 ) / ( U2 - U1 )
      END DO

C
C     Cycle for scaled nodes XX(I)
C
      DO I = 1, INST
         ARG = XX( I )

C
C        Computing Chebyshev polynomial values for argument ARG and
C        orders of the polynomials from 0 to M-1:
C
         CALL CNP ( Y, ARG, M-1 )

C
C        Cycle for polynomials orders (order is equal J-1)
C
         DO J = 1, M

C
C           Approximated function values
C
            FC( J ) = VDOTG (A, Y, J )

C
C           Errors of approximation
C
            DF( I, J ) = F( I ) - FC( J )

         END DO

      END DO

C
C     Calculation of mean square root deviations for different orders
C     of Chebyshev's polynomials
C

C
C     Cycle for orders (order is equal to J-1)
C
      DO J = 1, M

         DFRES = 0.D0

         DO J1 = 1, INST

            DFRES = DFRES + DF( J1, J ) * DF( J1, J )

         END DO

C
C        D(J) = DSQRT(VDOTG(DF(1,J),DF(1,J),INST))/INST
C
         D( J ) = DSQRT ( DFRES ) / INST

      END DO

C
C     Calculation of maximum module of deviation from the given input
C     function values
C
      DO J = 1, M

         DO  I = 1, INST

            R( I ) = DABS ( DF( I, J ) )

         END DO

C
C        Calculate the worst deviation
C
         CALL MAXAD ( R, INST, DMAX, IMAX )

C
C        The worst deviation for the order (J-1)
C
         DM( J ) = DMAX

      END DO
      
      CALL CHKOUT ( 'CHEBPL' )
      RETURN
      END
