function V = operatorV2dLaplace(P,E,rad)

load GaussQuadratureForLineSegments

G = @(r) (-log(r)/2/pi);

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';
    elseif length(isc)==1
      P0i = P(:,isc);
      P1i = P(:,setdiff(Ei,isc));
      P1j = P(:,setdiff(Ej,isc));
      u = P1i-P0i;
      v = P1j-P0i;
      V(i,j) = 1/4/pi;
      xy1 = u(1)-ksi*v(1);
      xy2 = u(2)-ksi*v(2);
      normxy = sqrt(xy1.^2+xy2.^2);
      V(i,j) = V(i,j) + (-1/2/pi) * (1/2) * w * G(normxy)';
      xy1 = ksi*u(1)-v(1);
      xy2 = ksi*u(2)-v(2);
      normxy = sqrt(xy1.^2+xy2.^2);
      V(i,j) = V(i,j) + (-1/2/pi) * (1/2) * w * G(normxy)';
      V(i,j) = Si*Sj * V(i,j);
    else
      V(i,i) = - (log(Si)-3/2) * Si^2 / 2 / pi;
    end
  end
end
