function M = edge_mass_matrix(V,E,F)

  n = length(V);
  f = size(F,1);
  k = size(E,1);
  M = 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';
  
  for i=1:f
    
    % Vertex indices
    iV0 = F(i, 1);
    iV1 = F(i, 2);
    iV2 = F(i, 3);
  
    % Edges
    E0 = [iV1, iV2];
    E1 = [iV2, iV0];
    E2 = [iV0, iV1];
    
    % Edge indices
    iE0 = reverse_edges(iV1,iV2);
    iE1 = reverse_edges(iV2,iV0);
    iE2 = reverse_edges(iV0,iV1);
    
    % Face Area
    A = norm(cross([V(iV1, :) - V(iV0, :), 0], [V(iV2, :) - V(iV0, :),0]));
    
    M(iE0, iE0) = M(iE0, iE0) + A/3.0;
    M(iE1, iE1) = M(iE0, iE0) + A/3.0;
    M(iE2, iE2) = M(iE0, iE0) + A/3.0;
    
    
    
  endfor
  
  
endfunction
