function bc = baryCoords( vxs, n, v )
%bc = baryCoords( vxs, n, v )
%   Calculate the barycentric coordinates of point v with respect to
%   triangle vxs, in three dimensions.  vxs contains the three vertexes of
%   the triangle as rows. n is a normal vector to the triangle.
%   The component of v perpendicular to the plane of the
%   triangle is ignored.
%   If n is the empty array, the normal vector will be computed.  If
%   supplied, the normal vector may be of any non-zero length.
%   n and v must be row vectors, and the result is a 3-element row vector.
%   v can also be a matrix of row vectors; the result will then be a matrix
%   of 3-element row vectors.
%   This procedure works in 3 or more dimensions.
%   For 2 dimensions, use baryCoordsPlanar( vxs, v ).

    if isempty(n)
        n = trinormal( vxs );
    end
    
    v1 = vxs(1,:);
    v2 = vxs(2,:);
    v3 = vxs(3,:);
    vv = zeros(size(v));
    nsq = n*n';
    for i=1:size(v,1)
        vv(i,:) = v(i,:) + (((v1-v(i,:))*n')/nsq)*n;  % vv is in the plane of the triangle.
    end
    if false
        try
            vv = v + (((v1-v)*n')/nsq)*n;  % vv is in the plane of the triangle.
        catch
            xxx = 0;
        end
    end
  % vv = v + (dot3(v1-v,n)/(n*n'))*n;  % v is now in the plane of the triangle.
    
    a1 = cross3( v2-v3, n );
    oneval = dot3(a1,v1);
    zeroval = dot3(a1,v2);
    vval = vv*a1';
    b1 = (vval-zeroval)/(oneval-zeroval);
    
    a2 = cross3( v3-v1, n );
    oneval = dot3(a2,v2);
    zeroval = dot3(a2,v3);
    vval = vv*a2';
    b2 = (vval-zeroval)/(oneval-zeroval);
    
    a3 = cross3( v1-v2, n );
    oneval = dot3(a3,v3);
    zeroval = dot3(a3,v1);
    vval = vv*a3';
    b3 = (vval-zeroval)/(oneval-zeroval);
    
    bc = [ b1, b2, b3 ];
    
  % baryCoords_err = norm(bc*vxs-vv)
  % vv
end

function d = dot3( v1, v2 )
%v3 = dot3( v1, v2 )
%   Dot product of two 3-element row vectors.

    d = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3);
end

function v3 = cross3( v1, v2 )
%v3 = dot3( v1, v2 )
%   Dot product of two 3-element row vectors.

    v3 = [ v1(2)*v2(3) - v1(3)*v2(2), ...
           v1(3)*v2(1) - v1(1)*v2(3), ...
           v1(1)*v2(2) - v1(2)*v2(1) ];
end
