function [Vt,Wu,G] = BEM2dLaplace_potentials(nodes,edges,order,t,u,volPoints,varargin);

load GaussQuadratureForLineSegments.mat

tau  = points{order};
etau = 1-tau;
w = weights{order};

m = size(edges,2);
n = size(volPoints,2);

Vt = zeros(1,n);
Wu = zeros(1,n);

G = zeros(1,n);
if length(varargin)==2
    f = varargin{1};
    nf = length(f);
    xf = varargin{2};
else
    nf = 0;
end

for i=1:n
    xi1 = volPoints(:,i);

    if nf>0
        for j=1:nf
            G(i) = G(i) + f(j) * log(norm(xf(:,j)-xi1));
        end
    end

    for j=1:m
        edgej = edges(:,j);
        nodesj = nodes(:,edgej);
        xj1 = nodesj(:,1);
        xj2 = nodesj(:,2);
        tj = xj2-xj1;
        Sj = norm(tj);
        nj = [tj(2);-tj(1)]/Sj;

        a = xi1-xj1;

        AB1 = a(1)-tj(1)*tau;
        AB2 = a(2)-tj(2)*tau;
        norm1 = sqrt(AB1.^2+AB2.^2);
        F = log(norm1);
        Vt(i) = Vt(i) + t(j) * Sj*w*F';

        F = - (nj(1)*AB1+nj(2)*AB2) ./ norm1.^2;
        Wu(i) = Wu(i) + u(edgej(1)) * Sj*w*(F.*etau)';
        Wu(i) = Wu(i) + u(edgej(2)) * Sj*w*(F.*tau)';
        
    end
    
end

Vt = Vt/(-2*pi);
Wu = Wu/(-2*pi);
G = G/(-2*pi);
