SUBROUTINE CDORTH (
*
*  inputs
*
     :                  X, WEIGHT, ORDER, NPMAX, NPTS,
*
*  outputs
*
     :                  ORTHO)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  calculate the coefficients for orthogonal polynomials for grid x
*
*  Description:
*  ------------
*  maximum order of polynomial is 20.
*
*  FORTRAN name: CDORTH.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:
*                                       --input X array
        DOUBLE PRECISION        X(1),
*                                       --weights
     :                          WEIGHT(1)
*                                       --order of fitting polynomial
        INTEGER                 ORDER,
*                                       --maximum order
     :                          NPMAX,
*                                       --number of input points
     :                          NPTS
*
*== output:
*                                       --orthogonal polynomial construct
        DOUBLE PRECISION        ORTHO(0:NPMAX, 0:NPMAX)
*
*== local:
*                                       --work array
        DOUBLE PRECISION        Z(0:40), TERM
*                                       --loop indices
        INTEGER                 I, J, L
*
*------------------------------------------------------------------------------
*
*  first put the q's in ortho
*
*  calculate the q's
*
*  first calculate Z(L) = norm for X**L , L = 0,... 2 * ORDER
*
        DO 10 L = 0, 2*ORDER
            Z(L) = 0.D0
   10   CONTINUE
*
        DO 30 I = 1, NPTS
            TERM = WEIGHT(I)
*
            DO 20 L = 0, 2*ORDER
                Z(L) = Z(L) + TERM
                TERM = TERM * X(I)
   20       CONTINUE
   30   CONTINUE
*
        DO 40 L = 0, 2*ORDER
            Z(L) = Z(L) / DBLE(NPTS)
   40   CONTINUE
*
*  now use recursion relations to calculate the q's
*
        DO 70 I = 0, ORDER
            DO 60 L = 0, I
                ORTHO(L, I) = Z(L+I)
*
                DO 50 J = 0, L-1
                    ORTHO(L, I) = ORTHO(L, I) - ORTHO(J, L) *
     :                                  ORTHO(J, I) / ORTHO(J, J)
   50           CONTINUE
   60       CONTINUE
   70   CONTINUE
*
*  now convert from q's to coefficients
*
        DO 90 I = 0, ORDER
            DO 80 L = 0, I-1
                ORTHO(L, I) = -ORTHO(L, I) / ORTHO(L, L)
   80       CONTINUE
   90   CONTINUE
*
        RETURN
        END