function L = edge_edge_laplacian(V,E,F)

  n = length(V);
  f = size(F,1);
  k = length(E);
  L = sparse(k, k);
  
%  reverse_edges=sparse(n,n);
%  for i=1:k
%    reverse_edges(E(i,1),E(i,2)) = i;
%    reverse_edges(E(i,2),E(i,1)) = i;
%  end
  
  reverse_edges = sparse(E(:,1), E(:,2), 1:k, k, k);
  reverse_edges = reverse_edges + reverse_edges';
  
  C = 4 * cotangent(V, F);
  
  for i=1:f
    
    % Vertex indices
    iV1 = F(i, 1);
    iV2 = F(i, 2);
    iV3 = F(i, 3);
  
    % Edges
    E1 = [iV2, iV3];
    E2 = [iV3, iV1];
    E3 = [iV1, iV2];
    
    % Edge indices
    iE1 = reverse_edges(iV2,iV3);
    iE2 = reverse_edges(iV3,iV1);
    iE3 = reverse_edges(iV1,iV2);
    
%    % Edge vector
%    D0 = V(iV3, :) - V(iV2, :);
%    D1 = V(iV1, :) - V(iV3, :);
%    D2 = V(iV2, :) - V(iV1, :);
%    
%    % Edge angles
%    costheta2 = -D2 * D1' / norm(D2) / norm(D1);
%    costheta3 = -D0 * D2' / norm(D0) / norm(D2);
%    cosTheta3 = -D1 * D0' / norm(D1) / norm(D0);
%    
%    theta1 = acos(costheta2);
%    theta2 = acos(costheta3);
%    theta3 = acos(cosTheta3);
%    
%    cotantheta1 = 1.0/tan(theta1);
%    cotantheta2 = 1.0/tan(theta2);
%    cotantheta3 = 1.0/tan(theta3);
%    
%    % Laplacian matrix filling
%    L(iE1, iE1) = L(iE1, iE1) + 2*cotantheta3 + 2*cotantheta2;
%    L(iE1, iE2) = L(iE1, iE2) - 2*cotantheta3;
%    L(iE1, iE3) = L(iE1, iE3) - 2*cotantheta2;
%    
%    L(iE2, iE2) = L(iE2, iE2) + 2*cotantheta1 + 2*cotantheta3;
%    L(iE2, iE3) = L(iE2, iE3) - 2*cotantheta1;
%    L(iE2, iE1) = L(iE2, iE1) - 2*cotantheta3;
%    
%    L(iE3, iE3) = L(iE3, iE3) + 2*cotantheta2 + 2*cotantheta1;
%    L(iE3, iE1) = L(iE3, iE1) - 2*cotantheta2;
%    L(iE3, iE2) = L(iE3, iE2) - 2*cotantheta1;
    
    
    % Laplacian matrix filling
    L(iE1, iE1) = L(iE1, iE1) + C(i,3) + C(i,2);
    L(iE1, iE2) = L(iE1, iE2) - C(i,3);
    L(iE1, iE3) = L(iE1, iE3) - C(i,2);
    
    L(iE2, iE2) = L(iE2, iE2) + C(i,1) + C(i,3);
    L(iE2, iE3) = L(iE2, iE3) - C(i,1);
    L(iE2, iE1) = L(iE2, iE1) - C(i,3);
    
    L(iE3, iE3) = L(iE3, iE3) + C(i,2) + C(i,1);
    L(iE3, iE1) = L(iE3, iE1) - C(i,2);
    L(iE3, iE2) = L(iE3, iE2) - C(i,1);
    
  endfor
  
  
endfunction
