function [areas,normals] = triangleareas( nodes, triangles, whichtriangles )
%[areas,normals] = triangleareas( nodes, triangles )
%   Compute the areas and normal vectors to a set of triangles specified by
%   the positions of their vertexes and the triples of vertexes that make
%   the triangles.

    if nargin < 3
        v1 = nodes( triangles(:,1), : );
        v2 = nodes( triangles(:,2), : );
        v3 = nodes( triangles(:,3), : );
    else
        v1 = nodes( triangles(whichtriangles,1), : );
        v2 = nodes( triangles(whichtriangles,2), : );
        v3 = nodes( triangles(whichtriangles,3), : );
    end
    if size(nodes,2) < 3
        z = zeros( size(v1,1), 1 );
        v1 = [ v1, z ];
        v2 = [ v2, z ];
        v3 = [ v3, z ];
    end
    v12 = v2-v1;
    v13 = v3-v1;
    normals = cross( v12, v13, 2 );
    areas = sqrt( sum( normals.*normals, 2 ) );
    for i=1:size(normals,2)
        normals(:,i) = normals(:,i) ./ areas;
    end
    areas = areas/2;
end
