function Vkappa = operatorV2dHelmholtz(kappa,P,E,rad)

V0 = operatorV2dLaplace(P,E,rad);

load GaussQuadratureForLineSegments

Gtilde = @(r) ((i/4)*besselh(0,kappa*r)+log(r)/2/pi);

w = weights{rad};
ksi = points{rad};
w2 = weights{rad+1};
ksi2 = points{rad+1};
[ksiAlfa,ksiBeta] = meshgrid(ksi,ksi2);

n = size(E,2);
Vrem = zeros(n,n);
for i=1:n
  Ei = E(:,i);
  Pi0 = P(:,Ei(1));
  Pi1 = P(:,Ei(2));
  Si = norm(Pi1-Pi0);
  for j=1:n
    Ej = E(:,j);
    Pj0 = P(:,Ej(1));
    Pj1 = P(:,Ej(2));
    Sj = norm(Pj1-Pj0);
    X1 = Pi0(1) + (Pi1(1)-Pi0(1))*ksiAlfa;
    X2 = Pi0(2) + (Pi1(2)-Pi0(2))*ksiAlfa;
    Y1 = Pj0(1) + (Pj1(1)-Pj0(1))*ksiBeta;
    Y2 = Pj0(2) + (Pj1(2)-Pj0(2))*ksiBeta;
    Vrem(i,j) = Si*Sj * w2 * Gtilde(sqrt((X1-Y1).^2+(X2-Y2).^2)) * w';
  end
end

Vkappa = V0 + Vrem;
