#  CDORTH -- calculate the orthogonal polynomials coefficients for the grid x
#
#  date		author		description
#  12-01-87	J.-C. Hsu	design and coding
#  03-09-90	J.-C. Hsu	rewrite in SPP
#-------------------------------------------------------------------------------

procedure cdorth (x, weight, order, npts, ortho)

double 	x[npts]			# input: independent variable array
double	weight[npts]		# input: weights
int	order			# input: order of the fitting polynomial
int	npts			# input: number of input points
double 	ortho[0:order, 0:order]	# output: orthogonal polynomial construct

pointer	sp, z
double 	term
int	i, j, k			# loop indices
#------------------------------------------------------------------------------
begin
	call smark (sp)
	call salloc (z, 2*order+1, TY_DOUBLE)

#  first put the q's in ortho
#  calculate the q's
#  first calculate z(k) = norm for x**k , k = 0,... 2 * order

        do k = 0, 2*order
            Memd[z+k] = 0.d0

        do i = 1, npts {
            term = weight[i]
            do k = 0, 2*order {
                Memd[z+k] = Memd[z+k] + term
                term = term * x[i]
  	    } 
   	}

        do k = 0, 2*order
            Memd[z+k] = Memd[z+k] / double(npts)

#  now use recursion relations to calculate the q's
        do i = 0, order {
            do k = 0, i {
                ortho[k, i] = Memd[z+k+i]
                do j = 0, k-1
                    ortho(k, i) = ortho(k, i) - ortho(j, k) *
                                       ortho(j, i) / ortho(j, j)
	    }
	}

#  now convert from q's to coefficients
        do i = 0, order {
            do k = 0, i-1
                ortho(k, i) = -ortho(k, i) / ortho(k, k)
	}

	call sfree(sp)

end