function order = rule_order ( p )

%*****************************************************************************80
%
%% rule_order() returns the order of a pyramid quadrature rule of given precision.
%
%  Discussion:
%
%    The unit pyramid with square base is the region
%
%      -1 <= X <= 1
%      -1 <= Y <= 1
%       0 <= Z <= 1.
%
%    The "order" of a quadrature rule is the number of points involved.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 April 2023
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jan Jaskowiec, Natarajan Sukumar,
%    High order symmetric cubature rules for tetrahedra and pyramids,
%    International Journal for Numerical Methods in Engineering,
%    Volume 122, Number 1, pages 148-171, 24 August 2020.
%
%  Input:
%
%    integer p: the precision, 0 <= p <= 20.
%
%  Output:
%
%    integer order: the order of the rule.
%
  if ( p < 0 ) 
    error ( 'rule_order(): Input "p" < 0.\n' );
  end

  if ( 20 < p )
    error ( 'rule_order(): Input 20 < "p".\n' );
  end

  order_vec = [ ...
      1, ...
      1,   5,   6,  10,  15, ...
     23,  31,  47,  62,  80, ...
    103, 127, 152, 184, 234, ...
    285, 319, 357, 418, 489 ];

  order = order_vec(p+1);

  return
end

