PROGRAM verify
*     verify the discrete form of Laplace's equation
*     for a point charge on a square lattice
      IMPLICIT NONE
      REAL*8 L, dx, dy
      INTEGER IC1_, IC2_, IC3_, IR_
      CHARACTER stop_plot*80
      LOGICAL Ltmp1_
      DATA IC1_/1/
      DATA IC2_/2/
      DATA IC3_/3/
      CALL GWopen(IR_, 0)
      CALL initial(L,dx,dy,IC1_,IC2_)
      CALL draw_grid(L,dx,dy,IC3_)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL potential(L,dx,dy,stop_plot,IC1_,IC2_,IC3_)
         Ltmp1_ = stop_plot .NE. 'yes'
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L,dx,dy,IC1_,IC2_)
      IMPLICIT NONE
      REAL*8 L, dx, dy
      REAL GWaspect
      INTEGER IR_, IC1_, IC2_
      WRITE(*,'(A,$)') 'lattice dimension = '
      READ(*,*) L
      WRITE(*,'(A,$)') 'grid spacing in x direction = '
      READ(*,*) dx
      WRITE(*,'(A,$)') 'grid spacing in y direction = '
      READ(*,*) dy
      CALL GWvport(IR_, 0.0, 0.8, 0.25*GWaspect(1), 1.0)
      CALL GWindow(IR_, -1.2, -1.0, 2.0, 2.0)
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWvport(IR_, 0.0, 0.0, 0.25*GWaspect(1), 0.6)
      CALL GWindow(IR_, 0.0, 24.0, 80.0, 0.0);
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      END
C
      SUBROUTINE draw_grid(L,dx,dy,IC1_)
      IMPLICIT NONE
      REAL*8 x, y, L, dx, dy
      REAL a, xwin, ywin, GWaspect
      INTEGER IR_, IC1_
      CALL GWvport(IR_, 0.25*GWaspect(1), 0.0, GWaspect(1), 1.0)
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      a = 0.2*dx                            ! "radius" of visual image of charge
      CALL GWsetmrk(IR_, 6, 2*a, 16, -1, -1)
      CALL GWputmrk(IR_, 0.0, 0.0)
      CALL GWsetmrk(IR_, 6, a/2, 0, -1, -1)
      DO y = -L, L, dy
         DO x = -L, L, dx
            CALL GWputmrk(IR_, REAL(x), REAL(y))
         END DO
      END DO
      CALL GWrect(IR_, REAL(-L), REAL(-L), REAL(L), REAL(L))
      END
C
      SUBROUTINE potential(L,dx,dy,stop_plot,IC1_,IC2_,IC3_)
      IMPLICIT NONE
      REAL*8 V0, V1, V2, V3, V4, L, dx, dy
      REAL xs, ys
      INTEGER IR_, IC1_, IC2_, IC3_
      CHARACTER*80 stop_plot, Stmp1_
      CALL GWselvp(IR_, IC3_)
      CALL GWsetogn(IR_, 0, IC3_)
      CALL GWsetmsg(IR_, 'click on mouse to choose site')
      CALL GWcappnt(IR_, xs, ys, '')
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWerase(IR_, IC1_, -1)
      stop_plot = ''
      IF(abs(xs) .LE. L .AND. abs(ys) .LE. L) THEN
         CALL showpotential(DBLE(xs), DBLE(ys), 0, 0, V0)
         CALL showpotential(xs+dx, DBLE(ys), 1, 0, V1)
         CALL showpotential(xs-dx, DBLE(ys), -1, 0, V2)
         CALL showpotential(DBLE(xs), ys+dy, 0, 1, V3)
         CALL showpotential(DBLE(xs), ys-dy, 0, -1, V4)
         CALL GWselvp(IR_, IC2_)
         CALL GWsetogn(IR_, 0, IC2_)
         CALL GWerase(IR_, IC2_, -1)
         CALL GWputtxt(IR_, 1.0, 1.0, ' average potential')
         CALL GWputtxt(IR_, 1.0, 2.0, ' of four neighbors:')
         WRITE(Stmp1_,'(F8.4)') 0.25*(V1+V2+V3+V4)
         CALL GWputtxt(IR_, 1.0, 3.0, Stmp1_)
      ELSE
         stop_plot = 'yes'
      END IF
      END
C
      SUBROUTINE showpotential(x,y,xp,yp,V)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 x, y, V
      INTEGER IR_, xp, yp
      V = 1/sqrt(x*x + y*y)                 ! potential of a point charge
      WRITE(Stmp1_,'(F8.4)') V
      CALL GWputtxt(IR_, REAL(xp), REAL(yp), Stmp1_)
      END