function [l,angle] = len(node,elem)
%[l,angle] = len(node,elem)
%
%     l: Array with the elemental lengths
% angle: Array with the elemental angles
%  node: Nodal positions
%   elem: Element connectivities

%The number of elements
nel=size(elem,1);

for i=1:nel
    %Find the node numbers
    node1=elem(i,1);
    node2=elem(i,2);
    
    %Nodal positions
    x1=node(node1,1);
    y1=node(node1,2);
    x2=node(node2,1);
    y2=node(node2,2);
    
    %The length of the element
    l(i)=sqrt((x2-x1)^2+(y2-y1)^2);
    
    %The angle of the element, in radians
    
    
    if and(x1<x2,y1==y2) %+ horizontal
        angle(i)=0;        
    elseif and(x1>x2,y1==y2) %- horizontal
        angle(i)=pi;
    elseif and(x1==x2,y1<y2) %+ vertical
        angle(i)=pi/2;
    elseif and(x1==x2,y1>y2) %- vertical
        angle(i)=3*pi/2;
    elseif and(x1<x2,y1>y2) % ccw left boundary
        angle(i)=2*pi-asin((y1-y2)/l(i));
    elseif and(x1<x2,y1<y2) % ccw right boundary
        angle(i)=asin((y1-y2)/l(i));
    elseif and(x1>x2,y1<y2) % cw left boundary
        angle(i)=pi-asin((y1-y2)/l(i));  
    elseif and(x1>x2,y1>y2) % cw right boundary
        angle(i)=pi+asin((y1-y2)/l(i));        
    end
    
end
