function [V,K,D,M] = BEM2dHelmholtz(kappa,nodes,edges,order)

load GaussQuadratureForLineSegments.mat

tau  = points{order};
w = weights{order};

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

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

V = zeros(m,m);
K = zeros(m,n);
D = zeros(n,n);
M = K;
Curl = K;

for ii=1:m
    edgei = edges(:,ii);
    nodesi = nodes(:,edgei);
    xi1 = nodesi(:,1);
    xi2 = nodesi(:,2);
    ti = xi2-xi1;
    Si = norm(ti);
    ni = [ti(2);-ti(1)]/Si;
    
    M(ii,edgei) = Si*[0.5,0.5];
    Curl(ii,edgei) = [-1,1]/Si;
    
    for j=1:m
        if ii==j % identical segments
            F = 2*besselh(0,kappa*Si*S.*(1-T));
            V(ii,ii) = Si^2 * (w*(S.*F)*w');
            
            F11 = (1-S).*(1-S.*T);
            F12 = 0.5*((1-S).*S.*T+(1-S.*T).*S);
            F22 = S.*S.*T;
            D(edgei(1),edgei(1)) = D(edgei(1),edgei(1)) + Si^2 * (w*(S.*F.*F11)*w');
            D(edgei(1),edgei(2)) = D(edgei(1),edgei(2)) + Si^2 * (w*(S.*F.*F12)*w');
            D(edgei(2),edgei(1)) = D(edgei(2),edgei(1)) + Si^2 * (w*(S.*F.*F12)*w');
            D(edgei(2),edgei(2)) = D(edgei(2),edgei(2)) + Si^2 * (w*(S.*F.*F22)*w');

        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;
            ninj = ni'*nj;
            
            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);

                AB1a = a(1)-b(1)*T;
                AB2a = a(2)-b(2)*T;
                norma = sqrt(AB1a.^2+AB2a.^2);
                F1 = besselh(0,kappa*S.*norma);
                AB1b = a(1)*T-b(1);
                AB2b = a(2)*T-b(2);
                normb = sqrt(AB1b.^2+AB2b.^2);
                F2 = besselh(0,kappa*S.*normb);
                V(ii,j) = Si*Sj * (w*(S.*(F1+F2))*w');
                
                F11 = (1-S).*(1-S.*T);
                F12 = (1-S).*S.*T;
                F21 = S.*(1-S.*T);
                F22 = S.*S.*T;
                D(edgei(idxi1),edgei(idxi1)) = D(edgei(idxi1),edgei(idxi1)) + Si*Sj*ninj * (w*(S.*(F1+F2).*F11)*w');
                D(edgei(idxi1),edgej(idxj2)) = D(edgei(idxi1),edgej(idxj2)) + Si*Sj*ninj * (w*(S.*(F1.*F12+F2.*F21))*w');
                D(edgei(idxi2),edgei(idxi1)) = D(edgei(idxi2),edgei(idxi1)) + Si*Sj*ninj * (w*(S.*(F1.*F21+F2.*F12))*w');
                D(edgei(idxi2),edgej(idxj2)) = D(edgei(idxi2),edgej(idxj2)) + Si*Sj*ninj * (w*(S.*(F1+F2).*F22)*w');
                
                F1 = besselh(1,kappa*S.*norma) .* (nj(1)*AB1a+nj(2)*AB2a)./norma;
                F2 = besselh(1,kappa*S.*normb) .* (nj(1)*AB1b+nj(2)*AB2b)./normb;
                K(ii,edgei(idxi1)) = K(ii,edgei(idxi1)) + Si*Sj* ( w*(S.*(F1.*(1-S.*T)+F2.*(1-S)))*w' );
                K(ii,edgej(idxj2)) = K(ii,edgej(idxj2)) + Si*Sj* ( w*(S.*(F1.*S.*T+F2.*S))*w' );

            else % disjoint segments
                a = xi1-xj1;
  
                AB1 = a(1)+ti(1)*S-tj(1)*T;
                AB2 = a(2)+ti(2)*S-tj(2)*T;
                norm1 = sqrt(AB1.^2+AB2.^2);
                F = besselh(0,kappa*norm1);
                V(ii,j) = Si*Sj * (w*F*w');
                
                F11 = (1-S).*(1-T);
                F12 = (1-S).*T;
                F21 = S.*(1-T);
                F22 = S.*T;
                D(edgei(1),edgej(1)) = D(edgei(1),edgej(1)) + Si*Sj*ninj * (w*(F.*F11)*w');
                D(edgei(1),edgej(2)) = D(edgei(1),edgej(2)) + Si*Sj*ninj * (w*(F.*F12)*w');
                D(edgei(2),edgej(1)) = D(edgei(2),edgej(1)) + Si*Sj*ninj * (w*(F.*F21)*w');
                D(edgei(2),edgej(2)) = D(edgei(2),edgej(2)) + Si*Sj*ninj * (w*(F.*F22)*w');
                
                F = besselh(1,kappa*norm1) .* (nj(1)*AB1+nj(2)*AB2) ./ norm1;
                K(ii,edgej(1)) = K(ii,edgej(1)) + Si*Sj * w*(F.*(1-T))*w';
                K(ii,edgej(2)) = K(ii,edgej(2)) + Si*Sj * w*(F.*T)*w';

            end

        end
        
    end
    
end

V = (i/4) * V;
K = (i/4) * kappa*K;
D = Curl'*V*Curl - (i/4)*kappa^2*D;
