!=======================================================================
!
!      Course       : CPS615
!
!      Program Name : laplace.f
!
!      Author       : Kim, Dongmin (096-08-3342)
!
!      In this program, we solve the Laplace Heat equation using
!
!      Jacobi iteration algorithm.
!
!      The ranges are 0 <= x <= 2 and 0 <= y <= 1 and
!
!      the theoretical solution is
!
!      U(x, y) = -4.7 + 0.2x - 0.3y + 0.09(x**2 - y**2) +
!                ln(sqrt(((x + 0.2)**2) + ((y + 0.2)**2)))
!
!      The Jacobi iteration is executed two time,
!
!      one is with 0 initial value, and the other is with the average
!
!      of the all boundary values.
!
!       In this program, 1 dimension is used as a column and
!
!       2 dimension is used as a row e.g.
!
!      A(1,1)    A(1,2)   . . .  A(1,Max)
!      A(2,1)    A(2,2)   . . .  A(2,Max)
!        .         .               .
!        .         .               .
!        .         .               .
!
!      A(Max,1)  A(Max,2) . . .  A(Max,Max)
!
!=======================================================================

program LaplaceHeat

      implicit none
      integer, parameter :: Max = 40            ! The maximun size of the problem
      double precision, parameter :: Xmin = 0.0, Xmax = 2.0 ! The range of X
      double precision, parameter :: Ymin = 0.0, Ymax = 1.0 ! The range of Y
      double precision, dimension( 0 : Max, 0 : Max ) :: &
                          PrevTable, CurrTable  ! The tables
      double precision :: scale, U
      double precision :: x, y, temp, average
      integer :: indxX, indxY
      integer :: i, j

      ! read the scale from consol
      print *, '/ Scale > '
      read *, scale

      PrevTable = 0.0

      indxX = INT( ( Xmax - Xmin ) / scale )
      indxY = INT( ( Ymax - Ymin ) / scale )

      ! setup the boundary value
      average = 0
      do i = 0, indxX
         PrevTable( i, 0 )     = U( i * scale, Ymin )
         PrevTable( i, indxY ) = U( i * scale, Ymax )
         average = average + PrevTable( i, 0 ) + PrevTable( i, indxY )
      end do

      do j = 1, indxY - 1
         PrevTable( 0, j )     = U( Xmin, j * scale )
         PrevTable( indxX, j ) = U( Xmax, j * scale )
         average = average + PrevTable( 0, j ) + PrevTable( indxX, j )
      end do

      ! find the average of boundary values
      average = average / ( 2 * indxX + 2 * indxY )

      CurrTable = PrevTable

      call Jacobi( PrevTable(0:indxX, 0:indxY), CurrTable(0:indxX, 0:indxY), indxX, indxY )

!      do i = 0, indxX
!         print *, CurrTable( i, j )
!      end do

      temp = 0
      do i = 1, indxX - 1
         do j = 1, indxY - 1
            temp = temp + ABS( CurrTable( i, j ) - U( i * scale, j * scale ) )
         end do
      end do

      print *, 'cumulative error : ', temp


      PrevTable( 1 : indxX - 1, 1 : indxY - 1 ) = average

      CurrTable = PrevTable

      call Jacobi( PrevTable(0:indxX, 0:indxY), CurrTable(0:indxX, 0:indxY), indxX, indxY )

!      do i = 0, indxX
!         print *, CurrTable( i, j )
!      end do

      temp = 0
      do i = 1, indxX - 1
         do j = 1, indxY - 1
            temp = temp + ABS( CurrTable( i, j ) - U( i * scale, j * scale ) )
         end do
      end do

      print *, 'cumulative error : ', temp

end program LaplaceHeat


function U( x, y ) result ( boundary )

      implicit none
      double precision :: x, y, boundary

      boundary = -4.7 + 0.2 * x - 0.3 * y + 0.09 * (x**2 - y**2) + &
                 LOG( SQRT( ( (x + 0.2)**2 ) + ( (y + 0.2)**2 ) ) )

end function U


subroutine Jacobi( PrevTable, CurrTable, indxX, indxY )

      implicit none

      double precision, parameter :: CONVERGE = 0.000003
      integer :: indxX, indxY
      double precision, dimension( 0 : indxX, 0 : indxY ) :: &
                          PrevTable, CurrTable  ! The tables
      double precision :: x, y, temp
      real :: test, test2
      integer :: i, j, iteration

      iteration = 0

      do ! 
         test = 0.0
         do i = 1, indxX - 1
            do j = 1, indxY - 1
               temp = 0.25 * &
                  ( PrevTable( i, j + 1 ) + PrevTable( i, j - 1 ) + &
                    PrevTable( i + 1, j ) + PrevTable( i - 1, j ) )
               test2 = ABS( temp - PrevTable( i, j ) )
               test = AMAX1( test, test2 )
               CurrTable( i, j ) = temp
            end do
         end do

         iteration = iteration + 1

         PrevTable = CurrTable

         if ( test .LE. CONVERGE ) then
            print *, 'iteration : ', iteration
            EXIT
         end if

      end do

end subroutine Jacobi
