function [lgood, lbad] = flatmesh(meshele,meshnod,frac)
% [lgood, lbad] = flatmesh(meshele,meshnod,frac);
% This checks the elements of meshele
%  Display all triangle elements which are too flat
%  Too flat means containing an angle which is 
%  either less than frac*pi  or greater (1-frac)*pi
disp(['mesh/meshtools/flatmesh'])


x = meshnod(meshele(:,2),2) - meshnod(meshele(:,3),2);
y = meshnod(meshele(:,2),3) - meshnod(meshele(:,3),3);
a2 = x.*x + y.*y;
a = sqrt(a2);
x = meshnod(meshele(:,3),2) - meshnod(meshele(:,4),2);
y = meshnod(meshele(:,3),3) - meshnod(meshele(:,4),3);
b2 = x.*x + y.*y;
b = sqrt(b2);
x = meshnod(meshele(:,4),2) - meshnod(meshele(:,2),2);
y = meshnod(meshele(:,4),3) - meshnod(meshele(:,2),3);
c2 = x.*x + y.*y;
c = sqrt(c2);
alpha = acos((b2+c2-a2)./(2*b.*c));
beta = acos((a2+c2-b2)./(2*a.*c));
gam = acos((a2+b2-c2)./(2*a.*b));
ang = (1-frac)*pi;
angs = frac*pi;

% No angles greater than .98pi or less than .02pi 
lbad = find(alpha>ang | beta>ang | gam>ang |...
            alpha<angs | beta<angs | gam<angs);
lgood = 1:length(meshele(:,1));
lgood(lbad) = lgood(lbad)*NaN; 
lgood = lgood(find(~isnan(lgood)));


%ltri = meshele(lbad(1),[2:4 2])
%plot(meshnod(ltri,2),meshnod(ltri,3),'r')

