program cfit4
c
c -- linear least-squares fit
c -- figure 5.10
c     
      integer length, out, max, lines
      real x(20), y(20), a, b, ycalc(20)
      character*1 ch
      common /inout/ out, max
c
      out = 6
      max = 20
10    call input(x, y, length)
      call linfit(x, y, ycalc, a, b, length)
      call output(x, y, ycalc, length, a, b)
      write(*,*) 'Press RETURN for plot'
      read(*,'(a1)') ch
      lines = 2*(length - 1) + 1
      call plot(x, y, ycalc, length, out, lines)
      goto 10
      end
      subroutine input(x, y, n)
c -- get values for n and arrays x and y
c
      integer n, out, i, j, max
      real x(n), y(n), fudge
      common /inout/ out, max
      data a/2.0/, b/5.0/
c
      write(out, 101)
      read(*,103) fudge
      if (fudge .lt. 0.0) goto 99
5     write(out, 102)
      read(*,104) n
      if (n .gt. max) goto 5
      if (n .lt. 0) goto 99
      do 10 i = 1, n
       j = n + 1 - i
       x(i) = j
       y(i) = (a + b*j)*(1.0+(2.0*rand(0) - 1.0)*fudge)
10    continue
      return
99    stop
101   format(/' fudge? ' )
102   format(' how many points? ' )
103   format(f10.0)
104   format(i2)
      end
      subroutine output(x, y, ycalc, n, a, b)
c
c -- print out the answers
c
      integer n, out, i
      real x(1), y(1), ycalc(1)
      common /inout/ out, max
c
      write(out, 101)
      write(out, 102) (i, x(i), y(i), ycalc(i), i = 1, n)
      write(out, 103) a
      write(out, 104) b
      return
101   format('   i      x       y     y calc')
102   format(i4, f8.1, 2f9.2)
103   format(/' intercept is', f7.2)
104   format('     slope is', f7.2)
      end
      subroutine linfit(x, y, ycalc, a, b, n)
c
c -- fit a straight line for (ycalc) through
c -- n pairs of x-y points
c
      integer n, i
      real x(1), y(1), ycalc(1), a, b, sxy, syy, sxx
      real sumx, sumy, sumxy, sumx2, sumy2, xi, yi
c
      sumx = 0.0
      sumy = 0.0
      sumxy = 0.0
      sumx2 = 0.0
      sumy2 = 0.0
      do 10 i = 1, n
       xi = x(i)
       yi = y(i)
       sumx = sumx + xi
       sumy = sumy + yi
       sumxy = sumxy + xi * yi
       sumx2 = sumx2 + xi * xi
       sumy2 = sumy2 + yi * yi
10    continue
      sxx = sumx2 - sumx * sumx / n
      sxy = sumxy - sumx * sumy / n
      syy = sumy2 - sumy * sumy / n
      b = sxy / sxx
      a = (sumy - b * sumx) / n
      do 20 i = 1, n
       ycalc(i) = a + b * x(i)
20    continue
      return
      end