function pwl_interp_2d_test01 ( prob, n )

%*****************************************************************************80
%
%% pwl_interp_2d_test01 tests pwl_interp_2d.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer PROB, the problem number.
%
%    integer N, the grid size in each dimension.
%
  nxd = n;
  nyd = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_2D_TEST01:\n' );
  fprintf ( 1, '  Interpolate data from TEST_INTERP_2D problem #%d\n', prob );
  fprintf ( 1, '  Using polynomial interpolant of degree %d x %d\n', nxd, nyd );

  nd = nxd * nyd;
  fprintf ( 1, '  Number of data points = %d\n', nd );

  xd_1d = linspace ( 0.0, 1.0, nxd );
  yd_1d = linspace ( 0.0, 1.0, nyd );

  xd = zeros (nxd*nyd,1);
  yd = zeros (nxd*nyd,1);

  ij = 0;
  for j = 1 : nyd
    for i = 1 : nxd
      ij = ij + 1;
      xd(ij) = xd_1d(i);
      yd(ij) = yd_1d(j);
    end
  end

  zd = f00_f0 ( prob, nd, xd, yd );

  if ( nd <= 20 )
    r8vec3_print ( xd, yd, zd, '  X, Y, Z data:' );
  end
%
%  #1:  Does interpolant match function at data points?
%
  ni = nd;
  xi(1:ni) = xd(1:ni);
  yi(1:ni) = yd(1:ni);
  zd = reshape ( zd, nxd, nyd );
  zi = pwl_interp_2d ( nxd, nyd, xd_1d, yd_1d, zd, ni, xi, yi );

  if ( ni <= 20 )
    r8vec3_print ( xi, yi, zi, '  X, Y, Z interpolation:' );
  end

  zd = reshape ( zd, nxd * nyd, 1 );
  int_error = r8vec_norm_affine ( ni, zi, zd ) / ni;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  RMS data interpolation error = %g\n', int_error );
%
%  #2:  Does interpolant approximate data at midpoints?
%
  if ( 1 < nd )

    xi_1d(1:nxd-1) = 0.5 * ( xd_1d(1:nxd-1) + xd_1d(2:nxd) );
    yi_1d(1:nyd-1) = 0.5 * ( yd_1d(1:nyd-1) + yd_1d(2:nyd) );

    ni = ( nxd - 1 ) * ( nyd - 1 );

    xi = zeros ( ni, 1 );
    yi = zeros ( ni, 1 );
   
    ij = 0;
    for j = 1 : nyd - 1
      for i = 1 : nxd - 1
        ij = ij + 1;
        xi(ij) = xi_1d(i);
        yi(ij) = yi_1d(j);
      end
    end

    zdm = f00_f0 ( prob, ni, xi, yi );

    zd = reshape ( zd, nxd, nyd );
    zi = pwl_interp_2d ( nxd, nyd, xd_1d, yd_1d, zd, ni, xi, yi );
    zd = reshape ( zd, nxd * nyd, 1 );
    
    app_error = r8vec_norm_affine ( ni, zi, zdm ) / ni;

    fprintf ( 1, '  RMS data approximation error = %g\n', app_error );

  end

  return
end
 
