function [ boundary_path, boundary_edge ] = boundary_edge_to_path ( ...
  boundary_edge )

%*****************************************************************************80
%
%% boundary_edge_to_path() constructs a path of nodes around a triangulation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 August 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer BOUNDARY_EDGE(BOUNDARY_EDGE_NUM,2), the pairs of nodes
%    that make up the boundary edges.
%
%  Output:
%
%    integer BOUNDARY_PATH(BOUNDARY_EDGE_NUM,1), a list of nodes that
%    form a counterclockwise path around the boundary.
%
%    integer BOUNDARY_EDGE(BOUNDARY_EDGE_NUM,2), reordered to conform
%    to the path.
%
  boundary_edge_num = size ( boundary_edge, 1 );

  boundary_path = zeros ( boundary_edge_num, 1 );

  k = 1;
  e = 1;
  boundary_path(k) = boundary_edge(1,1);
 
  while ( boundary_edge(e,2) ~= boundary_path(1) )
    k = k + 1;
    boundary_path(k) = boundary_edge(e,2);

    found = false;
    for e2 = e + 1 : boundary_edge_num
      if ( boundary_edge(e2,1) == boundary_edge(e,2) )
        found = true;
        break;
      end
    end

    if ( ~ found )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'boundary_edge_to_path - Fatal error!\n' );
      fprintf ( 1, '  Stuck at node %d\n', boundary_path(k) );
      error ( 'boundary_edge_to_path - Fatal error!\n' );
    end

    boundary_edge(e+1:boundary_edge_num,1:2) = [ ...
      boundary_edge(e2,1:2);
      boundary_edge(e+1:e2-1,1:2);
      boundary_edge(e2+1:boundary_edge_num,1:2)];

    e = e + 1;

  end

  return
end

