procedure polylq (x, y, sigmay, npts, nterms, mode, result, chisqr)
 
#    polylq -- Performs Least Squares fit to a polynomial
# 
#    Description:
#    ------------
#    This is a general Polynomial least squares fitting subroutine.
#    It is based on the program POLFIT shown in Bevington (p. 140ff).
#    The returned parameters are an array of the fitted polynomial
#    coefficients and the reduced Chi-squared. 
#    The value of Chi-squared will be zero for a degenerate fit.
# 
#
#    Subroutines, Functions called:
#    ------------------------------
#    R = DETERM (ARRAY, SIZE, NORDER) ... Returns determinant of subarray
#                                        NORDER squared from ARRAY which
#                                        is really dimensioned SIZE.
#                                        SQUARE ARRAYS ONLY!!!!!!!!!!!!!
#
#        Version    Date        Author     Description
#	    	 27-JAN-1993  R. Williamson Converted to SPP
#           6    29-NOV-1984  Z.G. Levay   Updated PDL
#           5     1-NOV-1983  PAT MURPHY   Submitted for inspection
#           4    22-SEP-1983  PAT MURPHY   Put in Bevington's code
#           3    19-SEP-1983  PAT MURPHY   Insert Fortran code
#           2    19-SEP-1983  PAT MURPHY   Insert PDL
#           1    19-SEP-1983  PAT MURPHY   Prolog Generation

int npts    	    	    	    #  (I) Number of points in data arrays
int nterms  	    	    	    #  (I) Degree of the fit
int mode    	    	    	    #  (I) Mode of sigma handling
int nmax, i, j, k, l, n	    	    #  (L) all local variables

double x[npts], y[npts]	    	    #  (I) Arrays to be fitted
real sigmay[npts]   	    	    #  (I) Sigma of these
double result[nterms]	    	    #  (O) Result of the fit
double chisqr	    	    	    #  (O) Chi-squared
double weight	    	    	    #  (L) local variable
double xi, yi	    	    	    #  (L) Local variables
double determ()    	    	    #  (S) Determinant function
double sumx[19], sumy[19], xterm    
double yterm, array[10,10], chisq  
double delta

begin

	#  if there aren't enough points or
        if (npts   < nterms) {
        call error (1,"Not enough points for fit")
        }

#  the number of terms is not in the
#  range 2-10, we want to quit now.
        if (nterms <= 1  || nterms > 10)    {
        call error (1, "Not the right number of terms for fit")
        }



        chisq = 0.
        nmax = 2 * nterms - 1

        do i = 1, nmax {
                sumx[i] = 0.0d0
                sumy[i] = 0.0d0
        }

#       Figure out the weights to apply for each point

        do i = 1, npts {
          xi = x[i]
          yi = y[i]
          weight = 1.
          if (mode < 0 && yi != 0.) weight = 1./yi
          if (mode > 0 && sigmay(i) != 0.) weight = 1./sigmay[i]
          xterm = weight
          do j = 1, nmax {
            sumx [j] = sumx [j] + xterm
            xterm = xterm * xi
          }
          yterm = weight * yi
          do j = 1, nterms {
            sumy [j] = sumy [j] + yterm
            yterm = yterm * xi
          }
          chisq = chisq + weight * yi ** 2
        }

#       Now construct Matrices

        do j = 1, nterms {
          do k = 1, nterms {
            n = j + k - 1
            array[j,k] = sumx[n]
          }
        }

#       Get determinant of matrix

        delta = determ (array, 10, nterms)
        if (delta == 0.) {
        call eprintf("ERROR: Determinant of matrix is zero ")
          do j = 1, nterms {
            result[j] = 0.0d0
          }
          chisqr = 0.
        }

#       Calculate the polynomial coefficients

        do l = 1, nterms {
          do j = 1, nterms {
            do k = 1, nterms {
              array [j,k] = sumx[j + k - 1]
            }
            array [j,l] = sumy [j]
          }
          result [l] = determ (array, 10, nterms) / delta
        }

#       Now calculate chi squared

        chisqr = 0.
        if (nterms != npts) {
          do j = 1, nterms {
            chisq = chisq - 2. * result[j] * sumy[j]
            do k = 1, nterms {
              chisq = chisq + result[j] * result[k] * sumx [j + k - 1]
            }
          }
          chisqr = chisq / (npts - nterms)
        }

end