c Demonstration program for data fitting (uses a polynomial)
      program polyfit

      logical fail

      parameter (maxpts=50,maxco=10) ! array sizes
      parameter (eps=1.0e-6) ! numerical bandwidth

c     main data and parameter arrays
      dimension x(maxpts),y(maxpts),w(maxpts),yco(maxpts)

c     working space
      dimension wrkptsco(maxpts,maxco)

c     read data
      read *,npts
      if (npts.gt.maxpts) stop 'Too many points!'
      do i=1,npts
         read *,x(i),y(i),w(i) ! ordinate, data, weight
         yco(i) = y(i) ! copy to fitting array
      enddo

c     choose number of parameters for fitting polynomial
      read *,nco ! coefficients in polynomial (= order + 1)
      if (nco.gt.maxco) stop 'Polynomial order too large!'

c     call fitting routine: polynomial: c1 + c2.x + c3.x^2 + c4.x^3 + ...
c     supply y-values in yco, returns coefficients in first nco elements
      call bm_polyfitg(x,w,npts,nco,eps,yco,fail,res,wrkptsco)
      if (fail) print *,'Singular matrix!'

c     print coefficients
      print *,'polynomial coefficients (increasing powers of x): '
      print *,(yco(i),i=1,nco)
      print *,'residue=',res

c     graphics - uncomment calls if you have the Bubblesoft Graphics Lubrary
      read *,xmin,xmax,ymin,ymax ! display limits
      call sk_init
      call sk_simplegrid(xmin,xmax,ymin,ymax)

c     draw polynomial
      read *,nlines ! polynomial lines
      dx=(xmax-xmin)/nlines
      x1=xmin
      y1=calcpoly(nco,yco,x1)
      call sk_setforecol(0.0,1.0,1.0) ! turquoise
      do i=1,nlines
         x2=xmin+i*dx
         y2=calcpoly(nco,yco,x2)
         call sk_line(x1,y1,x2,y2)
         x1=x2
         y1=y2
      enddo

c     plot data
      call sk_setforecol(0.0,1.0,0.0) ! green
      do i=1,npts
         call sk_plotsym(x(i),y(i),'X')
      enddo

      call sk_beep
      call sk_waitkey
      call sk_end

      end

c Polynomial
      function calcpoly(n,co,x)
      dimension co(*)

      cum=co(n)
      do i=n-1,1,-1
	 cum=cum*x+co(i)
      enddo
      calcpoly=cum

      return
      end