function V = operatorV2d(G,P,E,rad)

load GaussQuadratureForLineSegments

w = weights{rad};
ksi = points{rad};
N = length(w);
[ksiAlfa,ksiBeta] = meshgrid(ksi,ksi);

n = size(E,2);
V = 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);
    isc = intersect(Ei,Ej);
    if length(isc)==0
      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;
      V(i,j) = Si*Sj * w * G(sqrt((X1-Y1).^2+(X2-Y2).^2)) * w';
%    else if length(isc)==1
%      ...
    end
  end
end
