function [V,K,D,M,G,dG,GP] = BEM2dLaplace(nodes,edges,order,varargin)

load GaussQuadratureForLineSegments.mat

tau  = points{order};
w = weights{order};
nG = length(w);

[S,T] = meshgrid(tau,tau);

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

V = zeros(m,m);
K = zeros(m,n);
Curl = K;
Mdata = zeros(1,2*m);

G = zeros(m,1);
dG = zeros(n,1);
GP = zeros(n,1);
if length(varargin)==2
    f = varargin{1};
    nf = length(f);
    xf = varargin{2};
    for j=1:nf
        GP = GP - f(j) * (1/2/pi) * log(sqrt((xf(1,j)-nodes(1,:)).^2+(xf(2,j)-nodes(2,:)).^2))';
    end
else
    nf = 0;
end

for i=1:m
    edgei = edges(:,i);
    nodesi = nodes(:,edgei);
    xi1 = nodesi(:,1);
    xi2 = nodesi(:,2);
    ti = xi2-xi1;
    Si = norm(ti);
    
    if nf>0
        yi = xi1*ones(1,nG)+ti*tau;
        ni = [ti(2);-ti(1)]/Si;
        for j=1:nf
            xfjyi = xf(:,j)*ones(1,nG)-yi;
            normxfjyi = sqrt(xfjyi(1,:).^2+xfjyi(2,:).^2);
            G(i) = G(i) + f(j) * Si * w * log(normxfjyi)';
            dG(edgei(1)) = dG(edgei(1)) - f(j) * Si * w * ( (1-tau).*((ni(1)*xfjyi(1,:)+ni(2)*xfjyi(2,:)))./normxfjyi.^2 )';
            dG(edgei(2)) = dG(edgei(2)) - f(j) * Si * w * ( tau.*((ni(1)*xfjyi(1,:)+ni(2)*xfjyi(2,:)))./normxfjyi.^2 )';
        end
    end
    
%     M(i,edgei) = Si*[0.5,0.5];
    Curl(i,edgei) = [-1,1]/Si;
    idx = 2*(i-1)+[1,2];
    Mdata(idx) = Si/2;
    
    for j=1:m
        if i==j % identical segments
            V(i,i) = Si^2 * (log(Si)-1.5);
        else
            edgej = edges(:,j);
            nodesj = nodes(:,edgej);
            xj1 = nodesj(:,1);
            xj2 = nodesj(:,2);
            tj = xj2-xj1;
            Sj = norm(tj);
            nj = [tj(2);-tj(1)]/Sj;
            
            commonVertex = intersect(edgei,edgej);
            if ~isempty(commonVertex) % common vertex
                idxi1 = find(edgei==commonVertex);
                idxi2 = find(edgei~=commonVertex);
                idxj2 = find(edgej~=commonVertex);
                a = nodesi(:,idxi2)-nodesi(:,idxi1);
                b = nodesj(:,idxj2)-nodesi(:,idxi1);
                
                F1 = (a(1)-b(1)*tau).^2+(a(2)-b(2)*tau).^2;
                F2 = (a(1)*tau-b(1)).^2+(a(2)*tau-b(2)).^2;
                F = 0.5*log( F1 .* F2 );
                V(i,j) = Si*Sj * 0.5 * ( -1 + w*F' );

                K(i,edgei(idxi1)) = K(i,edgei(idxi1)) - Si*Sj*(a'*nj) * ( w*((1-0.5*tau)./F1 + 0.5*tau./F2)' );
                K(i,edgej(idxj2)) = K(i,edgej(idxj2)) - Si*Sj*(a'*nj) * ( w*(0.5*tau./F1 + 0.5*tau./F2)' );
                
            else % disjoint segments
                a = xi1-xj1;

                AB1 = a(1)+ti(1)*S-tj(1)*T;
                AB2 = a(2)+ti(2)*S-tj(2)*T;
                norm2 = AB1.^2+AB2.^2;
                F = 0.5 * log( norm2 );
                V(i,j) = Si*Sj * w*F*w';
                
                F = (AB1*nj(1)+AB2*nj(2)) ./ norm2;
                K(i,edgej(1)) = K(i,edgej(1)) - Si*Sj * w*(F.*(1-T))*w';
                K(i,edgej(2)) = K(i,edgej(2)) - Si*Sj * w*(F.*T)*w';
                
            end
            
        end
        
    end
    
end

V = V / (-2*pi);
K = K / (-2*pi);

D = Curl'*V*Curl;

rows = reshape([1;1]*[1:m],2*m,1);
M = sparse(rows,edges,Mdata,m,n);

G = G / (-2*pi);
dG = dG / (-2*pi);

