program Laplace_Equation
     
       implicit none

       integer, parameter ::  m1 = 20, n1 = 10
       integer, parameter ::  m2 = 40, n2 = 20

       double precision, dimension (0:m1, 0:n1):: T1, error_matrix1
       double precision, dimension (0:m2, 0:n2):: T2, error_matrix2

       integer i, j
       integer iter
       double precision ERROR_SUM, boundary_sum, avg_boundary
       double precision scale, xmin, xmax, ymin, ymax
       double precision denom, init
       xmin = 0.0
       xmax = 2.0
       ymin = 0.0
       ymax = 1.0

       init = 0.0
       scale = 0.1

       print *
       print *, "================================"
       print *, "Problem 2.1a)"
       print *, "================================"

       do i=0, m1
          do j=0, n1
            T1(i, j) = init
          end do
       end do

       do i=0, m1
          call uu(i*scale, ymin, T1(i, 0))
          call uu(i*scale, ymax, T1(i, n1))
       end do

       do j=1, n1
          call uu(xmin, j*scale, T1(0, j))
          call uu(xmax, j*scale, T1(m1, j))
       end do

       do i=0, m1
          do j=0, n1
             call uu(i*scale, j*scale, error_matrix1(i, j))
          end do
       end do       

       call Laplace_Solver(T1, error_matrix1, m1, n1, iter, ERROR_SUM)
       print *
       print *, "Number of iterations:" , iter
       print *, "Culmulative error:   ", ERROR_SUM
       print *

       print *
       print *, "================================"
       print *, "Problem 2.1b)"
       print *, "================================"

       do i=0, m2
          do j=0, n2
            T1(i, j) = init
          end do
       end do
       scale = 0.05

       do i=0, m2
          call uu(i*scale, ymin, T2(i, 0))
          call uu(i*scale, ymax, T2(i, n2))
       end do

       do j=1, n2
          call uu(xmin, j*scale, T2(0, j))
          call uu(xmax, j*scale, T2(m2, j))
       end do

       do i=0, m2
          do j=0, n2
             call uu(i*scale, j*scale, error_matrix2(i, j))
          end do
       end do     

       call Laplace_Solver(T2, error_matrix2, m2, n2, iter, ERROR_SUM)
       print *
       print *, "Number of iterations:" , iter
       print *, "Culmulative error:   ", ERROR_SUM
       print *

       print *
       print *, "================================"
       print *, "Problem 2.2a)"
       print *, "================================"

       boundary_sum = 0.0

       scale = 0.1
       do i=0, m1
          call uu(i*scale, ymin, T1(i, 0))
          call uu(i*scale, ymax, T1(i, n1))
          boundary_sum = boundary_sum +  T1(i, 0) + T1(i, n1)
       end do

       do j=1, n1-1
          call uu(xmin, j*scale, T1(0, j))
          call uu(xmax, j*scale, T1(m1, j))
          boundary_sum = boundary_sum + T1(0, j) + T1(m1, j)
       end do

       denom = ((n1+1)+(m1-1))*2
       avg_boundary = boundary_sum/denom

       do i=1, m1-1
          do j=1, n1-1
            T1(i, j) = avg_boundary
          end do
       end do
       do i=0, m1
          do j=0, n1
             call uu(i*scale, j*scale, error_matrix1(i, j))
          end do
       end do          
       call Laplace_Solver(T1, error_matrix1,m1, n1, iter, ERROR_SUM)
       print *
       print *, "Number of iterations:" , iter
       print *, "Culmulative error:   ", ERROR_SUM
       print *

end program Laplace_Equation


subroutine Laplace_Solver(PHIOLD, ERR_MAT, MTOT, NTOT, ITER, ERROR_SUM)

      implicit none

      integer MTOT,NTOT
      double precision, dimension (0:MTOT, 0:NTOT):: PHINEW, PHIOLD, ERR_MAT
      integer ITER
      double precision ERROR_SUM

      integer i, j
      double precision TEMP,  CONVG
      real TESTS, tmp

      CONVG = 0.000003
      ERROR_SUM = 0
      ITER = 0

     print *
     print *, "Initial internal state"
     do 15 i=0, MTOT
 15  write(6, 200) (PHIOLD(i, j), j=0, NTOT)

 10    TESTS = 0
       ITER = ITER+1
 
       do i=1, MTOT-1
         do j=1, NTOT-1
           TEMP=0.25*(PHIOLD(i,j+1)+PHIOLD(i,j-1)+PHIOLD(i+1,j)+PHIOLD(i-1,j))
           tmp =  abs(TEMP - PHIOLD(i,j))
           TESTS = amax1(TESTS, tmp)
           PHINEW(i,j) = TEMP

         end do
       end do

       do i=1, MTOT-1
         do j=1, NTOT-1

           PHIOLD(i,j) = PHINEW(i,j)

         end do
       end do

       if (TESTS > CONVG) go to 10

       do i=1, MTOT-1
          do j=1, NTOT-1
            ERROR_SUM = ERROR_SUM + abs(PHIOLD(i, j) - ERR_MAT(i, j))
          end do
       end do      

     print *
     print * 
     print *, "Result"
     do 20 i=0, MTOT
 20  write(6, 200) (PHIOLD(i, j), j=0, NTOT)
 200 format(11F9.4)

end subroutine Laplace_Solver

subroutine uu(x, y, result)
       
  double precision, intent(in) :: x, y
  double precision, intent(out):: result

  result = -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 subroutine uu


