function [x,y,meshele,meshnod]=trigrid(dx,dy,lx,ly,angle)
% [x,y,meshele,meshnod]=trigrid(dx,dy,lx,ly,angle)
% make a parallelogram of connected triangles spanning
% lx,ly  with spacing dx,dy
% optional argument of angle is defaulted to 30
% angle = 0-90 degrees.  
% 0 square;  30 equilaterals; 90 is flattened beyond reason ugly
%  [x,y,meshele,meshnod]=trigrid(800,800,20800,20800,30);
if nargin==4
angle = 30;
end
cang = cos(angle*pi/180);
sang = sin(angle*pi/180);

ly = lx*sqrt(3)/2
dy = dx;

n = fix(lx/dx) ;
m = fix(ly/dy) ;
m = n;

x = ones(n*m,1);
y = ones(n*m,1);
meshele=ones(n*m*2,4);
mk = 0;

for i = 1:(n+1)
for j = 1:(m+1)
e1= (i-1)*(n+1)+j;
%x(e1) = (i-1)*dx + (j-1)*dy*sang;
%y(e1) = (j-1)*dy ;
x(e1) = (i-1)*dx + (j-1)*dy*sang;
y(e1) = -(i-1)*dx*0 + (j-1)*dy*cang ;
end
end
meshnod = [[1:length(x)]' x y];

for i = 1:n
for j = 1:m
e1= (i-1)*(n+1)+j;
e2= (i+1-1)*(n+1)+j ;
e3= (i-1)*(n+1)+j+1 ;
e4= (i+1-1)*(n+1)+j+1 ;

mk = mk+1;
meshele(mk,:)= [ mk e1 e2 e3];
mk = mk+1;
meshele(mk,:)=  [ mk e4 e3 e2];

end
end


[n,m]=size(meshnod);

slope =  (meshnod(1,3)-meshnod(2,3))/(meshnod(1,2)-meshnod(2,2)) 

for nod1 = 1:n
%   function to define bad points
%  h^2+ly^2 = 4h^2     h^2 = ly^2/3  h = sqrt*(1/3) ly/2
%   y = -sqrt*(2/3) *xn  + ly
%    Line going up parallel to coast (45deg) and intersecting X at -1.6e4
yn=meshnod(nod1,3);
xn = meshnod(nod1,2);   
nodebad = (yn < (-slope *xn  + ly -dy/2 ) )  | ...
          (yn > (-slope *xn  + 3*ly +dy/2 ) ) > 0   ;
if nodebad
%   Eliminate all meshele which use nod1
%      (or keep all meshele which do not use nod1)
   l =find(meshele(:,2)~=nod1 & meshele(:,3)~=nod1 & meshele(:,4)~=nod1 ) ;
   meshele=meshele(l,:);
   meshele(:,1) = [1:length(l)]';
end
end
xnl = 0:dx:(.55*lx);
xnu = (.95*lx):dx:(1.3*lx);
slope =  (meshnod(1,3)-meshnod(2,3))/(meshnod(1,2)-meshnod(2,2)) 
ynl = (-slope *xnl  + ly -dy/2 );
ynu = (-slope *xnu  + 3*ly +dy/2 );
x3 = meshele(:,2:4);
x3(:)=x(meshele(:,2:4));
y3 = meshele(:,2:4);
y3(:)=y(meshele(:,2:4));
fill(x3',y3','w')
hold on; plot(xnl,ynl,'r',xnu,ynu,'b');hold off
axis('equal')
grid
