!****************************************************************
!Example Usage:
!gfortran LaplaceEq.f90 -o l
!(echo 31;echo 100 -100;echo 1e-6)|./l
!gnuplot> splot "data" w l not
!gnuplot> set hidden3d;set pm3d;replot
!Check convergence: (no. iters) ~ L^2
!L=31 ; for e in 1e-4 1e-6 1e-8 1e-10 1e-12 1e-14 1e-16;do (echo $L;echo 100 -100;echo $e)|./l | tail -n 1;done
!L=61 ; for e in 1e-4 1e-6 1e-8 1e-10 1e-12 1e-14 1e-16;do (echo $L;echo 100 -100;echo $e)|./l | tail -n 1;done
!L=128; for e in 1e-4 1e-6 1e-8 1e-10 1e-12 1e-14 1e-16;do (echo $L;echo 100 -100;echo $e)|./l | tail -n 1;done
!****************************************************************
program laplace_em
 use, intrinsic                       :: iso_fortran_env
 implicit none
 integer ,parameter                   :: dp = real64
 logical ,allocatable, dimension(:,:) :: isConductor
 real(dp),allocatable, dimension(:,:) :: V
 real(dp)                             :: V1,V2,epsilon
 integer                              :: L
 real(dp),parameter                   :: ZERO=0.0_dp,QUARTER=0.25_dp

 print *,'Enter L:'
 read  *, L
 print *,'Enter V1,V2:'
 read  *, V1,V2
 print *,'Enter epsilon:'
 read  *, epsilon
 print *,'Starting Laplace:'
 print *,'Grid Size= ',L
 print *,'Conductors set at V1= ',V1,' V2= ',V2
 print *,'Relaxing with accuracy epsilon= ',epsilon
!The arrays V and isConductor are initialized
 call initialize_lattice
!We enter initialized V,isConductor. On exit the
!routine gives the solution V
 call laplace
!We print V in a file.
 call print_results

!****************************************************
!----------------------------------------------------
contains
!----------------------------------------------------
!****************************************************
 subroutine     initialize_lattice
  integer                :: i,j

  ALLOCATE(isConductor(L,L), V(L,L))

!Initialize to 0 and .FALSE (default values for boundary and
!interior sites).
  V                          = ZERO
  isConductor                = .FALSE.
!We set the boundary to be a conductor: (V=0 by default)
  isConductor(      1,    :) = .TRUE.
  isConductor(      :,    1) = .TRUE.
  isConductor(      L,    :) = .TRUE.
  isConductor(      :,    L) = .TRUE.

  V          (  L/3+1,5:L-5) = V1
  isConductor(  L/3+1,5:L-5) = .TRUE.

  V          (2*L/3+1,5:L-5) = V2
  isConductor(2*L/3+1,5:L-5) = .TRUE.
  
 end subroutine initialize_lattice
!****************************************************************
 subroutine     laplace
  integer                :: i,j,icount
  real(dp)                :: Vav,error,dV

  icount  = 0                !counts number of sweeps
  do while (.TRUE.)          !an infinite loop:
   error  = ZERO             !Exit when error<epsilon
   do  j  = 2,L-1
    do i  = 2,L-1
!We change V only for non conductors:
     if( .NOT. isConductor(i,j))then
      Vav = ( V(i-1,j)+V(i+1,j)+V(i,j+1)+V(i,j-1)) * QUARTER
      dV  = DABS(V(i,j)-Vav)
      if(error .LT. dV ) error = dV !maximum error
      V(i,j) = Vav          ! we immendiately update V(i,j)
     end if
    end do
   end do
   icount = icount + 1
   print *,icount,' err= ',error
   if( error .LT. epsilon) return !return to main program
  end do

 end subroutine laplace
!****************************************************************
 subroutine     print_results
  integer                :: i,j, unit

  open(newunit=unit,file="data")
  do  i = 1,L
   do j = 1,L
    write(unit,*)i,j,V(i,j)
   end do
   write (unit,*)'' !print empty line for gnuplot, separate isolines
  end do

end subroutine print_results
!****************************************************
!----------------------------------------------------
end program laplace_em
!----------------------------------------------------
!****************************************************