program FIT_GALILEO

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  Fitting done with the method of least squares.  The chi^2 minimization
c  is carried out with the MINUIT routines from the CERN program library.
c  Program must be linked with MINUIT and with FCN_GALILEO.

      implicit        NONE

c  local variables

      integer            n_dof

      double precision   chi2

c  variables for MINUIT

      integer            npar
      parameter          (npar = 2)

      character*10       chnam(npar)

      integer            ierr
      integer            ird
      integer            isav
      integer            istat
      integer            ivarbl
      integer            iwr
      integer            npari, nparx

      double precision   arglis (10)
      double precision   bnd1, bnd2
      double precision   deriv (npar)
      double precision   dpar (npar)
      double precision   fmin, fedm, errdef
      double precision   covmat (npar, npar)
      external           FCN
      double precision   par (npar)

c  common block for data 

      integer            max_points
      parameter          (max_points = 100)  ! must be same as in FCN
      integer            num_points
      real               x(max_points), y(max_points), dy(max_points)
      common             / data_block / num_points, x, y, dy

c  Initialize MINUIT

      ird = 5         ! unit number for input to Minuit (keyboard)
      iwr = 6         ! unit number for output from Minuit (screen)
      isav = 7        ! unit number for use of the SAVE command
      call MNINIT (ird, iwr, isav)

c  Set print level to -1 (no output except SHOW commands)

      arglis(1) = -1.d0
      call MNEXCM (FCN, 'SET PRIN', arglis, 1, ierr, 0)

c  Define parameters alpha and beta, give initial values and step sizes.

      call MNPARM (1, 'alpha', 50.d0, 1.d0, 0.d0, 0.d0, ierr)
      call MNPARM (2, 'beta', 0.5d0, 0.1d0, 0.d0, 0.d0, ierr)

c  Get input data by calling FCN with iflag=1

      arglis(1) = 1.d0
      call MNEXCM (FCN, 'CALL', arglis, 1, ierr, 0)

c  Minimize chi2 using SIMPLEX and MIGRAD 

      call MNEXCM (FCN, 'SIMPLEX', arglis, 0, ierr, 0)
      call MNEXCM (FCN, 'MIGRAD', arglis, 0, ierr, 0)

c  Compute covariance matrix with HESSE

      call MNEXCM (FCN, 'HESSE', arglis, 0, ierr, 0)

c  Get results of fit

      call MNSTAT (fmin, fedm, errdef, npari, nparx, istat)
      chi2 = fmin
      n_dof = num_points - npari
      call MNPOUT (1, chnam(1), par(1), dpar(1), bnd1, bnd2, ivarbl)
      call MNPOUT (2, chnam(2), par(2), dpar(2), bnd1, bnd2, ivarbl)
      call MNEMAT (covmat, npar)

      write (*, *)
      write (*, *) 'Fit results:'
      write (*, *)
      write (*, *) 'alpha           = ', par(1), ' +- ',  dpar(1)
      write (*, *) 'beta            = ', par(2), ' +- ',  dpar(2)
      write (*, *) 'cov[alpha,beta] = ', covmat(1,2)
      write (*, *) 'rho[alpha,beta] = ', covmat(1,2)/(dpar(1)*dpar(2))
      write (*, *) 'chi2            = ', chi2
      write (*, *) 'n_dof           = ', n_dof
      write (*, *) 'chi2/n_dof      = ', chi2/FLOAT(n_dof)

      stop
      END