function xy = ellipse_grid ( n, r, c, ng )

%*****************************************************************************80
%
%% ELLIPSE_GRID generates grid points inside an ellipse.
%
%  Discussion:
%
%    The ellipse is specified as
%
%      ( ( X - C1 ) / R1 )^2 + ( ( Y - C2 ) / R2 )^2 = 1
%
%    The user supplies a number N.  There will be N+1 grid points along
%    the shorter axis.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, real R(2), the half axis lengths.
%
%    Input, real C(2), the center of the ellipse.
%
%    Input, integer NG, the number of grid points inside the ellipse.
%
%    Output, real XY(2,NG), the grid points.
%
  if ( r(1) < r(2) )
    h = 2 * r(1) / ( 2 * n + 1 );
    ni = n;
    nj = ceil ( r(2) / r(1) ) * n;
  else
    h = 2 * r(2) / ( 2 * n + 1 );
    nj = n;
    ni = ceil ( r(1) / r(2) ) * n;
  end

  p = 0;

  for j = 0 : nj
    i = 0;
    x = c(1);
    y = c(2) + j * h;
    p = p + 1;
    xy(1:2,p) = [ x, y ]';
    if ( 0 < j )
      p = p + 1;
      xy(1:2,p) = [ x, 2 * c(2) - y ]';
    end
    while ( 1 )
      i = i + 1;
      x = c(1) + i * h;
      if ( 1 < ( ( x - c(1) ) / r(1) ).^2 + ( ( y - c(2) ) / r(2) ).^2 )
        break
      end
      p = p + 1;
      xy(1:2,p) = [ x, y ]';
      p = p + 1;
      xy(1:2,p) = [ 2 * c(1) - x, y ]';
      if ( 0 < j )
        p = p + 1;
        xy(1:2,p) = [ x, 2 * c(2) - y ]';
        p = p + 1;
        xy(1:2,p) = [ 2 * c(1) - x, 2 * c(2) - y ]';
      end
    end
  end

  return
end