      PROGRAM laplace
*     長方形領域内のラプラス方程式をヤコビの緩和法で解く
*    
      IMPLICIT NONE
      REAL*8 V(0:100,0:100), min_change
      INTEGER iterations, nx, ny
      CALL assign(V,nx,ny,min_change)
      CALL iterate(V,nx,ny,min_change,iterations)
      END
C
      SUBROUTINE assign(V,nx,ny,min_change)
*     長方形領域の1辺の長さ
      IMPLICIT NONE
      REAL*8 V0, V(0:100,0:100), min_change
      INTEGER x, y, nx, ny
      nx = 9                                ! 領域内の x 方向の点の数
      ny = 9                                ! 領域内の y 方向の点の数
      V0 = 10                               ! 長方形の境界のポテンシャル
      WRITE(*,'(A,$)') 'percentage change = '
      READ(*,*) min_change
      min_change = min_change/100
*     長方形の境界の電位を指定
      DO x = 0, nx + 1
         V(x,0) = V0
         V(x,ny+1) = V0
      END DO
      DO y = 0, ny + 1
         V(0,y) = V0
         V(nx+1,y) = V0
      END DO
*     内部の点の電位の値を予想
      DO y = 1, ny
         DO x = 1, nx
            V(x,y) = 0.9D0*V0
         END DO
      END DO
      CALL show_output(V,nx,ny,0)
      END
C
      SUBROUTINE iterate(V,nx,ny,min_change,iterations)
      IMPLICIT NONE
      REAL*8 Vave(100,100), change, dV, V(0:100,0:100), min_change
      INTEGER x, y, nx, ny, iterations
      LOGICAL Ltmp1_
      iterations = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
*        反復計算における格子点の電位の値の最大差
         change = 0
         DO y = 1, ny
            DO x = 1, nx
*              隣接したセルの電位の値を平均
               Vave(x,y) = V(x+1,y) + V(x-1,y)
               Vave(x,y) = Vave(x,y) + V(x,y+1) + V(x,y-1)
               Vave(x,y) = 0.25D0*Vave(x,y)
*              電位の値の変化の割合を計算
               IF(Vave(x,y) .NE. 0) THEN
                  dV = abs((V(x,y) - Vave(x,y))/Vave(x,y))
                  IF(dV .GT. change) change = dV
               END IF
            END DO
         END DO
         DO y = 1, ny
            DO x = 1, nx
               V(x,y) = Vave(x,y)
            END DO
         END DO
         iterations = iterations + 1
         CALL show_output(V,nx,ny,iterations)
         Ltmp1_ = change .GT. min_change
      END DO
      END
C
      SUBROUTINE show_output(V,nx,ny,iterations)
      IMPLICIT NONE
      REAL*8 V(0:100,0:100)
      INTEGER x, y, nx, ny, iterations
      WRITE(*,*)
      WRITE(*,*) 'iteration =', iterations
      DO y = 0, ny + 1
         DO x = nx + 1, 0, - 1
            WRITE(*, '(F7.2,$)') V(x,y)
         END DO
         WRITE(*,*)
      END DO
      END
