PROGRAM verify
! verify the discrete form of Laplace's equation
! for a point charge on a square lattice
LIBRARY "csgraphics"
CALL initial(L,dx,dy,#1,#2)
CALL draw_grid(L,dx,dy,#3)
DO
   CALL potential(L,dx,dy,stop_plot$,#1,#2,#3)
LOOP until stop_plot$ = "yes"
END

SUB initial(L,dx,dy,#1,#2)
    INPUT prompt "lattice dimension = ": L
    INPUT prompt "grid spacing in x direction = ": dx
    INPUT prompt "grid spacing in y direction = ": dy
    OPEN #1: screen 0,0.25,0.8,1
    SET WINDOW -1.2,2,-1,2
    OPEN #2: screen 0,0.25,0,0.6
END SUB

SUB draw_grid(L,dx,dy,#1)
    OPEN #1: screen 0.25,1,0,1
    CALL compute_aspect_ratio(L,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    LET a = 0.2*dx                ! "radius" of visual image of charge
    SET COLOR "blue"
    BOX CIRCLE -a,a,-a,a
    FLOOD 0,0
    SET COLOR "black"
    FOR y = -L to L step dy
        FOR x = -L to L step dx
            PLOT POINTS: x,y
        NEXT x
    NEXT y
    BOX LINES -L,L,-L,L
END SUB

SUB potential(L,dx,dy,stop_plot$,#1,#2,#3)
    WINDOW #3
    SET CURSOR 1,20
    PRINT "click on mouse to choose site"
    LET s = 0
    DO
       GET MOUSE xs,ys,s
    LOOP until s = 2
    WINDOW #1
    CLEAR
    LET stop_plot$ = ""
    IF abs(xs) <= L and abs(ys) <= L then
       CALL showpotential(xs,ys,0,0,V0)
       CALL showpotential(xs+dx,ys,1,0,V1)
       CALL showpotential(xs-dx,ys,-1,0,V2)
       CALL showpotential(xs,ys+dy,0,1,V3)
       CALL showpotential(xs,ys-dy,0,-1,V4)
       WINDOW #2
       SET CURSOR 1,1
       PRINT " average potential"
       PRINT " of four neighbors:"
       PRINT truncate(0.25*(V1+V2+V3+V4),4)
    ELSE
       LET stop_plot$ = "yes"
    END IF
END SUB

SUB showpotential(x,y,xp,yp,V)
    LET V = 1/sqr(x*x + y*y)      ! potential of a point charge
    LET V = truncate(V,4)
    PLOT TEXT, AT xp,yp: str$(V)
END SUB