function c = triangle_circumcenter ( tri, x, y )

%*****************************************************************************80
%
%% triangle_circumcenter() returns the centers of circumcircles of triangles.
%
%  Modified:
%
%    12 June 2006
%
%  Reference: 
%
%    Watson, p32.
%
%  Input:
%
%    integer TRI(N,3), the indices of nodes that make up the triangles.
%
%    real X(*), Y(*), the coordinates of nodes.
%
%  Output:
%
%    real C(N,2), the coordinates of the centers of the
%    circles that circumscribe each triangle.
%
  x1 = x(tri(:,1));
  x2 = x(tri(:,2));
  x3 = x(tri(:,3));

  y1 = y(tri(:,1));
  y2 = y(tri(:,2));
  y3 = y(tri(:,3));
%
%  Set equation for center of each circumcircle: 
%  [a11 a12;a21 a22]*[x;y] = [b1;b2] * 0.5;
%
  a11 = x2 - x1;
  a12 = y2 - y1;
  a21 = x3 - x1;
  a22 = y3 - y1;
%
%  Solve the 2-by-2 equation explicitly.
%
  idet = a11 .* a22 - a21 .* a12;
%
%  Add small random displacement to points that are either the same
%  or on a line.
%
  d = find ( idet == 0 );
%
%  Add small random displacement to points.
%
  if ( ~isempty ( d ) )

    delta = sqrt ( eps );
    x1(d) = x1(d) + delta * ( rand ( size ( d ) ) - 0.5 );
    x2(d) = x2(d) + delta * ( rand ( size ( d ) ) - 0.5 );
    x3(d) = x3(d) + delta * ( rand ( size ( d ) ) - 0.5 );
    y1(d) = y1(d) + delta * ( rand ( size ( d ) ) - 0.5 );
    y2(d) = y2(d) + delta * ( rand ( size ( d ) ) - 0.5 );
    y3(d) = y3(d) + delta * ( rand ( size ( d ) ) - 0.5 );

    a11 = x2 - x1;
    a12 = y2 - y1;
    a21 = x3 - x1;
    a22 = y3 - y1;

    idet = a11 .* a22 - a21 .* a12;

  end

  b1 = a11 .* (x2+x1) + a12 .* (y2+y1);
  b2 = a21 .* (x3+x1) + a22 .* (y3+y1);

  idet = 0.5 ./ idet;

  xcenter = ( a22 .* b1 - a12 .* b2) .* idet;
  ycenter = (-a21 .* b1 + a11 .* b2) .* idet;

  c = [ xcenter ycenter ];

  return
end
