function [ nval, ival ] = points_hull_2d ( node_num, node_xy )

%% POINTS_HULL_2D computes the convex hull of a set of points in 2D.
%
%  Modified:
%
%    03 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Output, integer NVAL, the number of nodes that lie on the convex hull.
%
%    Output, integer IVAL(NODE_NUM).  The first NVAL entries of IVAL contain
%    the indices of the nodes that form the convex hull, in order.
%    These indices are 1-based, not 0-based!
%
  dim_num = 2;

  if ( node_num < 1 )
    nval = 0;
    return
  end
%
%  If NODE_NUM = 1, the hull is the point.
%
  if ( node_num == 1 )
    nval = 1;
    ival(1) = 1;
    return
  end
%
%  If NODE_NUM = 2, then the convex hull is either the two distinct points,
%  or possibly a single (repeated) point.
%
  if ( node_num == 2 )

    nval = 1;
    ival(1) = 1;

    if ( any ( node_xy(1:dim_num,1) ~= node_xy(1:dim_num,2) ) )
      nval = nval + 1;
      ival(2) = 2;
    end

    return

  end
%
%  Find the leftmost point, and take the bottom-most in a tie.
%  Call it "Q".
%
  iq = 1;
  for i = 2 : node_num
    if ( node_xy(1,i) < node_xy(1,iq) | ...
       ( node_xy(1,i) == node_xy(1,iq) & node_xy(2,i) < node_xy(2,iq) ) )
      iq = i;
    end
  end
%
%  Remember the starting point.
%
  istart = iq;
  nval = 1;
  ival(1) = iq;
%
%  For the first point, make a dummy previous point, 1 unit south,
%  and call it "P".
%
  pp(1:2) = [ node_xy(1,iq), node_xy(2,iq) - 1.0 ];
%
%  Now, having old point P, and current point Q, find the new point R
%  so the angle PQR is maximal.
%
%  Watch out for the possibility that the two nodes are identical.
%
  while ( 1 )

    ir = 0;
    angmax = 0.0;

    for i = 1 : node_num

      if ( i ~= iq & ( node_xy(1,i) ~= node_xy(1,iq) | node_xy(2,i) ~= node_xy(2,iq) ) )

        ai = angle_rad_2d ( node_xy(1:dim_num,i), node_xy(1:dim_num,iq), pp );

        if ( ir == 0 | angmax < ai )

          ir = i;
          angmax = ai;
%
%  In case of ties, choose the nearer point.
%
        elseif ( ir ~= 0 & ai == angmax )

          di = sqrt ( sum ( ( node_xy(1:dim_num,i)  - node_xy(1:dim_num,iq) ).^2 ) );
          dr = sqrt ( sum ( ( node_xy(1:dim_num,ir) - node_xy(1:dim_num,iq) ).^2 ) );

          if ( di < dr )
            ir = i;
            angmax = ai;
          end

        end

      end

    end
%
%  If we've returned to our starting point, exit.
%
    if ( ir == istart )
      break
    end

    if ( node_num < nval + 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'POINTS_HULL_2D - Fatal error!\n' );
      fprintf ( 1, '  The algorithm failed.\n' );
      fprintf ( 1, '  We have not returned to the starting point,\n' );
      fprintf ( 1, '  and we have already used all the points.\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, ' NODE_NUM = NVAL = %d\n', nval );
      error ( 'POINTS_HULL_2D - Fatal error!' );
    end
%
%  Set Q := P, P := R, and repeat.
%
    nval = nval + 1;
    ival(nval) = ir;

    pp(1:dim_num) = node_xy(1:dim_num,iq);

    iq = ir;

  end