program least1
c
c -- parabolic least-squares fit
c -- function determ and
c -- subroutines cramer, setup, and plot required
c -- figure 7.1
c     
      integer out, maxr, maxc, lines, nrow, ncol
      real x(20), y(20), ycalc(20), coef(3), correl
      character*1 ch
      common /inout/ out
c
      out = 6
      maxr = 20
      maxc = 3
      call input(x, y, nrow)
      call linfit(x, y, ycalc, coef, nrow, ncol, correl)
      call output(x, y, ycalc, coef, nrow, ncol, correl)
      lines = 2*(nrow - 1) + 1
      write(out,*) 'Press RETURN for plot'
      read(*,'(a1)') ch
      call plot(x, y, ycalc, nrow, out, lines)
      stop
      end
      subroutine input(x, y, nrow)
c -- get values for nrow and arrays x and y
c
      integer nrow, i
      real x(1), y(1)
c
      nrow = 9
      do 10 i = 1, nrow
       x(i) = i
10    continue
      y(1) =  2.07
      y(2) =  8.6
      y(3) = 14.42
      y(4) = 15.8
      y(5) = 18.92
      y(6) = 17.96
      y(7) = 12.98
      y(8) =  6.45
      y(9) =  0.27
      return
      end
      subroutine output(x, y, ycalc, coef, 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
      common /inout/ out
c
      write(out, 101)
      write(out, 102) (i, x(i), y(i), ycalc(i), i = 1, nrow)
      write(out, 103)
      write(out, 104) (coef(i), i = 1, ncol)
      write(out, 105) correl
      return
101   format('   i      x       y     y calc')
102   format(i4, f8.1, 2f9.2)
103   format(/' coefficients')
104   format(f9.4)
105   format(/' correlation coefficient is', f7.4)
      end
      subroutine linfit(x, y, ycalc, coef, nrow, ncol, cor)
c
c -- least squares fit to a parabola
c -- nrow pairs of x-y points
c
      logical error
      integer nrow, ncol, i
      real x(1), y(1), ycalc(1), coef(1)
      real a(3,3), g(3)
      real sumx, sumy, sumxy, sumx2, sumy2, xi, yi
      real sumx3, sumx4, srs, x2, res, cor
c
      ncol = 3
      sumx = 0.0
      sumy = 0.0
      sumxy = 0.0
      sumx2 = 0.0
      sumy2 = 0.0
      sumx3 = 0.0
      sumx4 = 0.0
      sum2y = 0.0
      do 10 i = 1, nrow
       xi = x(i)
       yi = y(i)
       x2 = xi * xi
       sumx = sumx + xi
       sumy = sumy + yi
       sumxy = sumxy + xi * yi
       sumx2 = sumx2 + x2
       sumy2 = sumy2 + yi * yi
       sumx3 = sumx3 + xi * x2
       sumx4 = sumx4 + x2 * x2
       sum2y = sum2y + x2 * yi
10    continue
      a(1,1) = nrow
      a(2,1) = sumx
      a(1,2) = sumx
      a(3,1) = sumx2
      a(1,3) = sumx2
      a(2,2) = sumx2
      a(3,2) = sumx3
      a(2,3) = sumx3
      a(3,3) = sumx4
      g(1) = sumy
      g(2) = sumxy
      g(3) = sum2y
      call cramer(a, g, coef, error)
      srs = 0.0
      do 20 i = 1, nrow
       ycalc(i) = coef(1) + coef(2)*x(i) + coef(3)*x(i)*x(i)
       res = y(i) - ycalc(i)
       srs = srs + res*res
20    continue
      cor = sqrt(1.0 - srs/(sumy2 - sumy*sumy/nrow))
      return
      end