program least6
c
c -- least-squares fit to p-v-t of steam
c -- subroutines gausj and square required
c -- figure 7.11
c     
      integer in, out, maxr, maxc, nrow, ncol
      real y(20), ycalc(20), coef(3), correl
      real sig(3), resid(20), p(29), t(20), v(20)
      common /inout/ in, out, maxr, maxc
c
      in = 5
      out = 6
      maxr = 20
      maxc = 3
      call input(p, t, v, nrow)
10    call linfit(p, t, v, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      call output(p, t, v, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      goto 10
      end
      subroutine input(p, t, v, nrow)
c -- get values for nrow and arrays t and p
c
      integer nrow, i
      real t(1), p(1), v(1)
c
      nrow = 12
      t(1) = 400
      p(1) = 120
      v(1) = 4.079
      t(2) = 450
      p(2) = 120
      v(2) = 4.36
      t(3) = 500
      p(3) = 120
      v(3) = 4.633
      t(4) = 400
      p(4) = 140
      v(4) = 3.466
      t(5) = 450
      p(5) = 140
      v(5) = 3.713
      t(6) = 500
      p(6) = 140
      v(6) = 3.952
      t(7) = 400
      p(7) = 160
      v(7) = 3.007
      t(8) = 450
      p(8) = 160
      v(8) = 3.228
      t(9) = 500
      p(9) = 160
      v(9) = 3.44
      t(10) = 400
      p(10) = 180
      v(10) = 2.648
      t(11) = 450
      p(11) = 180
      v(11) = 2.85
      t(12) = 500
      p(12) = 180
      v(12) = 3.042
c -- convert t to rankine
      do 20 i = 1, nrow
       t(i) = t(i) + 460
20    continue
      return
      end
      subroutine output(p, t, v, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
c
c -- print out the answers
c
      integer in, out, nrow, ncol, i
      real p(1), t(1), v(1), y(1), ycalc(1), coef(1), correl
      real resid(1), sig(1), res(20)
      common /inout/ in, out, maxr, maxc
c
      write(out, 101)
      do 10 i = 1, nrow
       res(i) = 100 * resid(i) / y(i)
10    continue
      write(out, 102) (i, p(i), t(i), v(i), y(i),
     *  ycalc(i), res(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    p      t      v        y',
     *  '     y calc    % res')
102   format(i4, 2f7.0, f7.3, 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(p, t, v, y, ycalc, resid, coef, sig,
     *  nrow, ncol, cor)
c
c -- least squares fit to nrow sets of x-y points
c -- for properties of steam
c -- using gauss-jordan elimination
c -- subroutines square and gaussj needed
c
      logical error
      integer nrow, ncol, i, j, maxr, maxc, in, out
      integer index(20,3), nvec
      real p(1), t(1), v(1), y(1), ycalc(1), coef(1), resid(1)
      real a(3,3), xmatr(20,3), sig(1), gas, power
      real sumy, sumy2, yi, yc, res, cor, srs, see
      common /inout/ in, out, maxr, maxc
      data nvec/1/, gas/85.76/
c
      ncol = 2
      write(out, 101)
      read(in, 102) power
      if (power .lt. -10.0) goto 99
      do 10 i = 1, nrow
       xmatr(i,1) = p(i) / t(i)**power
       xmatr(i,2) = sqrt(p(i))
       y(i) = v(i) * p(i) - gas * t(i) / 144
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
99    stop
101   format(/' power? ' )
102   format(e10.0)
      end