c ----------------------------------------------------------------------
c Fit data with polynomials
c Program Note: Link with LFIT of "Numerical Recipes (2nd Edition)"
c Instructor: Nam Sun Wang
c ----------------------------------------------------------------------

c Declare arrays -------------------------------------------------------
      parameter (nmax=1000, mmax=10)
      dimension x(nmax), y(nmax), sig(nmax)
      dimension lista(mmax), a(mmax), covar(mmax,mmax)
      character fname*40
      external fcn

c Program header -------------------------------------------------------
      print *, 'Demonstrate how to fit data with a polynomial.'

c Read in data ---------------------------------------------------------
      print '(a\)', ' Filename containing data: '
      read(*,'(a)') fname
      open(1,file=fname,status='old')
      do i=1, nmax
19      read(1,*,err=19,end=21) x(i), y(i)
      enddo
21    close(1)
      ipt=i-1
      print *, ' Number of points read = ', ipt

c Assign the standard deviation of the data points to 1. ---------------
      sig = 1.

c Find all of the coefficients of the polynomial -----------------------
      lista = 1

c Find the order of fit ------------------------------------------------
100   print '(a\)', ' Order of fit (quit=0, linear=2, quad=3, etc): '
      read(*,*) m
      if(m .le. 0)stop
      if(m .gt. mmax)m=mmax

c Find the polynomial coefficients by calling LFIT of "Numerical Recipes"
c Usage:
c     call LFIT(x, y, sig, ndat, a, ia, ma, covar, npc, chisq, funcs)
      call LFIT(x, y, sig, ipt, a, lista, m, covar, mmax, chisq, fcn)

c Print the coefficients -----------------------------------------------
      print *, 'The coefficients are (from the lowest order up): '
      do i=1, m
        print '(a, i2, a, 1pe14.6)', ' a(', i, ')= ', a(i)
      enddo

c Print the residual error ---------------------------------------------
      print *, 'squared residual error = ', chisq
      goto 100

      end



c ----------------------------------------------------------------------
      subroutine fcn(x,basefcn,iorder)
c ----------------------------------------------------------------------
c Generate the base functions: basefcn(j), which are power series in x
c ----------------------------------------------------------------------
      dimension basefcn(1)

      do j=1, iorder
        basefcn(j) = x**(j-1)
      enddo

      end