!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: poisson.f90                                         !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: July 7, 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 POISSON
!
! Program for the one-dimensional Poisson equation.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N = 99, M = 2
  INTEGER :: I
  REAL :: PI, H, XM, XI, XP, X, MH
  REAL, DIMENSION (N) :: D, B, C, U
!
  PI  =  4.0*ATAN(1.0)
  H   =  1.0/(N+1)
!
! Evaluate the coefficient matrix elements
!
  DO I = 1, N
    D(I) =  2
    C(I) = -1
    XM = H*(I-1)
    XI = XM + H
    XP = XI + H
    B(I) = 2*SIN(PI*XI) -SIN(PI*XM) - SIN(PI*XP)
  END DO
!
! Obtain the solution
!
  CALL TRIDIAGONAL_LINEAR_EQ (N, D, C, C, B, U)
!
! Output the result
!
  X = H
  MH = M*H
  DO I=1, N, M
    WRITE (6,"(2F16.8)") X, U(I)
    X = X + MH
  END DO
!
END PROGRAM POISSON
!
SUBROUTINE TRIDIAGONAL_LINEAR_EQ (L, D, E, C, B, Z)
!
! Functione to solve the tridiagonal linear equation set.
!
  INTEGER, INTENT (IN) :: L
  INTEGER :: I
  REAL, INTENT (IN), DIMENSION (L):: D, E, C, B
  REAL, INTENT (OUT), DIMENSION (L):: Z
  REAL, DIMENSION (L):: Y, W
  REAL, DIMENSION (L-1):: V, T
!
! Evaluate the elements in the LU decomposition
!
  W(1) = D(1)
  V(1)  = C(1)
  T(1)  = E(1)/W(1)
  DO I = 2, L - 1
    W(I) = D(I)-V(I-1)*T(I-1)
    V(I) = C(I)
    T(I) = E(I)/W(I)
  END DO
  W(L) = D(L)-V(L-1)*T(L-1)
!
! Forward substitution to obtain y
!
  Y(1) = B(1)/W(1)
  DO I = 2, L
    Y(I) = (B(I)-V(I-1)*Y(I-1))/W(I)
  END DO
!
! Backward substitution to obtain z
  Z(L) = Y(L)
  DO I = L-1, 1, -1
    Z(I) = Y(I) - T(I)*Z(I+1)
  END DO
END SUBROUTINE TRIDIAGONAL_LINEAR_EQ