ccccccccccccccccccccccccc     Program 2.3     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM MILLIKAN
C
C Main program for a linear fit of the Millikan experimental
C data on the fundamental charge e_0 from e_n = e_0*n + de.
C
      PARAMETER (N = 15)
      PARAMETER (M =  2)
      DIMENSION X(N),F(N),A(M),U(M,N)
      DATA X /4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,
     *        12.0,13.0,14.0,15.0,16.0,17.0,18.0/
      DATA F /6.558,8.206,9.880,11.50,13.14,14.81,16.40,18.04,
     *        19.68,21.32,22.96,24.60,26.24,27.88,29.52/
C
      CALL PFIT(N,M,X,F,A,U)
      SUM0 = 0.
      SUMT = 0.
      DO      100  I = 1, N
        SUM0 = SUM0 + U(1,I)**2
        SUMT = SUMT + X(I)*U(1,I)**2
  100 CONTINUE
      E0   = A(2)
      DE   = A(1) - A(2)*SUMT/SUM0
      WRITE (6,999) E0,DE
      STOP
  999 FORMAT (2F16.8)
      END
C
      SUBROUTINE PFIT(N,M,X,F,A,U)
C
C Subroutine generating orthonormal polynomials U(N,M) up to
C (M-1)th order and coefficients A(M), for the least squares
C approximation of the function F(N) at X(N). Other variables
C used: G(K) for g_k, H(K) for h_k, S(K) for <u_k|u_k>.
C
      PARAMETER (NMAX = 101)
      PARAMETER (MMAX = 101)
      DIMENSION G(MMAX),H(MMAX),S(MMAX),X(N),F(N),U(M,N),A(M)
C
      IF(N.GT.NMAX) STOP 'Too many points'
      IF(M.GT.MMAX) STOP 'Order too high'
C
C Set up the zeroth order polynomial u_0
C
      DO      100  I = 1, N
        U(1,I) = 1.
  100 CONTINUE
      DO    200   I = 1, N
        TMP  = U(1,I)*U(1,I)
        S(1) = S(1) + TMP
        G(1) = G(1) + X(I)*TMP
        A(1) = A(1) + U(1,I)*F(I)
  200 CONTINUE
      G(1) = G(1)/S(1)
      H(1) = 0.
      A(1) = A(1)/S(1)
C
C Set up the first order polynomial u_1
C
      DO      300  I = 1, N
        U(2,I) = X(I)*U(1,I) - G(1)*U(1,I)
        S(2)   = S(2) + U(2,I)**2
        G(2)   = G(2) + X(I)*U(2,I)**2
        H(2)   = H(2) + X(I)*U(2,I)*U(1,I)
        A(2)   = A(2) + U(2,I)*F(I)
  300 CONTINUE
      G(2) = G(2)/S(2)
      H(2) = H(2)/S(1)
      A(2) = A(2)/S(2)
C
      IF(M.LT.3) RETURN
C
C Higher order polynomials u_k from the recursive relation
C
      DO      700  I = 2, M-1
        DO    500  J = 1, N
          U(I+1,J) = X(J)*U(I,J) - G(I)*U(I,J) - H(I)*U(I-1,J)
          S(I+1)   = S(I+1) + U(I+1,J)**2
          G(I+1)   = G(I+1) + X(J)*U(I+1,J)**2
          H(I+1)   = H(I+1) + X(J)*U(I+1,J)*U(I,J)
          A(I+1)   = A(I+1) + U(I+1,J)*F(J)
  500   CONTINUE
        G(I+1) = G(I+1)/S(I+1)
        H(I+1) = H(I+1)/S(I)
        A(I+1) = A(I+1)/S(I+1)
  700 CONTINUE
      RETURN
      END