function [ x, y ] = wathen_xy ( nx, ny, n )

%*****************************************************************************80
%
%% wathen_xy(): coordinates of Wathen nodes.
%
%  Discussion:
%
%    We will take the region to be the unit square.
%
%    The grid uses quadratic serendipity elements.
%
%    Here is an illustration of the node numbering for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 February 2020
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Input:
%
%    integer NX, NY, values which determine the size of the matrix.
%
%    integer N, the number of variables.
%
%  Output:
%
%    real X(N), Y(N), the node coordinates.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'wathen_xy - Fatal error!\n' );
    fprintf ( 1, '  Not enough input.\n' );
    error ( 'wathen_xy - Fatal error!' );
  end

  if ( nargin < 2 )
    ny = nx;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NY was not supplied.  Setting NY = NX = %d.\n', ny );
  end

  if ( nargin < 3 )
    n = wathen_order ( nx, ny );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  N was not supplied.  N = %d\n', n );
  end

  x = zeros ( n, 1 );
  y = zeros ( n, 1 );

  k = 0;
  for j = 1 : 2 * ny + 1
    if ( mod ( j, 2 ) == 1 )
      x(k+1:k+2*nx+1) = linspace ( 0.0, 1.0, 2*nx + 1 );
      y(k+1:k+2*nx+1) = ( j - 1 ) / ( 2 * ny ); 
      k = k + 2 * nx + 1; 
    else
      x(k+1:k+nx+1) = linspace ( 0.0, 1.0, nx + 1 ); 
      y(k+1:k+nx+1) = ( j - 1 ) / ( 2 * ny );
      k = k + nx + 1; 
    end
  end

  return
end

