# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "pyramid_witherden_rule.hpp"

//****************************************************************************80

void comp_next ( int n, int k, int a[], bool &more, int &h, int &t )

//****************************************************************************80
//
//  Purpose:
//
//    comp_next() computes the compositions of the integer N into K parts.
//
//  Discussion:
//
//    A composition of the integer N into K parts is an ordered sequence
//    of K nonnegative integers which sum to N.  The compositions (1,2,1)
//    and (1,1,2) are considered to be distinct.
//
//    The routine computes one composition on each call until there are no more.
//    For instance, one composition of 6 into 3 parts is
//    3+2+1, another would be 6+0+0.
//
//    On the first call to this routine, set MORE = FALSE.  The routine
//    will compute the first element in the sequence of compositions, and
//    return it, as well as setting MORE = TRUE.  If more compositions
//    are desired, call again, and again.  Each time, the routine will
//    return with a new composition.
//
//    However, when the LAST composition in the sequence is computed 
//    and returned, the routine will reset MORE to FALSE, signaling that
//    the end of the sequence has been reached.
//
//    This routine originally used a SAVE statement to maintain the
//    variables H and T.  I have decided that it is safer
//    to pass these variables as arguments, even though the user should
//    never alter them.  This allows this routine to safely shuffle
//    between several ongoing calculations.
//
//
//    There are 28 compositions of 6 into three parts.  This routine will
//    produce those compositions in the following order:
//
//     I         A
//     -     ---------
//     1     6   0   0
//     2     5   1   0
//     3     4   2   0
//     4     3   3   0
//     5     2   4   0
//     6     1   5   0
//     7     0   6   0
//     8     5   0   1
//     9     4   1   1
//    10     3   2   1
//    11     2   3   1
//    12     1   4   1
//    13     0   5   1
//    14     4   0   2
//    15     3   1   2
//    16     2   2   2
//    17     1   3   2
//    18     0   4   2
//    19     3   0   3
//    20     2   1   3
//    21     1   2   3
//    22     0   3   3
//    23     2   0   4
//    24     1   1   4
//    25     0   2   4
//    26     1   0   5
//    27     0   1   5
//    28     0   0   6
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    02 July 2008
//
//  Author:
//
//    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms for Computers and Calculators,
//    Second Edition,
//    Academic Press, 1978,
//    ISBN: 0-12-519260-6,
//    LC: QA164.N54.
//
//  Parameters:
//
//    Input, int N, the integer whose compositions are desired.
//
//    Input, int K, the number of parts in the composition.
//
//    Input/output, int A[K], the parts of the composition.
//
//    Input/output, bool &MORE.
//    Set MORE = FALSE on first call.  It will be reset to TRUE on return
//    with a new composition.  Each new call returns another composition until
//    MORE is set to FALSE when the last composition has been computed
//    and returned.
//
//    Input/output, int &H, &T, two internal parameters needed for the
//    computation.  The user should allocate space for these in the calling
//    program, include them in the calling sequence, but never alter them!
//
{
  int i;

  if ( !( more ) )
  {
    t = n;
    h = 0;
    a[0] = n;
    for ( i = 1; i < k; i++ )
    {
       a[i] = 0;
    }
  }
  else
  {
    if ( 1 < t )
    {
      h = 0;
    }
    h = h + 1;
    t = a[h-1];
    a[h-1] = 0;
    a[0] = t - 1;
    a[h] = a[h] + 1;
  }

  more = ( a[k-1] != n );

  return;
}
//****************************************************************************80

double *monomial_value ( int m, int n, int e[], double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_value() evaluates a monomial.
//
//  Discussion:
//
//    This routine evaluates a monomial of the form
//
//      product ( 1 <= i <= m ) x(i)^e(i)
//
//    The combination 0.0^0 is encountered is treated as 1.0.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    17 August 2014
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the spatial dimension.
//
//    int N, the number of evaluation points.
//
//    int E[M], the exponents.
//
//    double X[M*N], the point coordinates.
//
//  Output:
//
//    double MONOMIAL_VALUE[N], the monomial values.
//
{
  int i;
  int j;
  double *v;

  v = new double[n];
  for ( j = 0; j < n; j++)
  {
    v[j] = 1.0;
  }

  for ( i = 0; i < m; i++ )
  {
    if ( 0 != e[i] )
    {
      for ( j = 0; j < n; j++ )
      {
        v[j] = v[j] * pow ( x[i+j*m], e[i] );
      }
    }
  }

  return v;
}
/******************************************************************************/

double pyramid_unit_monomial_integral ( int expon[3] )

/******************************************************************************/
/*
  Purpose:

    pyramid_unit_monomial_integral(): monomial integral in a unit pyramid.

  Discussion:

    This function returns the value of the integral of X^ALPHA Y^BETA Z^GAMMA
    over the unit pyramid.

    The integration region is:

    - ( 1 - Z ) <= X <= 1 - Z
    - ( 1 - Z ) <= Y <= 1 - Z
              0 <= Z <= 1.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 March 2008

  Author:

    John Burkardt

  Reference:

    Arthur Stroud,
    Approximate Calculation of Multiple Integrals,
    Prentice Hall, 1971,
    ISBN: 0130438936,
    LC: QA311.S85.

  Input:

    int EXPON[3], the exponents.

  Output:

    double pyramid_unit_monomial_integral, the integral of the monomial
    over the pyramid.
*/
{
  int i;
  int i_hi;
  double mop;
  double value;

  value = 0.0;

  if ( ( expon[0] % 2 ) == 0 && ( expon[1] % 2 ) == 0 )
  {
    i_hi = 2 + expon[0] + expon[1];

    mop = 1.0;
    for ( i = 0; i <= i_hi; i++ )
    {
      value = value + mop * r8_choose ( i_hi, i )
      / ( double ) ( i + expon[2] + 1 );
      mop = - mop;
    }

    value = value
          * 2.0 / ( double ) ( expon[0] + 1 )
          * 2.0 / ( double ) ( expon[1] + 1 );
  }

  return value;
}
//****************************************************************************80

double pyramid_unit_volume ( )

//****************************************************************************80
//
//  Purpose:
//
//    pyramid_unit_volume(): volume of a unit pyramid with square base.
//
//  Discussion:
//
//    The volume of this unit pyramid is 4/3.
//
//    The integration region is:
//
//      - ( 1 - Z ) <= X <= 1 - Z
//      - ( 1 - Z ) <= Y <= 1 - Z
//                0 <= Z <= 1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    double pyramid_unit_volume: the volume of the unit pyramid.
//
{
  double volume;

  volume = 4.0 / 3.0;

  return volume;
}
//****************************************************************************80

double pyramid_volume ( double h, double s )

//****************************************************************************80
//
//  Purpose:
//
//    pyramid_volume() computes the volume of a pyramid with square base.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    15 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double h, s: the height of the pyramid, and the length of one
//    side of the square base.
//
//  Output:
//
//    double pyramid_volume: the volume of the pyramid.
//
{
  double value;

  value = s * s * h / 3.0;

  return value;
}
//****************************************************************************80

void pyramid_witherden_rule ( int p, int n, double x[], double y[], 
  double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    pyramid_witherden_rule() returns 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.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int p: the precision, 0 <= p <= 10.
//
//    int n: the order of the rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w(n): the quadrature weights.
//
{
  if ( p < 0 ) 
  {
    cout << "\n";
    cout << "pyramid_witherden_rule(): Fatal error!\n";
    cout << "  Input p < 0.\n";
    exit ( 1 );
  }

  if ( 20 < p )
  {
    cout << "\n";
    cout << "pyramid_witherden_rule(): Fatal error!\n";
    cout << "  Input 10 < p.\n";
    exit ( 1 );
  }

  if ( p == 0 )
  {
    rule00 ( n, x, y, z, w );
  }
  else if ( p == 1 )
  {
    rule01 ( n, x, y, z, w );
  }
  else if ( p == 2 )
  {
    rule02 ( n, x, y, z, w );
  }
  else if ( p == 3 )
  {
    rule03 ( n, x, y, z, w );
  }
  else if ( p == 4 )
  {
    rule04 ( n, x, y, z, w );
  }
  else if ( p == 5 )
  {
    rule05 ( n, x, y, z, w );
  }
  else if ( p == 6 )
  {
    rule06 ( n, x, y, z, w );
  }
  else if ( p == 7 )
  {
    rule07 ( n, x, y, z, w );
  }
  else if ( p == 8 )
  {
    rule08 ( n, x, y, z, w );
  }
  else if ( p == 9 )
  {
    rule09 ( n, x, y, z, w );
  }
  else if ( p == 10 )
  {
    rule10 ( n, x, y, z, w );
  }

  return;
}
//****************************************************************************80

double r8_choose ( int n, int k )

//****************************************************************************80
// 
//  Purpose:
//
//    r8_choose() computes the combinatorial coefficient C(N,K).
// 
//  Discussion:
// 
//    Real arithmetic is used, and C(N,K) is computed directly, via
//    Gamma functions, rather than recursively.
// 
//    C(N,K) is the number of distinct combinations of K objects
//    chosen from a set of N distinct objects.  A combination is
//    like a set, in that order does not matter.
// 
//    C(N,K) = N! / ( (N-K)! * K! )
// 
//  Example:
// 
//    The number of combinations of 2 things chosen from 5 is 10.
// 
//    C(5,2) = ( 5 * 4 * 3 * 2 * 1 ) / ( ( 3 * 2 * 1 ) * ( 2 * 1 ) ) = 10.
// 
//    The actual combinations may be represented as:
// 
//      (1,2), (1,3), (1,4), (1,5), (2,3),
//      (2,4), (2,5), (3,4), (3,5), (4,5).
// 
//  Licensing:
// 
//    This code is distributed under the MIT license.
// 
//  Modified:
// 
//    25 July 2011
// 
//  Author:
// 
//    John Burkardt
// 
//  Input:
// 
//    int N, the value of N.
// 
//    int K, the value of K.
//
//  Output:
//
//    double R8_CHOOSE, the value of C(N,K)
// 
{
  double arg;
  double fack;
  double facn;
  double facnmk;
  double value;

  if ( n < 0 )
  {
    value = 0.0;
  }
  else if ( k == 0 )
  {
    value = 1.0;
  }
  else if ( k == 1 )
  {
    value = ( double ) ( n );
  }
  else if ( 1 < k && k < n - 1 )
  {
    arg = ( double ) ( n + 1 );
    facn = lgamma ( arg );

    arg = ( double ) ( k + 1 );
    fack = lgamma ( arg );

    arg = ( double ) ( n - k + 1 );
    facnmk = lgamma ( arg );

    value = round ( exp ( facn - fack - facnmk ) );
  }
  else if ( k == n - 1 )
  {
    value = ( double ) ( n );
  }
  else if ( k == n )
  {
    value = 1.0;
  }
  else
  {
    value = 0.0;
  }

  return value;
}
//****************************************************************************80

void r8vec_copy ( int n, double a1[], double a2[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8vec_copy() copies an R8VEC.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of entries in the vectors.
//
//    double A1[N], the vector to be copied.
//
//  Output:
//
//    double A2[N], the copy of A1.
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return;
}
//****************************************************************************80

int rule_order ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    rule_order() returns the order of a pyramid quadrature rule of given precision.
//
//  Discussion:
//
//    The "order" of a quadrature rule is the number of points involved.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int p: the precision, 0 <= p <= 10.
//
//  Output:
//
//    int rule_order: the order of the rule.
//
{
  int order;
  static int order_save[] = {
      1, 
      1,   5,   6,  10,  15, 
     24,  31,  47,  62,  83 };

  if ( p < 0 ) 
  {
    cout << "\n";
    cout << "rule_order(): Fatal error!\n";
    cout << "  Input p < 0.\n";
    exit ( 1 );
  }

  if ( 10 < p )
  {
    cout << "\n";
    cout << "rule_order(): Fatal error!\n";
    cout << "  Input 10 < p.\n";
    exit ( 1 );
  }

  order = order_save[p];

  return order;
}
void rule00 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule00() returns the rule of precision 0.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000 };

  static double y_save[] = {
          0.0000000000000000 };

  static double z_save[] = {
          0.2500000000000000 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule01 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule01() returns the rule of precision 1.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000 };

  static double y_save[] = {
          0.0000000000000000 };

  static double z_save[] = {
          0.2500000000000000 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule02 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule02() returns the rule of precision 2.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.7189210558117962,
          0.0000000000000000,
         -0.7189210558117962,
          0.0000000000000000 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.7189210558117962,
          0.0000000000000000,
         -0.7189210558117962 };

  static double z_save[] = {
          0.6082910385597788,
          0.1453364835728557,
          0.1453364835728557,
          0.1453364835728557,
          0.1453364835728557 };

  double w_save[] = {
          0.2260773013241023,
          0.1934806746689745,
          0.1934806746689745,
          0.1934806746689745,
          0.1934806746689745 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule03 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule03() returns the rule of precision 3.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.5610836110587396,
          0.5610836110587396,
         -0.5610836110587396,
         -0.5610836110587396 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.5610836110587396,
         -0.5610836110587396,
          0.5610836110587396,
         -0.5610836110587396 };

  static double z_save[] = {
          0.5714285703860683,
          0.0000000056758500,
          0.1666666666666667,
          0.1666666666666667,
          0.1666666666666667,
          0.1666666666666667 };

  double w_save[] = {
          0.2522058839227606,
          0.1125000060660650,
          0.1588235275027936,
          0.1588235275027936,
          0.1588235275027936,
          0.1588235275027936 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule04 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule04() returns the rule of precision 4.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.6505815563982326,
          0.0000000000000000,
         -0.6505815563982326,
          0.0000000000000000,
          0.6579669971216900,
          0.6579669971216900,
         -0.6579669971216900,
         -0.6579669971216900 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.6505815563982326,
          0.0000000000000000,
         -0.6505815563982326,
          0.6579669971216900,
         -0.6579669971216900,
          0.6579669971216900,
         -0.6579669971216900 };

  static double z_save[] = {
          0.6772327888861374,
          0.1251369531087465,
          0.3223841495782137,
          0.3223841495782137,
          0.3223841495782137,
          0.3223841495782137,
          0.0392482838988154,
          0.0392482838988154,
          0.0392482838988154,
          0.0392482838988154 };

  double w_save[] = {
          0.1137418831706419,
          0.2068834025895523,
          0.1063245878893255,
          0.1063245878893255,
          0.1063245878893255,
          0.1063245878893255,
          0.0635190906706259,
          0.0635190906706259,
          0.0635190906706259,
          0.0635190906706259 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule05 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule05() returns the rule of precision 5.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.7065260315463245,
          0.0000000000000000,
         -0.7065260315463245,
          0.0000000000000000,
          0.7051171227788277,
          0.7051171227788277,
         -0.7051171227788277,
         -0.7051171227788277,
          0.4328828641035410,
          0.4328828641035410,
         -0.4328828641035410,
         -0.4328828641035410 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.7065260315463245,
          0.0000000000000000,
         -0.7065260315463245,
          0.7051171227788277,
         -0.7051171227788277,
          0.7051171227788277,
         -0.7051171227788277,
          0.4328828641035410,
         -0.4328828641035410,
          0.4328828641035410,
         -0.4328828641035410 };

  static double z_save[] = {
          0.7298578807825067,
          0.3004010208137690,
          0.0000000064917722,
          0.1250000000000000,
          0.1250000000000000,
          0.1250000000000000,
          0.1250000000000000,
          0.0611119070620230,
          0.0611119070620230,
          0.0611119070620230,
          0.0611119070620230,
          0.4236013371197248,
          0.4236013371197248,
          0.4236013371197248,
          0.4236013371197248 };

  double w_save[] = {
          0.0684353699091401,
          0.1693971144927240,
          0.0587045358285745,
          0.0764412931481202,
          0.0764412931481202,
          0.0764412931481202,
          0.0764412931481202,
          0.0396709015796455,
          0.0396709015796455,
          0.0396709015796455,
          0.0396709015796455,
          0.0597535502146247,
          0.0597535502146247,
          0.0597535502146247,
          0.0597535502146247 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule06 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule06() returns the rule of precision 6.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.8345953511147084,
          0.0000000000000000,
         -0.8345953511147084,
          0.0000000000000000,
          0.4339254093766991,
          0.0000000000000000,
         -0.4339254093766991,
          0.0000000000000000,
          0.5656808544256755,
          0.5656808544256755,
         -0.5656808544256755,
         -0.5656808544256755,
          0.4980790917807059,
          0.4980790917807059,
         -0.4980790917807059,
         -0.4980790917807059,
          0.9508994872144825,
          0.9508994872144825,
         -0.9508994872144825,
         -0.9508994872144825 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.8345953511147084,
          0.0000000000000000,
         -0.8345953511147084,
          0.0000000000000000,
          0.4339254093766991,
          0.0000000000000000,
         -0.4339254093766991,
          0.5656808544256755,
         -0.5656808544256755,
          0.5656808544256755,
         -0.5656808544256755,
          0.4980790917807059,
         -0.4980790917807059,
          0.4980790917807059,
         -0.4980790917807059,
          0.9508994872144825,
         -0.9508994872144825,
          0.9508994872144825,
         -0.9508994872144825 };

  static double z_save[] = {
          0.8076457976939595,
          0.0017638088528196,
          0.1382628064637306,
          0.4214239119356371,
          0.0974473410254620,
          0.0974473410254620,
          0.0974473410254620,
          0.0974473410254620,
          0.5660745906233009,
          0.5660745906233009,
          0.5660745906233009,
          0.5660745906233009,
          0.0294777308457207,
          0.0294777308457207,
          0.0294777308457207,
          0.0294777308457207,
          0.2649158632121295,
          0.2649158632121295,
          0.2649158632121295,
          0.2649158632121295,
          0.0482490706319360,
          0.0482490706319360,
          0.0482490706319360,
          0.0482490706319360 };

  double w_save[] = {
          0.0254628936626420,
          0.0160535131751913,
          0.1195795544525238,
          0.1030606701991518,
          0.0369505879363295,
          0.0369505879363295,
          0.0369505879363295,
          0.0369505879363295,
          0.0315875794881733,
          0.0315875794881733,
          0.0315875794881733,
          0.0315875794881733,
          0.0372001293894483,
          0.0372001293894483,
          0.0372001293894483,
          0.0372001293894483,
          0.0738823846769269,
          0.0738823846769269,
          0.0738823846769269,
          0.0738823846769269,
          0.0043401606367449,
          0.0043401606367449,
          0.0043401606367449,
          0.0043401606367449 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule07 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule07() returns the rule of precision 7.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.8640987597877147,
          0.0000000000000000,
         -0.8640987597877147,
          0.0000000000000000,
          0.6172133998483676,
          0.0000000000000000,
         -0.6172133998483676,
          0.0000000000000000,
          0.3541523808681161,
          0.3541523808681161,
         -0.3541523808681161,
         -0.3541523808681161,
          0.8027258501000878,
          0.8027258501000878,
         -0.8027258501000878,
         -0.8027258501000878,
          0.2541353468618572,
          0.2541353468618572,
         -0.2541353468618572,
         -0.2541353468618572,
          0.6143051077207853,
          0.6143051077207853,
         -0.6143051077207853,
         -0.6143051077207853,
          0.5248326982543755,
          0.5248326982543755,
         -0.5248326982543755,
         -0.5248326982543755 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.8640987597877147,
          0.0000000000000000,
         -0.8640987597877147,
          0.0000000000000000,
          0.6172133998483676,
          0.0000000000000000,
         -0.6172133998483676,
          0.3541523808681161,
         -0.3541523808681161,
          0.3541523808681161,
         -0.3541523808681161,
          0.8027258501000878,
         -0.8027258501000878,
          0.8027258501000878,
         -0.8027258501000878,
          0.2541353468618572,
         -0.2541353468618572,
          0.2541353468618572,
         -0.2541353468618572,
          0.6143051077207853,
         -0.6143051077207853,
          0.6143051077207853,
         -0.6143051077207853,
          0.5248326982543755,
         -0.5248326982543755,
          0.5248326982543755,
         -0.5248326982543755 };

  static double z_save[] = {
          0.0000454802082028,
          0.3935828767576542,
          0.8386827828168515,
          0.0666666666666667,
          0.0666666666666667,
          0.0666666666666667,
          0.0666666666666667,
          0.3333333333333334,
          0.3333333333333334,
          0.3333333333333334,
          0.3333333333333334,
          0.1292864095395495,
          0.1292864095395495,
          0.1292864095395495,
          0.1292864095395495,
          0.0801385301797781,
          0.0801385301797781,
          0.0801385301797781,
          0.0801385301797781,
          0.6055199301105999,
          0.6055199301105999,
          0.6055199301105999,
          0.6055199301105999,
          0.0000000000897583,
          0.0000000000897583,
          0.0000000000897583,
          0.0000000000897583,
          0.2905430754945976,
          0.2905430754945976,
          0.2905430754945976,
          0.2905430754945976 };

  double w_save[] = {
          0.0250005770341970,
          0.1005301826998931,
          0.0157071744270173,
          0.0266917592930029,
          0.0266917592930029,
          0.0266917592930029,
          0.0266917592930029,
          0.0287109375000000,
          0.0287109375000000,
          0.0287109375000000,
          0.0287109375000000,
          0.0659601871499773,
          0.0659601871499773,
          0.0659601871499773,
          0.0659601871499773,
          0.0147861379562357,
          0.0147861379562357,
          0.0147861379562357,
          0.0147861379562357,
          0.0295191096942239,
          0.0295191096942239,
          0.0295191096942239,
          0.0295191096942239,
          0.0133038449696241,
          0.0133038449696241,
          0.0133038449696241,
          0.0133038449696241,
          0.0357185398966593,
          0.0357185398966593,
          0.0357185398966593,
          0.0357185398966593 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule08 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule08() returns the rule of precision 8.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.7960960887742582,
          0.0000000000000000,
         -0.7960960887742582,
          0.0000000000000000,
          0.5648889587418399,
          0.0000000000000000,
         -0.5648889587418399,
          0.0000000000000000,
          0.7530226935489476,
          0.0000000000000000,
         -0.7530226935489476,
          0.0000000000000000,
          0.2465779424622870,
          0.2465779424622870,
         -0.2465779424622870,
         -0.2465779424622870,
          0.7338895325467361,
          0.7338895325467361,
         -0.7338895325467361,
         -0.7338895325467361,
          0.0000826826715128,
          0.0000826826715128,
         -0.0000826826715128,
         -0.0000826826715128,
          0.3989965570952277,
          0.3989965570952277,
         -0.3989965570952277,
         -0.3989965570952277,
          0.3653124144795900,
          0.3653124144795900,
         -0.3653124144795900,
         -0.3653124144795900,
          0.4834252440408382,
          0.4834252440408382,
         -0.4834252440408382,
         -0.4834252440408382,
          0.6314047976005609,
          0.8948561964915306,
          0.6314047976005609,
         -0.8948561964915306,
         -0.6314047976005609,
          0.8948561964915306,
         -0.6314047976005609,
         -0.8948561964915306 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.7960960887742582,
          0.0000000000000000,
         -0.7960960887742582,
          0.0000000000000000,
          0.5648889587418399,
          0.0000000000000000,
         -0.5648889587418399,
          0.0000000000000000,
          0.7530226935489476,
          0.0000000000000000,
         -0.7530226935489476,
          0.2465779424622870,
         -0.2465779424622870,
          0.2465779424622870,
         -0.2465779424622870,
          0.7338895325467361,
         -0.7338895325467361,
          0.7338895325467361,
         -0.7338895325467361,
          0.0000826826715128,
         -0.0000826826715128,
          0.0000826826715128,
         -0.0000826826715128,
          0.3989965570952277,
         -0.3989965570952277,
          0.3989965570952277,
         -0.3989965570952277,
          0.3653124144795900,
         -0.3653124144795900,
          0.3653124144795900,
         -0.3653124144795900,
          0.4834252440408382,
         -0.4834252440408382,
          0.4834252440408382,
         -0.4834252440408382,
          0.8948561964915306,
          0.6314047976005609,
         -0.8948561964915306,
          0.6314047976005609,
          0.8948561964915306,
         -0.6314047976005609,
         -0.8948561964915306,
         -0.6314047976005609 };

  static double z_save[] = {
          0.5848054341081260,
          0.2714503954191895,
          0.0714404302154540,
          0.0000000073236390,
          0.0000000073236390,
          0.0000000073236390,
          0.0000000073236390,
          0.4221837006010586,
          0.4221837006010586,
          0.4221837006010586,
          0.4221837006010586,
          0.1477639972965497,
          0.1477639972965497,
          0.1477639972965497,
          0.1477639972965497,
          0.6734528479254847,
          0.6734528479254847,
          0.6734528479254847,
          0.6734528479254847,
          0.2455708050362506,
          0.2455708050362506,
          0.2455708050362506,
          0.2455708050362506,
          0.8582503626505816,
          0.8582503626505816,
          0.8582503626505816,
          0.8582503626505816,
          0.0474032316194875,
          0.0474032316194875,
          0.0474032316194875,
          0.0474032316194875,
          0.4021857062205061,
          0.4021857062205061,
          0.4021857062205061,
          0.4021857062205061,
          0.1903148210864091,
          0.1903148210864091,
          0.1903148210864091,
          0.1903148210864091,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511,
          0.0454758042237511 };

  double w_save[] = {
          0.0494409951557794,
          0.0873546894411460,
          0.0384605774968570,
          0.0093760594001275,
          0.0093760594001275,
          0.0093760594001275,
          0.0093760594001275,
          0.0163110192164366,
          0.0163110192164366,
          0.0163110192164366,
          0.0163110192164366,
          0.0322397220710709,
          0.0322397220710709,
          0.0322397220710709,
          0.0322397220710709,
          0.0126350996218853,
          0.0126350996218853,
          0.0126350996218853,
          0.0126350996218853,
          0.0062155450261151,
          0.0062155450261151,
          0.0062155450261151,
          0.0062155450261151,
          0.0026132528336493,
          0.0026132528336493,
          0.0026132528336493,
          0.0026132528336493,
          0.0293249681950410,
          0.0293249681950410,
          0.0293249681950410,
          0.0293249681950410,
          0.0363593842641747,
          0.0363593842641747,
          0.0363593842641747,
          0.0363593842641747,
          0.0409820609459212,
          0.0409820609459212,
          0.0409820609459212,
          0.0409820609459212,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665,
          0.0100644114510665 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule09 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule09() returns the rule of precision 9.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.9258200958608426,
          0.0000000000000000,
         -0.9258200958608426,
          0.0000000000000000,
          0.6760969921163639,
          0.0000000000000000,
         -0.6760969921163639,
          0.0000000000000000,
          0.6163321778509252,
          0.0000000000000000,
         -0.6163321778509252,
          0.0000000000000000,
          0.4319711692851803,
          0.0000000000000000,
         -0.4319711692851803,
          0.0000000000000000,
          0.2418018136307699,
          0.2418018136307699,
         -0.2418018136307699,
         -0.2418018136307699,
          0.8440791629031899,
          0.8440791629031899,
         -0.8440791629031899,
         -0.8440791629031899,
          0.2279669387608253,
          0.2279669387608253,
         -0.2279669387608253,
         -0.2279669387608253,
          0.5027849771350101,
          0.5027849771350101,
         -0.5027849771350101,
         -0.5027849771350101,
          0.2605010749834311,
          0.2605010749834311,
         -0.2605010749834311,
         -0.2605010749834311,
          0.0926958730867243,
          0.0926958730867243,
         -0.0926958730867243,
         -0.0926958730867243,
          0.4832161680706445,
          0.4832161680706445,
         -0.4832161680706445,
         -0.4832161680706445,
          0.5671855056082324,
          0.5671855056082324,
         -0.5671855056082324,
         -0.5671855056082324,
          0.8463949915138761,
          0.8463949915138761,
         -0.8463949915138761,
         -0.8463949915138761,
          0.8471778681177453,
          0.4641907129661463,
          0.8471778681177453,
         -0.4641907129661463,
         -0.8471778681177453,
          0.4641907129661463,
         -0.8471778681177453,
         -0.4641907129661463 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.9258200958608426,
          0.0000000000000000,
         -0.9258200958608426,
          0.0000000000000000,
          0.6760969921163639,
          0.0000000000000000,
         -0.6760969921163639,
          0.0000000000000000,
          0.6163321778509252,
          0.0000000000000000,
         -0.6163321778509252,
          0.0000000000000000,
          0.4319711692851803,
          0.0000000000000000,
         -0.4319711692851803,
          0.2418018136307699,
         -0.2418018136307699,
          0.2418018136307699,
         -0.2418018136307699,
          0.8440791629031899,
         -0.8440791629031899,
          0.8440791629031899,
         -0.8440791629031899,
          0.2279669387608253,
         -0.2279669387608253,
          0.2279669387608253,
         -0.2279669387608253,
          0.5027849771350101,
         -0.5027849771350101,
          0.5027849771350101,
         -0.5027849771350101,
          0.2605010749834311,
         -0.2605010749834311,
          0.2605010749834311,
         -0.2605010749834311,
          0.0926958730867243,
         -0.0926958730867243,
          0.0926958730867243,
         -0.0926958730867243,
          0.4832161680706445,
         -0.4832161680706445,
          0.4832161680706445,
         -0.4832161680706445,
          0.5671855056082324,
         -0.5671855056082324,
          0.5671855056082324,
         -0.5671855056082324,
          0.8463949915138761,
         -0.8463949915138761,
          0.8463949915138761,
         -0.8463949915138761,
          0.4641907129661463,
          0.8471778681177453,
         -0.4641907129661463,
          0.8471778681177453,
          0.4641907129661463,
         -0.8471778681177453,
         -0.4641907129661463,
         -0.8471778681177453 };

  static double z_save[] = {
          0.0371459309315055,
          0.6942002631703261,
          0.0000000042251285,
          0.0000000042251285,
          0.0000000042251285,
          0.0000000042251285,
          0.2697317845200571,
          0.2697317845200571,
          0.2697317845200571,
          0.2697317845200571,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.5334178104457834,
          0.5334178104457834,
          0.5334178104457834,
          0.5334178104457834,
          0.6619880087153507,
          0.6619880087153507,
          0.6619880087153507,
          0.6619880087153507,
          0.1522047098054636,
          0.1522047098054636,
          0.1522047098054636,
          0.1522047098054636,
          0.4222655365894455,
          0.4222655365894455,
          0.4222655365894455,
          0.4222655365894455,
          0.3967643698988902,
          0.3967643698988902,
          0.3967643698988902,
          0.3967643698988902,
          0.2008207219371707,
          0.2008207219371707,
          0.2008207219371707,
          0.2008207219371707,
          0.8661453707796378,
          0.8661453707796378,
          0.8661453707796378,
          0.8661453707796378,
          0.0195040124924115,
          0.0195040124924115,
          0.0195040124924115,
          0.0195040124924115,
          0.2131686017065601,
          0.2131686017065601,
          0.2131686017065601,
          0.2131686017065601,
          0.0153477972480097,
          0.0153477972480097,
          0.0153477972480097,
          0.0153477972480097,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333,
          0.0833333333333333 };

  double w_save[] = {
          0.0320478691353425,
          0.0237932209450781,
          0.0041733723029944,
          0.0041733723029944,
          0.0041733723029944,
          0.0041733723029944,
          0.0222937913596180,
          0.0222937913596180,
          0.0222937913596180,
          0.0222937913596180,
          0.0291280227779532,
          0.0291280227779532,
          0.0291280227779532,
          0.0291280227779532,
          0.0115418065483010,
          0.0115418065483010,
          0.0115418065483010,
          0.0115418065483010,
          0.0101005748273597,
          0.0101005748273597,
          0.0101005748273597,
          0.0101005748273597,
          0.0020589782014206,
          0.0020589782014206,
          0.0020589782014206,
          0.0020589782014206,
          0.0342456887657178,
          0.0342456887657178,
          0.0342456887657178,
          0.0342456887657178,
          0.0129094136418972,
          0.0129094136418972,
          0.0129094136418972,
          0.0129094136418972,
          0.0371122811478335,
          0.0371122811478335,
          0.0371122811478335,
          0.0371122811478335,
          0.0022186981787511,
          0.0022186981787511,
          0.0022186981787511,
          0.0022186981787511,
          0.0159562091361951,
          0.0159562091361951,
          0.0159562091361951,
          0.0159562091361951,
          0.0230673181273276,
          0.0230673181273276,
          0.0230673181273276,
          0.0230673181273276,
          0.0045994005726981,
          0.0045994005726981,
          0.0045994005726981,
          0.0045994005726981,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138,
          0.0133170859459138 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule10 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule10() returns the rule of precision 10.
//
//  Discussion:
//
//    We suppose we are given a pyramid P with vertices
//      (-1,-1,0), (-1,+1,0), (+1,+1,0), (+1,-1,0), (0,0,1).
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    P is approximated by Q as follows:
//
//    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.9211301560146403,
          0.0000000000000000,
         -0.9211301560146403,
          0.0000000000000000,
          0.4546948844252537,
          0.0000000000000000,
         -0.4546948844252537,
          0.0000000000000000,
          0.3495469844967067,
          0.0000000000000000,
         -0.3495469844967067,
          0.0000000000000000,
          0.4036909359989291,
          0.4036909359989291,
         -0.4036909359989291,
         -0.4036909359989291,
          0.3714608909064742,
          0.3714608909064742,
         -0.3714608909064742,
         -0.3714608909064742,
          0.2995979966434122,
          0.2995979966434122,
         -0.2995979966434122,
         -0.2995979966434122,
          0.1582262545541785,
          0.1582262545541785,
         -0.1582262545541785,
         -0.1582262545541785,
          0.1819979349519966,
          0.1819979349519966,
         -0.1819979349519966,
         -0.1819979349519966,
          0.8455841498012413,
          0.8455841498012413,
         -0.8455841498012413,
         -0.8455841498012413,
          0.6468694059373429,
          0.6468694059373429,
         -0.6468694059373429,
         -0.6468694059373429,
          0.5875870294087208,
          0.5875870294087208,
         -0.5875870294087208,
         -0.5875870294087208,
          0.2237620599838161,
          0.2237620599838161,
         -0.2237620599838161,
         -0.2237620599838161,
          0.7455657084697210,
          0.2247174318849062,
          0.7455657084697210,
         -0.2247174318849062,
         -0.7455657084697210,
          0.2247174318849062,
         -0.7455657084697210,
         -0.2247174318849062,
          0.4056619527318771,
          0.7487267628018788,
          0.4056619527318771,
         -0.7487267628018788,
         -0.4056619527318771,
          0.7487267628018788,
         -0.4056619527318771,
         -0.7487267628018788,
          0.0707425682812096,
          0.6004569111268736,
          0.0707425682812096,
         -0.6004569111268736,
         -0.0707425682812096,
          0.6004569111268736,
         -0.0707425682812096,
         -0.6004569111268736,
          0.9422191470796681,
          0.6375631040785387,
          0.9422191470796681,
         -0.6375631040785387,
         -0.9422191470796681,
          0.6375631040785387,
         -0.9422191470796681,
         -0.6375631040785387 };

  static double y_save[] = {
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.0000000000000000,
          0.9211301560146403,
          0.0000000000000000,
         -0.9211301560146403,
          0.0000000000000000,
          0.4546948844252537,
          0.0000000000000000,
         -0.4546948844252537,
          0.0000000000000000,
          0.3495469844967067,
          0.0000000000000000,
         -0.3495469844967067,
          0.4036909359989291,
         -0.4036909359989291,
          0.4036909359989291,
         -0.4036909359989291,
          0.3714608909064742,
         -0.3714608909064742,
          0.3714608909064742,
         -0.3714608909064742,
          0.2995979966434122,
         -0.2995979966434122,
          0.2995979966434122,
         -0.2995979966434122,
          0.1582262545541785,
         -0.1582262545541785,
          0.1582262545541785,
         -0.1582262545541785,
          0.1819979349519966,
         -0.1819979349519966,
          0.1819979349519966,
         -0.1819979349519966,
          0.8455841498012413,
         -0.8455841498012413,
          0.8455841498012413,
         -0.8455841498012413,
          0.6468694059373429,
         -0.6468694059373429,
          0.6468694059373429,
         -0.6468694059373429,
          0.5875870294087208,
         -0.5875870294087208,
          0.5875870294087208,
         -0.5875870294087208,
          0.2237620599838161,
         -0.2237620599838161,
          0.2237620599838161,
         -0.2237620599838161,
          0.2247174318849062,
          0.7455657084697210,
         -0.2247174318849062,
          0.7455657084697210,
          0.2247174318849062,
         -0.7455657084697210,
         -0.2247174318849062,
         -0.7455657084697210,
          0.7487267628018788,
          0.4056619527318771,
         -0.7487267628018788,
          0.4056619527318771,
          0.7487267628018788,
         -0.4056619527318771,
         -0.7487267628018788,
         -0.4056619527318771,
          0.6004569111268736,
          0.0707425682812096,
         -0.6004569111268736,
          0.0707425682812096,
          0.6004569111268736,
         -0.0707425682812096,
         -0.6004569111268736,
         -0.0707425682812096,
          0.6375631040785387,
          0.9422191470796681,
         -0.6375631040785387,
          0.9422191470796681,
          0.6375631040785387,
         -0.9422191470796681,
         -0.6375631040785387,
         -0.9422191470796681 };

  static double z_save[] = {
          0.9167957817791272,
          0.3511368063403526,
          0.7419135679633000,
          0.0631615972799145,
          0.0631615972799145,
          0.0631615972799145,
          0.0631615972799145,
          0.1767304111617324,
          0.1767304111617324,
          0.1767304111617324,
          0.1767304111617324,
          0.6020938107079140,
          0.6020938107079140,
          0.6020938107079140,
          0.6020938107079140,
          0.5270494697681531,
          0.5270494697681531,
          0.5270494697681531,
          0.5270494697681531,
          0.3491104418985491,
          0.3491104418985491,
          0.3491104418985491,
          0.3491104418985491,
          0.0069339728754634,
          0.0069339728754634,
          0.0069339728754634,
          0.0069339728754634,
          0.7752643179331966,
          0.7752643179331966,
          0.7752643179331966,
          0.7752643179331966,
          0.5470650181463391,
          0.5470650181463391,
          0.5470650181463391,
          0.5470650181463391,
          0.0728984114756605,
          0.0728984114756605,
          0.0728984114756605,
          0.0728984114756605,
          0.2836223397977548,
          0.2836223397977548,
          0.2836223397977548,
          0.2836223397977548,
          0.0598130217927265,
          0.0598130217927265,
          0.0598130217927265,
          0.0598130217927265,
          0.0658308799806233,
          0.0658308799806233,
          0.0658308799806233,
          0.0658308799806233,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.0272502552746740,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.1618000129270508,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.3584705635890507,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847,
          0.0074406410768847 };

  double w_save[] = {
          0.0024012267625214,
          0.0469159279097637,
          0.0098466607485292,
          0.0056773923179313,
          0.0056773923179313,
          0.0056773923179313,
          0.0056773923179313,
          0.0412212597418143,
          0.0412212597418143,
          0.0412212597418143,
          0.0412212597418143,
          0.0092087532680317,
          0.0092087532680317,
          0.0092087532680317,
          0.0092087532680317,
          0.0070683961650691,
          0.0070683961650691,
          0.0070683961650691,
          0.0070683961650691,
          0.0298306858263750,
          0.0298306858263750,
          0.0298306858263750,
          0.0298306858263750,
          0.0058699928018168,
          0.0058699928018168,
          0.0058699928018168,
          0.0058699928018168,
          0.0049543654347689,
          0.0049543654347689,
          0.0049543654347689,
          0.0049543654347689,
          0.0159309372716684,
          0.0159309372716684,
          0.0159309372716684,
          0.0159309372716684,
          0.0050915274851829,
          0.0050915274851829,
          0.0050915274851829,
          0.0050915274851829,
          0.0068457933453088,
          0.0068457933453088,
          0.0068457933453088,
          0.0068457933453088,
          0.0150380203803664,
          0.0150380203803664,
          0.0150380203803664,
          0.0150380203803664,
          0.0156151631237942,
          0.0156151631237942,
          0.0156151631237942,
          0.0156151631237942,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0089093600668980,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0167917214598521,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0082057491053038,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804,
          0.0025215488592804 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}