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

c  Author:  Glen Cowan
c  Date:    8-Dec-1998

c  Generic FCN routine for use with MINUIT

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

c  Use KUIP routines to read in values:
c  KUPROR for real, KUPROI for integer, KUPROS for character string

      implicit        NONE

c  constants

      integer           max_points
      parameter         (max_points = 100)  ! same as in driver prog
      integer           nf
      parameter         (nf = 100)          ! number of points for curve
      integer           npar_max
      parameter         (npar_max = 100)

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             / 'ls_data.dat' /
      character*80      line
      character*80      title              / ' ' /

      integer           i
      integer           ios
      integer           istat
      integer           len
      integer           lun

      integer           sym_type

      logical           HEXIST
      logical           store_plot        / .false. /

      real              dx(max_points)    / max_points*0. /
      real              sym_size          / 0.15 /
      real              xf (nf), yf (nf)
      real              x_min, x_max
      real              y_min, y_max

      double precision  f                 ! fit function, defined below
      double precision  xdp               ! double precision version of x
      double precision  xfdp (nf)

c  begin

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

        call KUPROS ('enter file name for input data', infile, len)
        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), dy(i)
          endif
        end do
        close (lun)
        num_points = i

      endif

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

      endif

c  calculate chi2

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

      if ( iflag .eq. 3 ) then        ! comes here to make plots

        call SET_PLOT_OPTIONS         ! sets up sizes of various things

        call KUPROI ('store plot in ps file (0=no,1=yes)', store_plot)
        if ( store_plot ) call PS_OUTPUT

c  Set plot limits (hit carriage return to accept current values)
        
        call KUPROR ('enter x_min', x_min)
        call KUPROR ('enter x_max', x_max)
        call KUPROR ('enter y_min', y_min)
        call KUPROR ('enter y_max', y_max)
        call KUPROS ('enter x axis title', title, len)

c  Set up the plot with HBOOK/HPLOT routines

	if ( HEXIST (1) ) call HDELET (1)
	call HBOOK1 (1, title, 2, x_min, x_max, 0.)
        call HMINIM (1, y_min)
        call HMAXIM (1, y_max)
	call HPLSET ('DMOD', 1.)                    ! solid line
	call HPLOT (1, ' ', 'HIST', 0)

c  plot the data values

        sym_type = 20                               ! circle
        sym_size = 0.15 
        call HPLERR (x, y, dx, dy, num_points, ' ', sym_type, sym_size)

c  plot the fitted function (evaluate at nf points)

        do i = 1, nf
          xf(i) = (x_max - x_min)*FLOAT(i-1)/FLOAT(nf-1) + x_min
          xfdp (i) = xf(i)
          yf(i) = f (xfdp(i), par, npar)
        end do
        call HPLFUN (xf, yf, nf, ' ')

c  show parameter values and chi2 to left of plot

        call SHOW_PAR
        call IGTERM                                 ! force to screen

        if ( store_plot ) call PS_OUTPUT            ! closes ps file

      endif

      return
      END

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      double precision function f (x, par, npar)

c  A typical fit function, here a polynomial of order m=npar-1.  

      implicit NONE

      integer             i, npar

      double precision    par (*)
      double precision    x

c  begin

      f = par(1)              !  zeroth order term
      do i = 2, npar
        if ( x .ne. 0. ) then
          f = f + par(i) * x**(i-1)
        else
          f = f + par(i)
        endif
      end do

      return
      END