
?- op(700, xfx, ->).

isN(16).
isEps(1.0E-6).

main(Res) :- gridDomain(UInit), initialVals(UInit), loop(0, UInit, Res).

% `gridDomain' asserts that a mapping (``array'') has an N x N grid as
% domain (``index space'').  NB: it imposes a particular ordering---
% reverse lexicographic---on sites in the mapping list.

gridDomain(U) :- isN(N), enumDomain(N, N, U).

enumDomain(0, 0, [site(0, 0) -> _])     :- !.
enumDomain(I, 0, [site(I, 0) -> _ | U]) :-
			!, IM is I - 1, isN(N), enumDomain(IM, N, U).
enumDomain(I, J, [site(I, J) -> _ | U]) :-
			JM is J - 1, enumDomain(I, JM, U).

% `initialVals' asserts that a mapping defines a suitable initial state.

initialVals([]).
initialVals([site(I, J) -> B | UTail]) :-
    isN(N),
    (I =:= 0 ; I =:= N ; J =:= 0 ; J =:= N), !,
    X is float(I) / float(N),
    Y is float(J) / float(N),
    bdry(X, Y, B),
    initialVals(UTail).
initialVals([_ -> 0.0 | UTail]) :- initialVals(UTail).

% The main loop...

loop(Niter, U, Res) :- gridDomain(UNew),
                       update(U, UNew),
                       Miter is Niter + 1,
                       loopTest(Miter, U, UNew, Res).

% `update' asserts that all mappings in 2nd arg equal the average of
% values of neighbour mappings in 1st arg.

update(_, []).
update(Uold, [site(I, J) -> X | UTail]) :-
    isN(N),
    (I =:= 0 ; I =:= N ; J =:= 0 ; J =:= N), !,
    element(site(I, J) -> X, Uold),
    update(Uold, UTail).
update(Uold, [site(I, J) -> X | UTail]) :-
    JM is J - 1, element(site(I, JM) -> A, Uold),
    JP is J + 1, element(site(I, JP) -> B, Uold),
    IM is I - 1, element(site(IM, J) -> C, Uold),
    IP is I + 1, element(site(IP, J) -> D, Uold),
    X  is (A + B + C + D) / 4.0,
    update(Uold, UTail).

% Check for convergence.

%loopTest(Niter, OldU, U, res(Niter, U)) :- converged(OldU, U), !.
loopTest(Niter, OldU, U, Niter) :- converged(OldU, U), !.
loopTest(Niter, _,    U, Res)           :- loop(Niter, U, Res).

converged([], []).
converged([_ -> X | U], [_ -> Y | V]) :- Dif is abs(X - Y),
					 isEps(Eps), Dif < Eps,
					 converged(U, V).

bdry(X, Y, B) :- B is X * X - Y * Y.

element(X, [X | _]) :- !.
element(X, [_ | L]) :- element(X, L).

