subroutine FCN (npar, grad, chi2, par, iflag, futil)

c  Author:  Glen Cowan
c  Date:    13-OCT-1997

c  Program to fit a function d = alpha*(h**beta) to measurements of
c  impact distance d of a ball rolling down a ramp of height h.

c  Input:  integer            npar         number of parameters to fit
c          double precision   par(npar)    parameter vector
c          integer            iflag        select what to do
c          double precision   futil        optional external function
c
c  Output: double precision   grad(npar)   gradient vector (not filled)
c          double precision   chi2         chi^2 to be minimized

      implicit        NONE

c  constants

      integer           max_points
      parameter         (max_points = 100)  ! same as in driver prog

c  arguments and common blocks

      integer           npar,iflag
      double precision  futil, chi2, par(*), grad(*)

      integer           num_points
      real              x(max_points), y(max_points), dy(max_points)
      common            / data_block / num_points, x, y, dy

c  local variables

      character*80      infile
      character*80      line

      integer           i
      integer           ios
      integer           lun

      double precision  alpha
      double precision  beta
      double precision  f
      double precision  sigma

c  begin

      if ( iflag .eq. 1 ) then        ! read input data

        write (*, *) 'enter file name for input data'
        read (*, '(a80)') infile
        lun = 30
        open (unit=lun, file = infile,
     1  form = 'formatted', status = 'old')

        i = 0
        ios = 0
        do while ( ios .eq. 0 )
          read (lun, fmt = '(a80)', iostat = ios) line
          if (ios .eq. 0 .and. .not. line(1:1) .eq. '!') then
            i = i + 1
            read (line, *) x(i), y(i)
          endif
        end do
        close (lun)
        num_points = i

c  Set measurement error (here take equal for all points)

        write (*, *) 'enter standard error of measurements'
        read (*, *) sigma
        do i = 1, num_points
          dy(i) = sigma
        end do

      endif

      if ( iflag .eq. 2 ) then        ! calculate gradient vector

      endif

c  calculate chi2

      alpha = par(1)
      beta  = par(2)

      chi2 = 0.
      do i = 1, num_points
        f = alpha * (x(i)**beta)
        chi2 = chi2 + (y(i) - f)**2 / dy(i)**2
      end do

      if ( iflag .eq. 3 ) then        ! finish up

      endif

      return
      END