program lexp
c
c -- least-squares fit for the diffusion of zn in cu
c -- subroutines gausj, square, and plot required
c -- figure 10.3
c     
      integer out, maxr, maxc, lines, nrow, ncol
      real x(20), y(20), ycalc(20), coef(4)
      real resid(20), t(20), d(20)
      character*1 ch
      common /inout/ out, maxr, maxc
c
      out = 6
      maxr = 20
      maxc = 4
      call input(x, y, t, d, nrow)
      call linfit(x, y, ycalc, resid, coef, srs,
     *  nrow, ncol)
      call output(t, d, ycalc, coef, srs, nrow, ncol)
      write(out,*) 'Press RETURN for plot'
      read(*,'(a1)') ch
      lines = 2*(nrow - 1) + 1
      call plot(x, y, y, -nrow, out, lines)
      stop
      end
      subroutine input(x, y, t, d, nrow)
c -- get values for nrow and arrays x and y
c
      integer nrow, i
      real x(1), y(1), t(1), d(1)
c
      nrow = 7
      d(1) =  1.4e-12
      d(2) =  5.5e-12
      d(3) =  1.8e-11
      d(4) =  6.1e-11
      d(5) =  1.6e-10
      d(6) =  4.4e-10
      d(7) =  1.2e-9
      do 10 i = 1, nrow
       t(i) =  550 + 50 * i
       x(i) = 1.0 / (t(i) + 273.15)
       y(i) = alog(d(i))
10    continue
      return
      end
      subroutine output(x, y, ycalc, coef, srs, nrow, ncol)
c
c -- print out the answers
c
      integer out, nrow, ncol, i
      real x(1), y(1), ycalc(1), coef(1), d0, q
      common /inout/ out, maxr, maxc
      data r/ 1.987/
c
      write(out, 101)
      write(out, 102) (i, x(i), y(i), ycalc(i), i = 1, nrow)
      write(out, 103)
      write(out, 106) coef(1)
      write(out, 104) (coef(i), i = 2, ncol)
      d0 = exp(coef(1))
      q = - r * coef(2) / 1000
      write(out, 107) d0, q, srs
      return
101   format('   i    t c       d         d calc')
102   format(i4, 0pf8.1, 1p2e12.3)
103   format(/' coefficients')
104   format(1pe12.3)
106   format(1pe12.3, '  constant term')
107   format(/' d0 =', f9.2, ' cm sq/sec.' /
     *  ' q  =', f9.2, ' kcal/mole' //'  srs = ', f6.2)
      end
      subroutine linfit(x, y, ycalc, resid, coef, srs,
     *  nrow, ncol)
c
c -- least squares fit nrow sets of x-y points
c -- using gauss-jordan elimination
c -- subroutines square and gaussj needed
c
      logical error
      integer nrow, ncol, i, maxr, maxc, out
      integer index(20,3), nvec
      real x(1), y(1), ycalc(1), coef(1), resid(1)
      real a(4,4), xmatr(20,4), a1
      real sumy, sumy2, res, srs
      common /inout/ out, maxr, maxc
      data nvec/1/
c
      ncol = 2
      do 10 i = 1, nrow
       xmatr(i,1) = 1.0
       xmatr(i,2) = x(i)
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
      a1 = exp(coef(1))
      do 20 i = 1, nrow
       ycalc(i) = a1 * exp(coef(2) * x(i))
       if (y(i) .eq. 0.0) res = 1.0
       if (y(i) .ne. 0.0) res = ycalc(i)/y(i) - 1.0
       resid(i) = res
       srs = srs + res*res
20    continue
      return
      end