function [Vt,Wu] = BEM2dHelmholtz_potentials(kappa,nodes,edges,order,t,u,volPoints);

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);

one = ones(1,m);
Tau = tau'*one;
Etau = etau'*one;
onetau = ones(1,length(tau));

for i=1:n
    xi1 = volPoints(:,i);
    
%    X1 = nodes(:,edges(1,:));
%    X2 = nodes(:,edges(2,:));
%    T = X2-X1;
%    S = sqrt(T(1,:).^2+T(2,:).^2);
%    N = [T(2,:)./S;-T(1,:)./S];
%    
%    A = xi1*one-X1;
%    
%    AB1 = onetau'*A(1,:)-tau'*T(1,:);
%    AB2 = onetau'*A(2,:)-tau'*T(2,:);
%    Norm1 = sqrt(AB1.^2+AB2.^2);
%    F = besselh(0,kappa*Norm1);
%    Vt(i) = (t.*S)*transpose(w*F);
%    
%    F = besselh(1,kappa*Norm1) .* ((onetau'*N(1,:)).*AB1+(onetau'*N(2,:)).*AB2) ./ Norm1;
%    Wu(i) = (transpose(u(edges(1,:))).*S)*transpose(w*(F.*Etau));
%    Wu(i) = Wu(i) + (transpose(u(edges(2,:))).*S)*transpose(w*(F.*Tau));
    
     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 = besselh(0,kappa*norm1);
         Vt(i) = Vt(i) + t(j) * Sj*w*F';
 
         F = besselh(1,kappa*norm1) .* (nj(1)*AB1+nj(2)*AB2) ./ norm1;
         Wu(i) = Wu(i) + u(edgej(1)) * Sj*w*transpose(F.*etau);
         Wu(i) = Wu(i) + u(edgej(2)) * Sj*w*transpose(F.*tau);        
     end
    
end

Vt = (i/4) * Vt;
Wu = (i/4) * kappa*Wu;
