function pyramid_unit_quad_test ( degree_max )

%*****************************************************************************80
%
%% pyramid_unit_quad_test() tests the rules for the unit pyramid.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer DEGREE_MAX, the maximum total degree of the
%    monomials to check.
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'pyramid_unit_quad_test():\n' );
  fprintf ( 1, '  for the unit pyramid,\n' );
  fprintf ( 1, '  we approximate monomial integrals with:\n' );
  fprintf ( 1, '  pyramid_unit_o01(),\n' );
  fprintf ( 1, '  pyramid_unit_o05(),\n' );
  fprintf ( 1, '  pyramid_unit_o06(),\n' );
  fprintf ( 1, '  pyramid_unit_o08(),\n' );
  fprintf ( 1, '  pyramid_unit_o08b(),\n' );
  fprintf ( 1, '  pyramid_unit_o09(),\n' );
  fprintf ( 1, '  pyramid_unit_o13(),\n' );
  fprintf ( 1, '  pyramid_unit_o18(),\n' );
  fprintf ( 1, '  pyramid_unit_o27(),\n' );
  fprintf ( 1, '  pyramid_unit_o48().\n' );

  expon = [];
  more = 0;
  h = 0;
  t = 0;

  while ( true )

    [ expon, more, h, t ] = subcomp_next ( degree_max, dim_num, expon, more, h, t );

    if ( mod ( expon(1), 2 ) == 1 || ...
         mod ( expon(2), 2 ) == 1 )
      continue
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Monomial exponents: ' );
    for ( dim = 1 : dim_num )
      fprintf ( 1, '  %2d', expon(dim) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    order = 1;
    [ w, xyz ] = pyramid_unit_o01 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) *  w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 5;
    [ w, xyz ] = pyramid_unit_o05 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 6;
    [ w, xyz ] = pyramid_unit_o06 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 8;
    [ w, xyz ] = pyramid_unit_o08 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 8;
    [ w, xyz ] = pyramid_unit_o08b ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 9;
    [ w, xyz ] = pyramid_unit_o09 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 13;
    [ w, xyz ] = pyramid_unit_o13 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 18;
    [ w, xyz ] = pyramid_unit_o18 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 27;
    [ w, xyz ] = pyramid_unit_o27 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    order = 48;
    [ w, xyz ] = pyramid_unit_o48 ( );
    v = monomial_value ( dim_num, order, expon, xyz );
    quad = pyramid_unit_volume ( ) * w' * v;
    fprintf ( 1, '  %6d  %14f\n', order, quad );

    fprintf ( 1, '\n' );
    quad = pyramid_unit_monomial ( expon );
    fprintf ( 1, '   Exact  %14f\n', quad );

    if ( ~ more )
      break;
    end

  end

  return
end
