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

using namespace std;

# include "tetrahedron_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;
}
//****************************************************************************80

double r8mat_det_4d ( double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_det_4d() computes the determinant of a 4 by 4 R8MAT.
//
//  Discussion:
//
//    Although in this special case it doesn't matter so much, the
//    matrix is now stored in row-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double a[4*4], the matrix whose determinant is desired.
//
//  Output:
//
//    double r8mat_det_4d: the determinant of the matrix.
//
{
  double det;

  det =
      a[0*4+0] * (
          a[1*4+1] * ( a[2*4+2] * a[3*4+3] - a[2*4+3] * a[3*4+2] )
        - a[1*4+2] * ( a[2*4+1] * a[3*4+3] - a[2*4+3] * a[3*4+1] )
        + a[1*4+3] * ( a[2*4+1] * a[3*4+2] - a[2*4+2] * a[3*4+1] ) )
    - a[0*4+1] * (
          a[1*4+0] * ( a[2*4+2] * a[3*4+3] - a[2*4+3] * a[3*4+2] )
        - a[1*4+2] * ( a[2*4+0] * a[3*4+3] - a[2*4+3] * a[3*4+0] )
        + a[1*4+3] * ( a[2*4+0] * a[3*4+2] - a[2*4+2] * a[3*4+0] ) )
    + a[0*4+2] * (
          a[1*4+0] * ( a[2*4+1] * a[3*4+3] - a[2*4+3] * a[3*4+1] )
        - a[1*4+1] * ( a[2*4+0] * a[3*4+3] - a[2*4+3] * a[3*4+0] )
        + a[1*4+3] * ( a[2*4+0] * a[3*4+1] - a[2*4+1] * a[3*4+0] ) )
    - a[0*4+3] * (
          a[1*4+0] * ( a[2*4+1] * a[3*4+2] - a[2*4+2] * a[3*4+1] )
        - a[1*4+1] * ( a[2*4+0] * a[3*4+2] - a[2*4+2] * a[3*4+0] )
        + a[1*4+2] * ( a[2*4+0] * a[3*4+1] - a[2*4+1] * a[3*4+0] ) );

  return det;
}
//****************************************************************************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 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:
//
//    25 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Jan Jaskowiec, Natarajan Sukumar,
//    High order symmetric cubature rules for tetrahedra and pyramids,
//    International Journal of Numerical Methods in Engineering,
//    Volume 122, Number 1, pages 148-171, 24 August 2020.
//
//  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,   4,   8,  14,  14, 
     24,  35,  46,  59,  81 };

  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 a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule00() returns the rule of precision 0.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.2500000000000000 };

  static double b_save[] = {
          0.2500000000000000 };

  static double c_save[] = {
          0.2500000000000000 };

  static double d_save[] = {
          0.2500000000000000 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule01 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule01() returns the rule of precision 1.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.2500000000000000 };

  static double b_save[] = {
          0.2500000000000000 };

  static double c_save[] = {
          0.2500000000000000 };

  static double d_save[] = {
          0.2500000000000000 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule02 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule02() returns the rule of precision 2.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.1381966011250105,
          0.1381966011250105,
          0.5854101966249685,
          0.1381966011250105 };

  static double b_save[] = {
          0.1381966011250105,
          0.5854101966249685,
          0.1381966011250105,
          0.1381966011250105 };

  static double c_save[] = {
          0.5854101966249685,
          0.1381966011250105,
          0.1381966011250105,
          0.1381966011250105 };

  static double d_save[] = {
          0.1381966011250104,
          0.1381966011250104,
          0.1381966011250104,
          0.5854101966249684 };

  double w_save[] = {
          0.2500000000000000,
          0.2500000000000000,
          0.2500000000000000,
          0.2500000000000000 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule03 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule03() returns the rule of precision 3.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3281633025163817,
          0.3281633025163817,
          0.0155100924508549,
          0.3281633025163817,
          0.1080472498984286,
          0.1080472498984286,
          0.6758582503047141,
          0.1080472498984286 };

  static double b_save[] = {
          0.3281633025163817,
          0.0155100924508549,
          0.3281633025163817,
          0.3281633025163817,
          0.1080472498984286,
          0.6758582503047141,
          0.1080472498984286,
          0.1080472498984286 };

  static double c_save[] = {
          0.0155100924508549,
          0.3281633025163817,
          0.3281633025163817,
          0.3281633025163817,
          0.6758582503047141,
          0.1080472498984286,
          0.1080472498984286,
          0.1080472498984286 };

  static double d_save[] = {
          0.3281633025163816,
          0.3281633025163817,
          0.3281633025163817,
          0.0155100924508549,
          0.1080472498984286,
          0.1080472498984286,
          0.1080472498984286,
          0.6758582503047141 };

  double w_save[] = {
          0.1362178425370874,
          0.1362178425370874,
          0.1362178425370874,
          0.1362178425370874,
          0.1137821574629126,
          0.1137821574629126,
          0.1137821574629126,
          0.1137821574629126 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule04 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule04() returns the rule of precision 4.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3108859192633006,
          0.3108859192633006,
          0.0673422422100982,
          0.3108859192633006,
          0.0927352503108912,
          0.0927352503108912,
          0.7217942490673264,
          0.0927352503108912,
          0.0455037041256496,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.4544962958743504 };

  static double b_save[] = {
          0.3108859192633006,
          0.0673422422100982,
          0.3108859192633006,
          0.3108859192633006,
          0.0927352503108912,
          0.7217942490673264,
          0.0927352503108912,
          0.0927352503108912,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.0455037041256496,
          0.4544962958743504 };

  static double c_save[] = {
          0.0673422422100982,
          0.3108859192633006,
          0.3108859192633006,
          0.3108859192633006,
          0.7217942490673264,
          0.0927352503108912,
          0.0927352503108912,
          0.0927352503108912,
          0.4544962958743504,
          0.4544962958743504,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.0455037041256496 };

  static double d_save[] = {
          0.3108859192633006,
          0.3108859192633005,
          0.3108859192633005,
          0.0673422422100981,
          0.0927352503108911,
          0.0927352503108911,
          0.0927352503108911,
          0.7217942490673263,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.4544962958743504,
          0.4544962958743504,
          0.0455037041256496 };

  double w_save[] = {
          0.1126879257180159,
          0.1126879257180159,
          0.1126879257180159,
          0.1126879257180159,
          0.0734930431163620,
          0.0734930431163620,
          0.0734930431163620,
          0.0734930431163620,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule05 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule05() returns the rule of precision 5.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3108859192633006,
          0.3108859192633006,
          0.0673422422100982,
          0.3108859192633006,
          0.0927352503108912,
          0.0927352503108912,
          0.7217942490673264,
          0.0927352503108912,
          0.0455037041256496,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.4544962958743504 };

  static double b_save[] = {
          0.3108859192633006,
          0.0673422422100982,
          0.3108859192633006,
          0.3108859192633006,
          0.0927352503108912,
          0.7217942490673264,
          0.0927352503108912,
          0.0927352503108912,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.0455037041256496,
          0.4544962958743504 };

  static double c_save[] = {
          0.0673422422100982,
          0.3108859192633006,
          0.3108859192633006,
          0.3108859192633006,
          0.7217942490673264,
          0.0927352503108912,
          0.0927352503108912,
          0.0927352503108912,
          0.4544962958743504,
          0.4544962958743504,
          0.4544962958743504,
          0.0455037041256496,
          0.0455037041256496,
          0.0455037041256496 };

  static double d_save[] = {
          0.3108859192633006,
          0.3108859192633005,
          0.3108859192633005,
          0.0673422422100981,
          0.0927352503108911,
          0.0927352503108911,
          0.0927352503108911,
          0.7217942490673263,
          0.0455037041256496,
          0.0455037041256496,
          0.4544962958743504,
          0.4544962958743504,
          0.4544962958743504,
          0.0455037041256496 };

  double w_save[] = {
          0.1126879257180159,
          0.1126879257180159,
          0.1126879257180159,
          0.1126879257180159,
          0.0734930431163620,
          0.0734930431163620,
          0.0734930431163620,
          0.0734930431163620,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815,
          0.0425460207770815 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule06 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule06() returns the rule of precision 6.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.0406739585346114,
          0.0406739585346114,
          0.8779781243961660,
          0.0406739585346114,
          0.3223378901422755,
          0.3223378901422755,
          0.0329863295731735,
          0.3223378901422755,
          0.2146028712591520,
          0.2146028712591520,
          0.3561913862225439,
          0.2146028712591520,
          0.6030056647916492,
          0.6030056647916492,
          0.0636610018750175,
          0.2696723314583158,
          0.0636610018750175,
          0.0636610018750175,
          0.2696723314583158,
          0.0636610018750175,
          0.0636610018750175,
          0.0636610018750175,
          0.2696723314583158,
          0.6030056647916492 };

  static double b_save[] = {
          0.0406739585346114,
          0.8779781243961660,
          0.0406739585346114,
          0.0406739585346114,
          0.3223378901422755,
          0.0329863295731735,
          0.3223378901422755,
          0.3223378901422755,
          0.2146028712591520,
          0.3561913862225439,
          0.2146028712591520,
          0.2146028712591520,
          0.0636610018750175,
          0.0636610018750175,
          0.0636610018750175,
          0.6030056647916492,
          0.2696723314583158,
          0.6030056647916492,
          0.0636610018750175,
          0.2696723314583158,
          0.0636610018750175,
          0.6030056647916492,
          0.0636610018750175,
          0.2696723314583158 };

  static double c_save[] = {
          0.8779781243961660,
          0.0406739585346114,
          0.0406739585346114,
          0.0406739585346114,
          0.0329863295731735,
          0.3223378901422755,
          0.3223378901422755,
          0.3223378901422755,
          0.3561913862225439,
          0.2146028712591520,
          0.2146028712591520,
          0.2146028712591520,
          0.2696723314583158,
          0.0636610018750175,
          0.6030056647916492,
          0.0636610018750175,
          0.6030056647916492,
          0.0636610018750175,
          0.6030056647916492,
          0.0636610018750175,
          0.2696723314583158,
          0.2696723314583158,
          0.0636610018750175,
          0.0636610018750175 };

  static double d_save[] = {
          0.0406739585346112,
          0.0406739585346113,
          0.0406739585346113,
          0.8779781243961657,
          0.3223378901422754,
          0.3223378901422754,
          0.3223378901422754,
          0.0329863295731734,
          0.2146028712591521,
          0.2146028712591521,
          0.2146028712591521,
          0.3561913862225440,
          0.0636610018750175,
          0.2696723314583158,
          0.2696723314583159,
          0.0636610018750175,
          0.0636610018750176,
          0.2696723314583159,
          0.0636610018750176,
          0.6030056647916493,
          0.6030056647916493,
          0.0636610018750176,
          0.6030056647916493,
          0.0636610018750175 };

  double w_save[] = {
          0.0100772110553206,
          0.0100772110553206,
          0.0100772110553206,
          0.0100772110553206,
          0.0553571815436547,
          0.0553571815436547,
          0.0553571815436547,
          0.0553571815436547,
          0.0399227502581675,
          0.0399227502581675,
          0.0399227502581675,
          0.0399227502581675,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857,
          0.0482142857142857 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule07 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule07() returns the rule of precision 7.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.2500000000000000,
          0.3157011497782028,
          0.3157011497782028,
          0.0528965506653916,
          0.3157011497782028,
          0.0504898225983964,
          0.4495101774016036,
          0.0504898225983964,
          0.0504898225983964,
          0.4495101774016036,
          0.4495101774016036,
          0.5751716375870000,
          0.5751716375870000,
          0.1888338310260010,
          0.0471607003609979,
          0.1888338310260010,
          0.1888338310260010,
          0.0471607003609979,
          0.1888338310260010,
          0.1888338310260010,
          0.1888338310260010,
          0.0471607003609979,
          0.5751716375870000,
          0.8108302410985486,
          0.8108302410985486,
          0.0212654725414833,
          0.1466388138184849,
          0.0212654725414833,
          0.0212654725414833,
          0.1466388138184849,
          0.0212654725414833,
          0.0212654725414833,
          0.0212654725414833,
          0.1466388138184849,
          0.8108302410985486 };

  static double b_save[] = {
          0.2500000000000000,
          0.3157011497782028,
          0.0528965506653916,
          0.3157011497782028,
          0.3157011497782028,
          0.4495101774016036,
          0.0504898225983964,
          0.0504898225983964,
          0.4495101774016036,
          0.0504898225983964,
          0.4495101774016036,
          0.1888338310260010,
          0.1888338310260010,
          0.1888338310260010,
          0.5751716375870000,
          0.0471607003609979,
          0.5751716375870000,
          0.1888338310260010,
          0.0471607003609979,
          0.1888338310260010,
          0.5751716375870000,
          0.1888338310260010,
          0.0471607003609979,
          0.0212654725414833,
          0.0212654725414833,
          0.0212654725414833,
          0.8108302410985486,
          0.1466388138184849,
          0.8108302410985486,
          0.0212654725414833,
          0.1466388138184849,
          0.0212654725414833,
          0.8108302410985486,
          0.0212654725414833,
          0.1466388138184849 };

  static double c_save[] = {
          0.2500000000000000,
          0.0528965506653916,
          0.3157011497782028,
          0.3157011497782028,
          0.3157011497782028,
          0.4495101774016036,
          0.4495101774016036,
          0.4495101774016036,
          0.0504898225983964,
          0.0504898225983964,
          0.0504898225983964,
          0.0471607003609979,
          0.1888338310260010,
          0.5751716375870000,
          0.1888338310260010,
          0.5751716375870000,
          0.1888338310260010,
          0.5751716375870000,
          0.1888338310260010,
          0.0471607003609979,
          0.0471607003609979,
          0.1888338310260010,
          0.1888338310260010,
          0.1466388138184849,
          0.0212654725414833,
          0.8108302410985486,
          0.0212654725414833,
          0.8108302410985486,
          0.0212654725414833,
          0.8108302410985486,
          0.0212654725414833,
          0.1466388138184849,
          0.1466388138184849,
          0.0212654725414833,
          0.0212654725414833 };

  static double d_save[] = {
          0.2500000000000000,
          0.3157011497782028,
          0.3157011497782029,
          0.3157011497782029,
          0.0528965506653917,
          0.0504898225983964,
          0.0504898225983963,
          0.4495101774016038,
          0.4495101774016036,
          0.4495101774016036,
          0.0504898225983963,
          0.1888338310260011,
          0.0471607003609980,
          0.0471607003609980,
          0.1888338310260010,
          0.1888338310260012,
          0.0471607003609980,
          0.1888338310260010,
          0.5751716375870001,
          0.5751716375870000,
          0.1888338310260011,
          0.5751716375870000,
          0.1888338310260011,
          0.0212654725414832,
          0.1466388138184849,
          0.1466388138184849,
          0.0212654725414833,
          0.0212654725414831,
          0.1466388138184849,
          0.0212654725414833,
          0.8108302410985485,
          0.8108302410985486,
          0.0212654725414832,
          0.8108302410985486,
          0.0212654725414832 };

  double w_save[] = {
          0.0954852894641308,
          0.0423295812099670,
          0.0423295812099670,
          0.0423295812099670,
          0.0423295812099670,
          0.0318969278328576,
          0.0318969278328576,
          0.0318969278328576,
          0.0318969278328576,
          0.0318969278328576,
          0.0318969278328576,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0372071307283346,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033,
          0.0081107708299033 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule08 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule08() returns the rule of precision 8.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.1079527249622109,
          0.1079527249622109,
          0.6761418251133675,
          0.1079527249622109,
          0.1851094877825866,
          0.1851094877825866,
          0.4446715366522403,
          0.1851094877825866,
          0.0423165436847673,
          0.0423165436847673,
          0.8730503689456981,
          0.0423165436847673,
          0.3141817091240390,
          0.3141817091240390,
          0.0574548726278830,
          0.3141817091240390,
          0.4355913285838302,
          0.0644086714161698,
          0.4355913285838302,
          0.4355913285838302,
          0.0644086714161698,
          0.0644086714161698,
          0.7174640634263083,
          0.7174640634263083,
          0.0214339301271306,
          0.2396680763194305,
          0.0214339301271306,
          0.0214339301271306,
          0.2396680763194305,
          0.0214339301271306,
          0.0214339301271306,
          0.0214339301271306,
          0.2396680763194305,
          0.7174640634263083,
          0.5837973783021444,
          0.5837973783021444,
          0.2041393338760291,
          0.0079239539457974,
          0.2041393338760291,
          0.2041393338760291,
          0.0079239539457974,
          0.2041393338760291,
          0.2041393338760291,
          0.2041393338760291,
          0.0079239539457974,
          0.5837973783021444 };

  static double b_save[] = {
          0.1079527249622109,
          0.6761418251133675,
          0.1079527249622109,
          0.1079527249622109,
          0.1851094877825866,
          0.4446715366522403,
          0.1851094877825866,
          0.1851094877825866,
          0.0423165436847673,
          0.8730503689456981,
          0.0423165436847673,
          0.0423165436847673,
          0.3141817091240390,
          0.0574548726278830,
          0.3141817091240390,
          0.3141817091240390,
          0.0644086714161698,
          0.4355913285838302,
          0.4355913285838302,
          0.0644086714161698,
          0.4355913285838302,
          0.0644086714161698,
          0.0214339301271306,
          0.0214339301271306,
          0.0214339301271306,
          0.7174640634263083,
          0.2396680763194305,
          0.7174640634263083,
          0.0214339301271306,
          0.2396680763194305,
          0.0214339301271306,
          0.7174640634263083,
          0.0214339301271306,
          0.2396680763194305,
          0.2041393338760291,
          0.2041393338760291,
          0.2041393338760291,
          0.5837973783021444,
          0.0079239539457974,
          0.5837973783021444,
          0.2041393338760291,
          0.0079239539457974,
          0.2041393338760291,
          0.5837973783021444,
          0.2041393338760291,
          0.0079239539457974 };

  static double c_save[] = {
          0.6761418251133675,
          0.1079527249622109,
          0.1079527249622109,
          0.1079527249622109,
          0.4446715366522403,
          0.1851094877825866,
          0.1851094877825866,
          0.1851094877825866,
          0.8730503689456981,
          0.0423165436847673,
          0.0423165436847673,
          0.0423165436847673,
          0.0574548726278830,
          0.3141817091240390,
          0.3141817091240390,
          0.3141817091240390,
          0.0644086714161698,
          0.0644086714161698,
          0.0644086714161698,
          0.4355913285838302,
          0.4355913285838302,
          0.4355913285838302,
          0.2396680763194305,
          0.0214339301271306,
          0.7174640634263083,
          0.0214339301271306,
          0.7174640634263083,
          0.0214339301271306,
          0.7174640634263083,
          0.0214339301271306,
          0.2396680763194305,
          0.2396680763194305,
          0.0214339301271306,
          0.0214339301271306,
          0.0079239539457974,
          0.2041393338760291,
          0.5837973783021444,
          0.2041393338760291,
          0.5837973783021444,
          0.2041393338760291,
          0.5837973783021444,
          0.2041393338760291,
          0.0079239539457974,
          0.0079239539457974,
          0.2041393338760291,
          0.2041393338760291 };

  static double d_save[] = {
          0.1079527249622108,
          0.1079527249622108,
          0.1079527249622108,
          0.6761418251133674,
          0.1851094877825867,
          0.1851094877825866,
          0.1851094877825866,
          0.4446715366522404,
          0.0423165436847674,
          0.0423165436847674,
          0.0423165436847673,
          0.8730503689456983,
          0.3141817091240390,
          0.3141817091240390,
          0.3141817091240391,
          0.0574548726278831,
          0.4355913285838302,
          0.4355913285838302,
          0.0644086714161698,
          0.0644086714161698,
          0.0644086714161698,
          0.4355913285838302,
          0.0214339301271306,
          0.2396680763194305,
          0.2396680763194305,
          0.0214339301271306,
          0.0214339301271306,
          0.2396680763194305,
          0.0214339301271306,
          0.7174640634263083,
          0.7174640634263083,
          0.0214339301271306,
          0.7174640634263083,
          0.0214339301271306,
          0.2041393338760291,
          0.0079239539457974,
          0.0079239539457974,
          0.2041393338760291,
          0.2041393338760291,
          0.0079239539457974,
          0.2041393338760291,
          0.5837973783021444,
          0.5837973783021444,
          0.2041393338760291,
          0.5837973783021444,
          0.2041393338760291 };

  double w_save[] = {
          0.0264266509084088,
          0.0264266509084088,
          0.0264266509084088,
          0.0264266509084088,
          0.0520317475637385,
          0.0520317475637385,
          0.0520317475637385,
          0.0520317475637385,
          0.0075252561535402,
          0.0075252561535402,
          0.0075252561535402,
          0.0075252561535402,
          0.0417637828569349,
          0.0417637828569349,
          0.0417637828569349,
          0.0417637828569349,
          0.0362809302613088,
          0.0362809302613088,
          0.0362809302613088,
          0.0362809302613088,
          0.0362809302613088,
          0.0362809302613088,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0071569028908444,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603,
          0.0154534861509603 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule09 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule09() returns the rule of precision 9.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.2500000000000000,
          0.0000000006198170,
          0.0000000006198170,
          0.9999999981405490,
          0.0000000006198170,
          0.1607745353952616,
          0.1607745353952616,
          0.5176763938142153,
          0.1607745353952616,
          0.3222765218214210,
          0.3222765218214210,
          0.0331704345357371,
          0.3222765218214210,
          0.0451089183454136,
          0.0451089183454136,
          0.8646732449637593,
          0.0451089183454136,
          0.1122965460043761,
          0.3877034539956239,
          0.1122965460043761,
          0.1122965460043761,
          0.3877034539956239,
          0.3877034539956239,
          0.0025545792330413,
          0.0025545792330413,
          0.4588714487524593,
          0.0797025232620401,
          0.4588714487524593,
          0.4588714487524593,
          0.0797025232620401,
          0.4588714487524593,
          0.4588714487524593,
          0.4588714487524593,
          0.0797025232620401,
          0.0025545792330413,
          0.7183503264420745,
          0.7183503264420745,
          0.0337758706853386,
          0.2140979321872483,
          0.0337758706853386,
          0.0337758706853386,
          0.2140979321872483,
          0.0337758706853386,
          0.0337758706853386,
          0.0337758706853386,
          0.2140979321872483,
          0.7183503264420745,
          0.0344159105781753,
          0.0344159105781753,
          0.1836413698099279,
          0.5983013498019689,
          0.1836413698099279,
          0.1836413698099279,
          0.5983013498019689,
          0.1836413698099279,
          0.1836413698099279,
          0.1836413698099279,
          0.5983013498019689,
          0.0344159105781753 };

  static double b_save[] = {
          0.2500000000000000,
          0.0000000006198170,
          0.9999999981405490,
          0.0000000006198170,
          0.0000000006198170,
          0.1607745353952616,
          0.5176763938142153,
          0.1607745353952616,
          0.1607745353952616,
          0.3222765218214210,
          0.0331704345357371,
          0.3222765218214210,
          0.3222765218214210,
          0.0451089183454136,
          0.8646732449637593,
          0.0451089183454136,
          0.0451089183454136,
          0.3877034539956239,
          0.1122965460043761,
          0.1122965460043761,
          0.3877034539956239,
          0.1122965460043761,
          0.3877034539956239,
          0.4588714487524593,
          0.4588714487524593,
          0.4588714487524593,
          0.0025545792330413,
          0.0797025232620401,
          0.0025545792330413,
          0.4588714487524593,
          0.0797025232620401,
          0.4588714487524593,
          0.0025545792330413,
          0.4588714487524593,
          0.0797025232620401,
          0.0337758706853386,
          0.0337758706853386,
          0.0337758706853386,
          0.7183503264420745,
          0.2140979321872483,
          0.7183503264420745,
          0.0337758706853386,
          0.2140979321872483,
          0.0337758706853386,
          0.7183503264420745,
          0.0337758706853386,
          0.2140979321872483,
          0.1836413698099279,
          0.1836413698099279,
          0.1836413698099279,
          0.0344159105781753,
          0.5983013498019689,
          0.0344159105781753,
          0.1836413698099279,
          0.5983013498019689,
          0.1836413698099279,
          0.0344159105781753,
          0.1836413698099279,
          0.5983013498019689 };

  static double c_save[] = {
          0.2500000000000000,
          0.9999999981405490,
          0.0000000006198170,
          0.0000000006198170,
          0.0000000006198170,
          0.5176763938142153,
          0.1607745353952616,
          0.1607745353952616,
          0.1607745353952616,
          0.0331704345357371,
          0.3222765218214210,
          0.3222765218214210,
          0.3222765218214210,
          0.8646732449637593,
          0.0451089183454136,
          0.0451089183454136,
          0.0451089183454136,
          0.3877034539956239,
          0.3877034539956239,
          0.3877034539956239,
          0.1122965460043761,
          0.1122965460043761,
          0.1122965460043761,
          0.0797025232620401,
          0.4588714487524593,
          0.0025545792330413,
          0.4588714487524593,
          0.0025545792330413,
          0.4588714487524593,
          0.0025545792330413,
          0.4588714487524593,
          0.0797025232620401,
          0.0797025232620401,
          0.4588714487524593,
          0.4588714487524593,
          0.2140979321872483,
          0.0337758706853386,
          0.7183503264420745,
          0.0337758706853386,
          0.7183503264420745,
          0.0337758706853386,
          0.7183503264420745,
          0.0337758706853386,
          0.2140979321872483,
          0.2140979321872483,
          0.0337758706853386,
          0.0337758706853386,
          0.5983013498019689,
          0.1836413698099279,
          0.0344159105781753,
          0.1836413698099279,
          0.0344159105781753,
          0.1836413698099279,
          0.0344159105781753,
          0.1836413698099279,
          0.5983013498019689,
          0.5983013498019689,
          0.1836413698099279,
          0.1836413698099279 };

  static double d_save[] = {
          0.2500000000000000,
          0.0000000006198171,
          0.0000000006198171,
          0.0000000006198171,
          0.9999999981405491,
          0.1607745353952614,
          0.1607745353952615,
          0.1607745353952615,
          0.5176763938142150,
          0.3222765218214210,
          0.3222765218214210,
          0.3222765218214210,
          0.0331704345357371,
          0.0451089183454134,
          0.0451089183454135,
          0.0451089183454135,
          0.8646732449637591,
          0.1122965460043760,
          0.1122965460043761,
          0.3877034539956238,
          0.3877034539956239,
          0.3877034539956239,
          0.1122965460043761,
          0.4588714487524594,
          0.0797025232620402,
          0.0797025232620402,
          0.4588714487524594,
          0.4588714487524593,
          0.0797025232620402,
          0.4588714487524593,
          0.0025545792330414,
          0.0025545792330414,
          0.4588714487524594,
          0.0025545792330414,
          0.4588714487524594,
          0.0337758706853385,
          0.2140979321872483,
          0.2140979321872483,
          0.0337758706853386,
          0.0337758706853386,
          0.2140979321872483,
          0.0337758706853386,
          0.7183503264420745,
          0.7183503264420745,
          0.0337758706853385,
          0.7183503264420745,
          0.0337758706853385,
          0.1836413698099280,
          0.5983013498019689,
          0.5983013498019688,
          0.1836413698099279,
          0.1836413698099278,
          0.5983013498019689,
          0.1836413698099279,
          0.0344159105781752,
          0.0344159105781752,
          0.1836413698099278,
          0.0344159105781753,
          0.1836413698099279 };

  double w_save[] = {
          0.0580105489124803,
          0.0000643192817593,
          0.0000643192817593,
          0.0000643192817593,
          0.0000643192817593,
          0.0231733384624255,
          0.0231733384624255,
          0.0231733384624255,
          0.0231733384624255,
          0.0295629123354293,
          0.0295629123354293,
          0.0295629123354293,
          0.0295629123354293,
          0.0080639799796162,
          0.0080639799796162,
          0.0080639799796162,
          0.0080639799796162,
          0.0381340801037025,
          0.0381340801037025,
          0.0381340801037025,
          0.0381340801037025,
          0.0381340801037025,
          0.0381340801037025,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0083844221982986,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0102345593527453,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881,
          0.0205249159679881 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule10 ( int n, double a[], double b[], double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//   rule10() returns the rule of precision 10.
//
//  Discussion:
//
//    We suppose we are given a tetrahedron T with vertices A, B, C, D.
//    We call a rule with n points, returning barycentric coordinates
//    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is 
//    approximated by Q as follows:
//
//    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D
//    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 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 a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.2500000000000000,
          0.3122500686951887,
          0.3122500686951887,
          0.0632497939144341,
          0.3122500686951887,
          0.1143096538573461,
          0.1143096538573461,
          0.6570710384279616,
          0.1143096538573461,
          0.1654860256196111,
          0.1654860256196111,
          0.4104307392189655,
          0.0136524959424580,
          0.4104307392189655,
          0.4104307392189655,
          0.0136524959424580,
          0.4104307392189655,
          0.4104307392189655,
          0.4104307392189655,
          0.0136524959424580,
          0.1654860256196111,
          0.9429887673452049,
          0.9429887673452049,
          0.0061380088247908,
          0.0447352150052137,
          0.0061380088247908,
          0.0061380088247908,
          0.0447352150052137,
          0.0061380088247908,
          0.0061380088247908,
          0.0061380088247908,
          0.0447352150052137,
          0.9429887673452049,
          0.4771903799042804,
          0.4771903799042804,
          0.1210501811455894,
          0.2807092578045408,
          0.1210501811455894,
          0.1210501811455894,
          0.2807092578045408,
          0.1210501811455894,
          0.1210501811455894,
          0.1210501811455894,
          0.2807092578045408,
          0.4771903799042804,
          0.5942562694800070,
          0.5942562694800070,
          0.0327794682164427,
          0.3401847940871077,
          0.0327794682164427,
          0.0327794682164427,
          0.3401847940871077,
          0.0327794682164427,
          0.0327794682164427,
          0.0327794682164427,
          0.3401847940871077,
          0.5942562694800070,
          0.8011772846583444,
          0.8011772846583444,
          0.0324852815648231,
          0.1338521522120095,
          0.0324852815648231,
          0.0324852815648231,
          0.1338521522120095,
          0.0324852815648231,
          0.0324852815648231,
          0.0324852815648231,
          0.1338521522120095,
          0.8011772846583444,
          0.6280718454753660,
          0.6280718454753660,
          0.1749793421839390,
          0.0219694701567559,
          0.1749793421839390,
          0.1749793421839390,
          0.0219694701567559,
          0.1749793421839390,
          0.1749793421839390,
          0.1749793421839390,
          0.0219694701567559,
          0.6280718454753660 };

  static double b_save[] = {
          0.2500000000000000,
          0.3122500686951887,
          0.0632497939144341,
          0.3122500686951887,
          0.3122500686951887,
          0.1143096538573461,
          0.6570710384279616,
          0.1143096538573461,
          0.1143096538573461,
          0.4104307392189655,
          0.4104307392189655,
          0.4104307392189655,
          0.1654860256196111,
          0.0136524959424580,
          0.1654860256196111,
          0.4104307392189655,
          0.0136524959424580,
          0.4104307392189655,
          0.1654860256196111,
          0.4104307392189655,
          0.0136524959424580,
          0.0061380088247908,
          0.0061380088247908,
          0.0061380088247908,
          0.9429887673452049,
          0.0447352150052137,
          0.9429887673452049,
          0.0061380088247908,
          0.0447352150052137,
          0.0061380088247908,
          0.9429887673452049,
          0.0061380088247908,
          0.0447352150052137,
          0.1210501811455894,
          0.1210501811455894,
          0.1210501811455894,
          0.4771903799042804,
          0.2807092578045408,
          0.4771903799042804,
          0.1210501811455894,
          0.2807092578045408,
          0.1210501811455894,
          0.4771903799042804,
          0.1210501811455894,
          0.2807092578045408,
          0.0327794682164427,
          0.0327794682164427,
          0.0327794682164427,
          0.5942562694800070,
          0.3401847940871077,
          0.5942562694800070,
          0.0327794682164427,
          0.3401847940871077,
          0.0327794682164427,
          0.5942562694800070,
          0.0327794682164427,
          0.3401847940871077,
          0.0324852815648231,
          0.0324852815648231,
          0.0324852815648231,
          0.8011772846583444,
          0.1338521522120095,
          0.8011772846583444,
          0.0324852815648231,
          0.1338521522120095,
          0.0324852815648231,
          0.8011772846583444,
          0.0324852815648231,
          0.1338521522120095,
          0.1749793421839390,
          0.1749793421839390,
          0.1749793421839390,
          0.6280718454753660,
          0.0219694701567559,
          0.6280718454753660,
          0.1749793421839390,
          0.0219694701567559,
          0.1749793421839390,
          0.6280718454753660,
          0.1749793421839390,
          0.0219694701567559 };

  static double c_save[] = {
          0.2500000000000000,
          0.0632497939144341,
          0.3122500686951887,
          0.3122500686951887,
          0.3122500686951887,
          0.6570710384279616,
          0.1143096538573461,
          0.1143096538573461,
          0.1143096538573461,
          0.0136524959424580,
          0.4104307392189655,
          0.1654860256196111,
          0.4104307392189655,
          0.1654860256196111,
          0.4104307392189655,
          0.1654860256196111,
          0.4104307392189655,
          0.0136524959424580,
          0.0136524959424580,
          0.4104307392189655,
          0.4104307392189655,
          0.0447352150052137,
          0.0061380088247908,
          0.9429887673452049,
          0.0061380088247908,
          0.9429887673452049,
          0.0061380088247908,
          0.9429887673452049,
          0.0061380088247908,
          0.0447352150052137,
          0.0447352150052137,
          0.0061380088247908,
          0.0061380088247908,
          0.2807092578045408,
          0.1210501811455894,
          0.4771903799042804,
          0.1210501811455894,
          0.4771903799042804,
          0.1210501811455894,
          0.4771903799042804,
          0.1210501811455894,
          0.2807092578045408,
          0.2807092578045408,
          0.1210501811455894,
          0.1210501811455894,
          0.3401847940871077,
          0.0327794682164427,
          0.5942562694800070,
          0.0327794682164427,
          0.5942562694800070,
          0.0327794682164427,
          0.5942562694800070,
          0.0327794682164427,
          0.3401847940871077,
          0.3401847940871077,
          0.0327794682164427,
          0.0327794682164427,
          0.1338521522120095,
          0.0324852815648231,
          0.8011772846583444,
          0.0324852815648231,
          0.8011772846583444,
          0.0324852815648231,
          0.8011772846583444,
          0.0324852815648231,
          0.1338521522120095,
          0.1338521522120095,
          0.0324852815648231,
          0.0324852815648231,
          0.0219694701567559,
          0.1749793421839390,
          0.6280718454753660,
          0.1749793421839390,
          0.6280718454753660,
          0.1749793421839390,
          0.6280718454753660,
          0.1749793421839390,
          0.0219694701567559,
          0.0219694701567559,
          0.1749793421839390,
          0.1749793421839390 };

  static double d_save[] = {
          0.2500000000000000,
          0.3122500686951886,
          0.3122500686951886,
          0.3122500686951886,
          0.0632497939144340,
          0.1143096538573461,
          0.1143096538573461,
          0.1143096538573461,
          0.6570710384279616,
          0.4104307392189654,
          0.0136524959424579,
          0.0136524959424580,
          0.4104307392189654,
          0.4104307392189655,
          0.0136524959424580,
          0.4104307392189655,
          0.1654860256196111,
          0.1654860256196111,
          0.4104307392189655,
          0.1654860256196111,
          0.4104307392189654,
          0.0061380088247907,
          0.0447352150052136,
          0.0447352150052136,
          0.0061380088247908,
          0.0061380088247907,
          0.0447352150052136,
          0.0061380088247908,
          0.9429887673452048,
          0.9429887673452049,
          0.0061380088247907,
          0.9429887673452049,
          0.0061380088247907,
          0.1210501811455894,
          0.2807092578045408,
          0.2807092578045408,
          0.1210501811455894,
          0.1210501811455894,
          0.2807092578045408,
          0.1210501811455894,
          0.4771903799042804,
          0.4771903799042804,
          0.1210501811455894,
          0.4771903799042804,
          0.1210501811455894,
          0.0327794682164427,
          0.3401847940871077,
          0.3401847940871078,
          0.0327794682164427,
          0.0327794682164427,
          0.3401847940871077,
          0.0327794682164427,
          0.5942562694800071,
          0.5942562694800071,
          0.0327794682164427,
          0.5942562694800071,
          0.0327794682164427,
          0.0324852815648231,
          0.1338521522120095,
          0.1338521522120094,
          0.0324852815648231,
          0.0324852815648230,
          0.1338521522120095,
          0.0324852815648230,
          0.8011772846583443,
          0.8011772846583443,
          0.0324852815648230,
          0.8011772846583443,
          0.0324852815648231,
          0.1749793421839391,
          0.0219694701567560,
          0.0219694701567561,
          0.1749793421839391,
          0.1749793421839392,
          0.0219694701567560,
          0.1749793421839392,
          0.6280718454753662,
          0.6280718454753662,
          0.1749793421839391,
          0.6280718454753662,
          0.1749793421839391 };

  double w_save[] = {
          0.0473997735560207,
          0.0269370599922687,
          0.0269370599922687,
          0.0269370599922687,
          0.0269370599922687,
          0.0098691597167934,
          0.0098691597167934,
          0.0098691597167934,
          0.0098691597167934,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0113938812201952,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0003619443443393,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0257397319804561,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0101358716797558,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0065761472770359,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620,
          0.0129070357988620 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, d_save, d );
  r8vec_copy ( n, w_save, w );

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

double tetrahedron_unit_monomial_integral ( int expon[] )

//****************************************************************************80
//
//  Purpose:
//
//    tetrahedron_unit_monomial_integral() integrates a monomial over the unit tetrahedron.
//
//  Discussion:
//
//    This routine evaluates a monomial of the form
//
//      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
//
//    where the exponents are nonnegative integers.  Note that
//    if the combination 0^0 is encountered, it should be treated
//    as 1.
//
//    Integral ( over unit tetrahedron ) x^l y^m z^n dx dy 
//    = l! * m! * n! / ( l + m + n + 3 )!
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    04 July 2007
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int EXPON[3], the exponents.
//
//  Output:
//
//    double tetrahedron_unit_monomial_integral, the value of the integral.
//
{
  int i;
  int k;
  double value;
//
//  The first computation ends with VALUE = 1.0;
//
  value = 1.0;

  k = 0;

  for ( i = 1; i <= expon[0]; i++ )
  {
    k = k + 1;
//  value = value * ( double ) ( i ) / ( double ) ( k );
  }

  for ( i = 1; i <= expon[1]; i++ )
  {
    k = k + 1;
    value = value * ( double ) ( i ) / ( double ) ( k );
  }

  for ( i = 1; i <= expon[2]; i++ )
  {
    k = k + 1;
    value = value * ( double ) ( i ) / ( double ) ( k );
  }

  k = k + 1;
  value = value / ( double ) ( k );

  k = k + 1;
  value = value / ( double ) ( k );

  k = k + 1;
  value = value / ( double ) ( k );

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

double tetrahedron_unit_volume ( )

//****************************************************************************80
//
//  Purpose:
//
//    tetrahedron_unit_volume() computes the volume of a unit tetrahedron.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    23 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    double tetrahedron_unit_volume: the volume.
//
{
  double volume;

  volume = 1.0 / 6.0;

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

double tetrahedron_volume ( double tetra[4*3] )

//****************************************************************************80
//
//  Purpose:
//
//    tetrahedron_volume() computes the volume of a tetrahedron.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    13 April 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double tetra[4*3]: the vertex coordinates.
//
//  Output:
//
//    double tetrahedron_volume: the volume.
//
{
  double a[4*4];
  int i;
  int j;
  double volume;
//
//  Stored by rows.
//  (Warning: I used to store 2D arrays by columns instead.)
//
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 3; j++ )
    { 
      a[i*4+j] = tetra[i*3+j];
    }
    a[i*4+3] = 1.0;
  }

  volume = fabs ( r8mat_det_4d ( a ) ) / 6.0;

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

void tetrahedron_witherden_rule ( int p, int n, double a[], double b[], 
  double c[], double d[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    tetrahedron_witherden_rule() returns a tetrahedron quadrature rule of given precision.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 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:
//
//    integer p: the precision, 0 <= p <= 10.
//
//    integer n: the order of the rule.
//
//  Output:
//
//    real a(n), b(n), c(n), d(n): the barycentric coordinates of quadrature points.
//
//    real w(n): the weights of quadrature points.
//
{
  if ( p < 0 )
  {
    cout << "\n";
    cout << "tetrahedron_witherden_rule(): Fatal error!\n";
    cout << "  Input p < 0.\n";
    exit ( 1 );
  }

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

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

  return;
}