program least4
c
c -- least-squares fit to the heat capacity equation
c -- subroutines gausj, square, and plot required
c -- figure 7.7
c     
      integer out, maxr, maxc, lines, nrow, ncol
      real x(20), y(20), ycalc(20), coef(3), correl
      real sig(3), resid(20)
      character*1 ch
      common /inout/ out, maxr, maxc
c
      out = 6
      maxr = 20
      maxc = 3
      call input(x, y, nrow)
      call linfit(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      call output(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      write(*,*) 'Press RETURN for plot'
      read(*,'(a1)') ch
      lines = 2*(nrow - 1) + 1
      call plot(x, y, ycalc, nrow, out, lines)
      stop
      end
      subroutine input(t, cp, nrow)
c -- get values for nrow and arrays t and cp
c
      integer nrow, i
      real t(1), cp(1)
c
      nrow = 10
      do 10 i = 1, nrow
       t(i) = (i + 2) * 100.0
10    continue
      cp(1) =  7.02
      cp(2) =  7.2
      cp(3) =  7.43
      cp(4) =  7.67
      cp(5) =  7.88
      cp(6) =  8.06
      cp(7) =  8.21
      cp(8) =  8.34
      cp(9) =  8.44
      cp(10) = 8.53
      return
      end
      subroutine output(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
c
c -- print out the answers
c
      integer out, nrow, ncol, i
      real x(1), y(1), ycalc(1), coef(1), correl
      real resid(1), sig(1)
      common /inout/ out, maxr, maxc
c
      write(out, 101)
      write(out, 102) (i, x(i), y(i), ycalc(i),
     *  resid(i), i = 1, nrow)
      write(out, 103)
      write(out, 106) coef(1), sig(1)
      write(out, 104) (coef(i), sig(i), i = 2, ncol)
      write(out, 105) correl
      return
101   format('   i      x       y     y calc    resid')
102   format(i4, f8.1, 3f9.2)
103   format(/' coefficients    errors')
104   format(1x, 1pe12.3, 3x, e12.3)
105   format(/' correlation coefficient is', f7.4)
106   format(1x, 1pe12.3, 3x, e12.3, '  constant term')
      end
      subroutine linfit(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, cor)
c
c -- least squares fit nrow sets of x-y points
c -- for heat capacity equation
c -- using gauss-jordan elimination
c -- subroutines square and gaussj needed
c
      logical error
      integer nrow, ncol, i, j, maxr, maxc, out
      integer index(20,3), nvec
      real x(1), y(1), ycalc(1), coef(1), resid(1)
      real a(3,3), xmatr(20,3), sig(1)
      real sumy, sumy2, xi, yi, yc, res, cor, srs, see
      common /inout/ out, maxr, maxc
      data nvec/1/
c
      ncol = 3
      do 10 i = 1, nrow
       xi = x(i)
       xmatr(i,1) = 1.0
       xmatr(i,2) = xi
       xmatr(i,3) = 1.0 / (xi * xi)
10    continue
      call square(xmatr, y, a, coef, nrow, ncol, maxr, maxc)
      call gaussj(a, coef, index, ncol, maxc, nvec, error, out)
      sumy = 0.0
      sumy2 = 0.0
      srs = 0.0
      do 20 i = 1, nrow
       yi = y(i)
       yc = 0.0
       do 15 j = 1, ncol
15      yc = yc + coef(j) * xmatr(i,j)
       ycalc(i) = yc
       res = yc - yi
       resid(i) = res
       srs = srs + res*res
       sumy = sumy + yi
       sumy2 = sumy2 + yi * yi
20    continue
      cor = sqrt(1.0 - srs/(sumy2 - sumy*sumy/nrow))
      if (nrow .eq.ncol) see = sqrt(srs) 
      if (nrow .ne.ncol) see = sqrt(srs / (nrow - ncol))
      do 30 i = 1, ncol
30     sig(i) = see * sqrt(a(i,i))
      return
      end