function plot_field(no2xy,el2no,el2ed,sol)
%  
el2no = sort(el2no);
% for interpolation: note sum k (phi_k(1:n)) = ones(1:n)
phi_1 = [4 3 2 1 0 3 2 1 0 2 1 0 1 0 0]'/4;
phi_2 = [0 1 2 3 4 0 1 2 3 0 1 2 0 1 0]'/4;
phi_3 = [0 0 0 0 0 1 1 1 1 2 2 2 3 3 4]'/4;

p1 = no2xy(:,el2no(1,:));
p2 = no2xy(:,el2no(2,:));
p3 = no2xy(:,el2no(3,:));

edge1 = p2 - p1;
edge2 = p3 - p1;
detJ  = edge1(1,:).*edge2(2,:) - edge1(2,:).*edge2(1,:);
grad_phi_2x =  edge2(2,:)./detJ;
grad_phi_2y = -edge2(1,:)./detJ;
grad_phi_3x = -edge1(2,:)./detJ;
grad_phi_3y =  edge1(1,:)./detJ;
grad_phi_1x = 0 - grad_phi_2x - grad_phi_3x; % phi1+phi2+phi3 = 1
grad_phi_1y = 0 - grad_phi_2y - grad_phi_3y;

sol1 = sol(el2ed(1,:)).';
sol2 = sol(el2ed(2,:)).';
sol3 = sol(el2ed(3,:)).';

% jwH = curl E; H = H ez (scalar)
% compare the formulas for the S matrix
Hz   = (sol1 - sol2 + sol3)./(detJ);

% E = sum sol(j) Nj, linearly interpolated, barycentric coord., 
% to finer mesh
Ex = phi_1*( grad_phi_2x.*sol1 + grad_phi_3x.*sol2) + ...
     phi_2*(-grad_phi_1x.*sol1 + grad_phi_3x.*sol3) + ...
     phi_3*(-grad_phi_1x.*sol2 - grad_phi_2x.*sol3);
%=== bug fixed: two last lines should have phi_2 and phi_3, not phi_1 !!
Ey = phi_1*( grad_phi_2y.*sol1 + grad_phi_3y.*sol2) + ...
     phi_2*(-grad_phi_1y.*sol1 + grad_phi_3y.*sol3) + ...
     phi_3*(-grad_phi_1y.*sol2 - grad_phi_2y.*sol3);

% refined mesh
psub = kron(p1,phi_1') + kron(p2,phi_2') + kron(p3,phi_3');
ih = ishold;
ax = newplot;
%shading flat
mrz = max(abs(Hz(:)));
caxis([-mrz,mrz])
patch('faces',el2no','vertices',no2xy','facevertexcdata',Hz(:),...
      'facecolor','flat',...
      'edgecolor',get(ax,'defaultsurfaceedgecolor'));
axis equal, hold on
quiver(psub(1,:),psub(2,:),Ex(:)',Ey(:)','k')
xy = [p1;p2;p3;p1;NaN*p1];
plot(xy(1:2:end),xy(2:2:end),'k')

%c = (0:64)'/64;
%d = [c c ones(size(c))];
%colormap([d; 1 1 1; d(end:-1:1,end:-1:1)])
if ~ih
    hold off
end





