#  CDPFIT -- fit a polynomial of order npoly to n data points x,y using 
#  orthogonal polynomials technique
#
#  date		author		description
#  12-01-87     J.-C. Hsu	design and coding
#  03-09-90     J.-C. Hsu	rewrite in SPP
#-------------------------------------------------------------------------------

procedure cdpfit (x, y, yerr, npts, wflag, order, ortho, a, sig, chisq)

double 	x[npts], y[npts]	# input: x and y arrays
double	yerr[npts]		# input: y's standard deviations
real	wflag			# input: weighting scheme
int	npts			# input: number of input data points
int	order			# input: order of polynomial
double 	ortho[0:order, 0:order]	# output: orthogonal polynomial construct
double	a[0:order]		# output: coefficients of orthogonal polynomials
				# defined by the input array x
double	sig[0:order]		# output: errors in a
double	chisq			# output: reduced chi-square

pointer	weight			# weights
pointer	work, ycal		# working buffer 
pointer	sp
double	ypin, p2in, chi
int	i, j			# loop index
double	cdpinn()
#------------------------------------------------------------------------------
begin

	call smark (sp)
	call salloc (work, npts, TY_DOUBLE)
	call salloc (weight, npts, TY_DOUBLE)
	call salloc (ycal, npts, TY_DOUBLE)

	# weights are determined by errors
	if (wflag == 0.) {
            do i = 1, npts
            	Memd[weight+i-1] = 1.d0
	} else {
            do i = 1, npts
            	Memd[weight+i-1] = yerr[i] ** wflag
	}

	# construct orthogonal polynomials for this set of x values
        call cdorth (x, Memd[weight], order, npts,     ortho)

	# now calculate inner products and use them to evaluate coefficients
        chi = cdpinn (y, y, Memd[weight], npts)

        do j = 1, npts
            Memd[ycal+j-1] = 0.d0

        do i = 0, order {

	    # evaluate this polynomial at x points and store in work
            call cdopfn (x, Memd[work], i, order, npts, ortho)
            ypin = cdpinn (y, Memd[work], Memd[weight], npts)
            p2in = cdpinn (Memd[work], Memd[work], Memd[weight], npts)
            a[i] = ypin / p2in

	    # calculate functional value
            do j = 1, npts
                Memd[ycal+j-1] = Memd[ycal+j-1] + Memd[work+j-1] * a[i]

	    # 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
	}

	# calculate chi-squared
        chisq = 0.d0
#        do i = 1, npts
#            chisq = chisq + ((Memd[ycal+i-1] - y[i]) / yerr[i]) ** 2

        chisq = chisq / dble(npts - order - 1)

	call sfree (sp)

end