!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: millikan.f90                                        !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: June 14, 2006                                           !
!                                                                         !
!  (1) This F90 program is created for the book, "An Introduction to      !
!      Computational Physics, 2nd Edition," written by Tao Pang and       !
!      published by Cambridge University Press on January 19, 2006.       !
!                                                                         !
!  (2) No warranties, express or implied, are made for this program.      !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM MILLIKAN
!
! An example of applying the least-squares approximation
! to the Millikan data on a straight line q=k*e+dq.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=14, M=21
  INTEGER :: I
  REAL :: SUM, E, DQ
  REAL, DIMENSION (N+1) :: K, F
  REAL, DIMENSION (M+1,N+2) :: U
  DATA K /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/
!
  CALL ORTHPOLY(N, M, K, F, U)
  SUM= 0
  DO I = 1, N+1
    SUM = SUM + K(I)
  END DO
  E = U(2,N+2)
  DQ   = U(1,N+2)-U(2,N+2)*SUM/(N+1)
  WRITE (6,"(2F16.8)") E, DQ
END PROGRAM MILLIKAN
!
SUBROUTINE ORTHPOLY(N, M, K, F, U)
!
! Subprogram to generate the orthogonal polynomials and the
! least-squares fitting coefficients.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: NMAX=101, MMAX=101
  INTEGER, INTENT (IN) :: N
  INTEGER, INTENT (INOUT) :: M
  INTEGER :: I,J
  REAL :: STMP
  REAL, INTENT (IN), DIMENSION (N+1) :: K, F
  REAL, INTENT (OUT), DIMENSION (M+1,N+2) :: U
  REAL, DIMENSION (N+1) :: S, G, H
!
  IF (M .GT. N) THEN
    M = N
    WRITE (6, "(A, I3)") "The highest power is adjusted to ", N
  END IF
!
! Set up the zeroth-order polynomial
!
  S(1) = 0
  G(1) = 0
  U(1,N+2) = 0
  H(1) = 0
  DO I = 1, N+1
    U(1,I) = 1
  END DO
  DO I = 1, N+1
    STMP = U(1,I)*U(1,I)
    S(1) = S(1) + STMP
    G(1) = G(1) + K(I)*STMP
    U(1,N+2) = U(1,N+2) + U(1,I)*F(I)
  END DO
  G(1) = G(1)/S(1)
  U(1,N+2) = U(1,N+2)/S(1)
!
! Set up the first-order polynomial
!
  S(2) = 0
  G(2) = 0
  U(2,N+2) = 0
  H(2) = 0
  DO I = 1, N+1
    U(2,I) = K(I)*U(1,I) - G(1)*U(1,I)
    S(2) = S(2) + U(2,I)**2
    G(2) = G(2) + K(I)*U(2,I)**2
    U(2,N+2) = U(2,N+2) + U(2,I)*F(I)
    H(2) = H(2) + K(I)*U(2,I)*U(1,I)
  END DO
  G(2) = G(2)/S(2)
  U(2,N+2) = U(2,N+2)/S(2)
  H(2) = H(2)/S(1)
!
! Higher order polynomials u_k from the recursive relation
!
  IF(M.GE.2) THEN
    DO I = 2, M
      S(I+1) = 0
      G(I+1) = 0
      U(I+1,N+2) = 0
      DO J = 1, N+1
        U(I+1,J) = K(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) + K(J)*U(I+1,J)**2
        U(I+1, N+2) = U(I+1,N+2) + U(I+1,J)*F(J)
        H(I+1) = H(I+1) + K(J)*U(I+1,J)*U(I,J)
      END DO
      G(I+1) = G(I+1)/S(I+1)
      U(I+1,N+2) = U(I+1,N+2)/S(I+1)
      H(I+1) = H(I+1)/S(I)
    END DO
  END IF
END SUBROUTINE ORTHPOLY