function L = edge_edge_normal_derivative(V,E,F)

  n = size(V,1);
  f = size(F,1);
  k = size(E,1);
  L = sparse(k, k);
%  L = sym('n', [k,k]);
%  L(:) = 0;
  
%  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



%  C = sym('c', size(F));
  %  C = 4*C;
  
  C = 4 * cotangent(V, F);

  [~,on_b] = on_boundary(F);
  [on_b_f,on_b_c] = find(on_b);
  
  reverse_edges = sparse(E(:,1), E(:,2), 1:k, k, k);
  reverse_edges = reverse_edges + reverse_edges';
  
  f_on_b=size(on_b_f,1);
  
%  R = sym('r',[k,k]);
%  R = R * 2;
%  EE=[   
%   0   0   1   0   1   0   0   0
%   0   0   1   0   0   0   1   0
%   1   1   0   0   1   0   1   0
%   0   0   0   0   1   0   0   1
%   1   0   1   1   0   0   0   1
%   0   0   0   0   0   0   1   1
%   0   1   1   0   0   1   0   1
%   0   0   0   1   1   1   1   0];
%   
%  R = R.*EE;
  
  for i=1:f_on_b
    
    % Vertex subscript
    iiV0 = mod(on_b_c(i)+1,3)+1;
    iiV1 = mod(on_b_c(i)+2,3)+1; % = on_b_c(i)
    iiV2 = mod(on_b_c(i)  ,3)+1;
    
    % Vertex indices
    iV0 = F(on_b_f(i), iiV0);
    iV1 = F(on_b_f(i), iiV1); % not on the boundary
    iV2 = F(on_b_f(i), iiV2);
  
    % Edges
    E0 = [iV1, iV2];
    E1 = [iV2, iV0]; % on the boundary
    E2 = [iV0, iV1];
    
    % Edge indices
    iE0 = reverse_edges(iV1,iV2);
    iE1 = reverse_edges(iV2,iV0);
    iE2 = reverse_edges(iV0,iV1);
    
%    % Edge vector
%    D0 = V(iV2, :) - V(iV1, :);
%    D1 = V(iV0, :) - V(iV2, :);
%    D2 = V(iV1, :) - V(iV0, :);
%    
%    % Edge angles
%    cosTheta0 = -D2 * D1' / norm(D2) / norm(D1);
%    cosTheta1 = -D0 * D2' / norm(D0) / norm(D2);
%    cosTheta2 = -D1 * D0' / norm(D1) / norm(D0);
%    
%    theta0 = acos(cosTheta0);
%    theta1 = acos(cosTheta1);
%    theta2 = acos(cosTheta2);
%    
%    cotanTheta0 = 1.0/tan(theta0);
%    cotanTheta1 = 1.0/tan(theta1);
%    cotanTheta2 = 1.0/tan(theta2);
%    
%    cotanTheta0 = 2*cotanTheta0;
%    cotanTheta1 = 2*cotanTheta1;
%    cotanTheta2 = 2*cotanTheta2;
    
    
    % Laplacian matrix filling

    L(iE1, iE0) = L(iE1, iE0) + C(on_b_f(i), iiV2);
    L(iE1, iE1) = L(iE1, iE1) - (C(on_b_f(i), iiV0) + C(on_b_f(i), iiV2));
    L(iE1, iE2) = L(iE1, iE2) + C(on_b_f(i), iiV0);
    
  endfor
  
endfunction
