real function FCN_VECFIT (n, x, y, dy, npar, iflag, npfit)

c  Glen Cowan
c  5 December, 1999
c  Sample FCN routine for use with PAW.  
c  User function YFIT (below) is currently a 4th order polynomial.

        implicit           NONE

c  arguments

        integer            n         ! number of data points
        integer            npar      ! number of fit parameters
        integer            iflag     ! tells what to do
        integer            npfit     ! number of data points used

        real               x(*), y(*), dy(*)

c  local variables

        integer            i
        integer            j
        integer            nplot

        real               chi2
        real               f
        real               fplot
        vector             fplot              ! makes visible to PAW
        real               xmin
        real               xmax
        real               xplot
        vector             xplot 
        real               xrange
        real               YFIT               ! fit function (defined below)

        double precision   covmat_dp(5,5)     ! dimension must be set by hand

c  Special for use with PAW

        real               fit_info
        vector             fit_info

        real               par, dpar, covmat
        vector             par, dpar, covmat
        
        double precision   FITPAD, FITFUN
        common  /HCFITD/   FITPAD(24), FITFUN

c  begin

        if (iflag .eq. 1) then
        
c  do any function initialization here

          nplot = 100
          do i = 1, n
            if ( i .eq. 1 ) then
              xmin = x(i)
              xmax = x(i)
            else
              xmin = MIN(xmin, x(i))
              xmax = MAX(xmax, x(i))
            endif
          end do
          xrange = xmax - xmin
          xmin = xmin - 0.1*xrange
          xmax = xmax + 0.1*xrange

        endif

c  compute chi2

        chi2 = 0.
        npfit = 0
        do i = 1, n
          f = YFIT (x(i))
          if ( dy(i) .gt. 0. ) then
            chi2 = chi2 + (y(i) - f)**2 / dy(i)**2
            npfit = npfit + 1
          endif
        end do
        FITFUN = chi2         ! passed back through common block HCFITD

c  Copy fitted parameters and errors into variables accessible to PAW.
c  Also make array of function values for plot.

        if ( iflag .eq. 3 ) then

          call MNEMAT (covmat_dp, npar)
          do i = 1, npar
            par(i) = FITPAD(i)
            dpar(i) = SQRT (covmat_dp(i,i))
            do j = 1, npar
              covmat(i,j) = covmat_dp(i,j)
            end do
          end do

          do i = 1, nplot
            xplot(i) = (xmax - xmin)*FLOAT(i-1)/FLOAT(nplot-1) + xmin
            fplot(i) = YFIT(xplot(i))
          end do

          fit_info(1) = chi2
          fit_info(2) = npfit
          fit_info(3) = npar

        endif

        return
        END

c*************************************************************************

        real function YFIT (x)

c  The fit function

        implicit           NONE

        real               x

c  function parameters

        double precision   FITPAD, FITFUN
        common  /HCFITD/   FITPAD(24), FITFUN

c  local variables

        real               a0, a1, a2, a3, a4

c  begin

        a0 = FITPAD(1)
        a1 = FITPAD(2)
        a2 = FITPAD(3)
        a3 = FITPAD(4)
        a4 = FITPAD(5)

        yfit = a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4

        return
        END