program nlexp
c
c -- least-squares fit for the diffusion of zn in cu
c -- subroutine newton is required
c -- figure 10.5
c     
      integer out, maxr, maxc, nrow, ncol
      real ycalc(20), coef(4)
      real resid(20), t(20), d(20)
      common /inout/ out, maxr, maxc
      common /fun/ nrow, a, x(20), y(20), ex(20)
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)
      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) = 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)
      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)
      a = coef(1)
      b = - r * coef(2) / 1000
      write(out, 107) a, b, 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 = ', f7.3)
      end
      subroutine linfit(x, y, ycalc, resid, coef, srs,
     *  nrow, ncol)
c
c -- least squares fit nrow sets of x-y points
c -- subroutine newton needed
c
      logical error
      integer nrow, ncol, i, maxr, maxc, out
      real x(1), y(1), ycalc(1), coef(1), resid(1)
      real sumy, xi, yi, res, srs
      external func
      common /inout/ out, maxr, maxc
      common /fun/ idum, a, dummy(40), ex(20)
c
      ncol = 2
      sumx  = 0.0
      sumy  = 0.0
      sumxy = 0.0
      sumx2 = 0.0
      do 10 i = 1, nrow
       xi = x(i)
       yi = alog(y(i))
       sumx  = sumx  + xi
       sumy  = sumy  + yi
       sumxy = sumxy + xi*yi
       sumx2 = sumx2 + xi*xi
10    continue
      b = (sumxy - sumx*sumy/nrow) / (sumx2 - sumx*sumx/nrow)
      call newton(b, func, error)
      coef(1) = a
      coef(2) = b
      srs = 0.0
      do 20 i = 1, nrow
       ycalc(i) = a * ex(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
      subroutine func(b, fb, dfb)
c
c -- calculate function and slope for newton
c
      common /fun/ nrow, a, x(20), y(20), ex(20)
c
      s1 = 0
      s2 = 0
      s3 = 0
      s4 = 0
      s5 = 0
      s6 = 0
      do 10 i = 1, nrow
       xi = x(i)
       x2 = xi*xi
       yi = y(i)
       y2 = yi*yi
       ex1 = exp(b * xi)
       ex(i) = ex1
       ex2 = ex1 * ex1
       s1 = s1 + xi * ex2/y2
       s2 = s2 + ex1/yi
       s3 = s3 + xi * ex1/yi
       s4 = s4 + ex2/y2
       s5 = s5 + 2 * x2 * ex2/y2
       s6 = s6 + x2 * ex1/yi
10    continue
      fb = s1*s2 - s3*s4
      dfb = s2*s5 - s1*s3 - s4*s6
      a = s2/s4
      return
      end