        SUBROUTINE CDPFIT(
*
*  inputs
*
     :                  X, Y, YERR, NPTS, WFLAG, ORDER, NPMAX,
*
*  outputs
*
     :                  ORTHO, A, SIG, CHISQ)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  fit a polynomial of order npoly to n data points x,y using orthogonal
*  polynomials technique
*
*  Description:
*  ------------
*  maximum number of input data points is 2000
*
*  FORTRAN name: CDPFIT.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       None
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       None
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       12-01-87     J.-C. HSU        design and coding
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --X and Y arrays
        DOUBLE PRECISION        X(1), Y(1),
*                                       --Y's standard deviations
     :                          YERR(1)
*                                       --weighting scheme
        REAL                    WFLAG
*                                       --number of input data points
        INTEGER                 NPTS,
*                                       --order of polynomial
     :                          ORDER,
*                                       --maximum order
     :                          NPMAX
*
*== output:
*                                       --orthogonal polynomial construct
        DOUBLE PRECISION        ORTHO(0:NPMAX, 0:NPMAX),
*                                       --coefficients of orthogonal polynomials
*                                       --define by X
     :                          A(0:NPMAX),
*                                       --errors in A
     :                          SIG(0:NPMAX),
*                                       --reduced chi-square
     :                          CHISQ
*
*== local:
*                                       --weights
        DOUBLE PRECISION        WEIGHT(2000),
*                                       --working buffer
     :                          WORK(2000), YCAL(2000), YPIN, P2IN,
     :                          CHI, CDPINN
*                                       --loop index
        INTEGER                 I, J
*
*------------------------------------------------------------------------------
*
*  weights are determined by errors
*
        DO 10 I = 1, NPTS
            WEIGHT(I) = YERR(I) ** WFLAG
   10   CONTINUE
*
*  construct orthogonal polynomials for this set of x values
*
        CALL CDORTH (X, WEIGHT, ORDER, NPMAX, NPTS,     ORTHO)
*
*  now calculate inner products and use them to evaluate coefficients
*
        CHI = CDPINN (Y, Y, WEIGHT, NPTS)
*
        DO 20 J = 1, NPTS
            YCAL(J) = 0.D0
   20   CONTINUE
*
        DO 40 I = 0, ORDER
*
*  evaluate this polynomial at x points and store in work
*
            CALL CDOPFN (X, I, NPMAX, NPTS, ORTHO,      WORK)
            YPIN = CDPINN (Y, WORK, WEIGHT, NPTS)
            P2IN = CDPINN (WORK, WORK, WEIGHT, NPTS)
            A(I) = YPIN / P2IN
*
*  calculate functional value
*
            DO 30 J = 1, NPTS
                YCAL(J) = YCAL(J) + WORK(J) * A(I)
   30       CONTINUE
*
*  sigma given by change in a which gives change in (non-reduced) chisq of 1
*
            SIG(I) = SQRT (1.D0 / P2IN)
*
*  reduced chi-square from fit of polynomial of degree i
*
            CHI = CHI - P2IN * A(I) ** 2
   40   CONTINUE
*
*  calculate chi-squared
*
        CHISQ = 0.D0
        DO 50 I = 1, NPTS
            CHISQ = CHISQ + ((YCAL(I) - Y(I)) / YERR(I)) ** 2
   50   CONTINUE
        CHISQ = CHISQ / DBLE(NPTS - ORDER - 1)
*
        RETURN
        END
