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

using namespace std;

# include "triangle_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

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:
//
//    23 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 <= 20.
//
//  Output:
//
//    int rule_order: the order of the rule.
//
{
  int order;
  static int order_save[] = {
      1,
      1,   3,   6,   6,   7, 
     12,  15,  16,  19,  25, 
     28,  33,  37,  42,  49, 
     55,  60,  67,  73,  79 };

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

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

  order = order_save[p];

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

void rule00 ( int n, double a[], double b[], double c[], double w[] )

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

  static double b_save[] = {
          0.3333333333333334 };

  static double c_save[] = {
          0.3333333333333333 };

  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, w_save, w );

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

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

  static double b_save[] = {
          0.3333333333333334 };

  static double c_save[] = {
          0.3333333333333333 };

  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, w_save, w );

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

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

  static double b_save[] = {
          0.6666666666666666,
          0.1666666666666667,
          0.1666666666666667 };

  static double c_save[] = {
          0.1666666666666666,
          0.1666666666666667,
          0.6666666666666666 };

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

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

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

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

  static double b_save[] = {
          0.1081030181680702,
          0.4459484909159649,
          0.4459484909159649,
          0.8168475729804585,
          0.0915762135097707,
          0.0915762135097707 };

  static double c_save[] = {
          0.4459484909159649,
          0.4459484909159649,
          0.1081030181680702,
          0.0915762135097707,
          0.0915762135097707,
          0.8168475729804585 };

  double w_save[] = {
          0.2233815896780115,
          0.2233815896780115,
          0.2233815896780115,
          0.1099517436553219,
          0.1099517436553219,
          0.1099517436553219 };

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

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

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

  static double b_save[] = {
          0.1081030181680702,
          0.4459484909159649,
          0.4459484909159649,
          0.8168475729804585,
          0.0915762135097707,
          0.0915762135097707 };

  static double c_save[] = {
          0.4459484909159649,
          0.4459484909159649,
          0.1081030181680702,
          0.0915762135097707,
          0.0915762135097707,
          0.8168475729804585 };

  double w_save[] = {
          0.2233815896780115,
          0.2233815896780115,
          0.2233815896780115,
          0.1099517436553219,
          0.1099517436553219,
          0.1099517436553219 };

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

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

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

  static double b_save[] = {
          0.3333333333333334,
          0.7974269853530873,
          0.1012865073234563,
          0.1012865073234563,
          0.0597158717897698,
          0.4701420641051151,
          0.4701420641051151 };

  static double c_save[] = {
          0.3333333333333333,
          0.1012865073234563,
          0.1012865073234563,
          0.7974269853530873,
          0.4701420641051151,
          0.4701420641051151,
          0.0597158717897698 };

  double w_save[] = {
          0.2250000000000000,
          0.1259391805448271,
          0.1259391805448271,
          0.1259391805448271,
          0.1323941527885062,
          0.1323941527885062,
          0.1323941527885062 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule06() returns the rule of precision 6.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.0630890144915022,
          0.8738219710169955,
          0.0630890144915022,
          0.2492867451709104,
          0.5014265096581791,
          0.2492867451709104,
          0.0531450498448169,
          0.6365024991213987,
          0.3103524510337844,
          0.6365024991213987,
          0.3103524510337844,
          0.0531450498448169 };

  static double b_save[] = {
          0.8738219710169955,
          0.0630890144915022,
          0.0630890144915022,
          0.5014265096581791,
          0.2492867451709104,
          0.2492867451709104,
          0.6365024991213987,
          0.0531450498448169,
          0.6365024991213987,
          0.3103524510337844,
          0.0531450498448169,
          0.3103524510337844 };

  static double c_save[] = {
          0.0630890144915022,
          0.0630890144915022,
          0.8738219710169955,
          0.2492867451709104,
          0.2492867451709104,
          0.5014265096581791,
          0.3103524510337844,
          0.3103524510337844,
          0.0531450498448169,
          0.0531450498448169,
          0.6365024991213987,
          0.6365024991213987 };

  double w_save[] = {
          0.0508449063702068,
          0.0508449063702068,
          0.0508449063702068,
          0.1167862757263794,
          0.1167862757263794,
          0.1167862757263794,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule07() returns the rule of precision 7.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.0337306485545878,
          0.9325387028908243,
          0.0337306485545878,
          0.2415773825954036,
          0.5168452348091929,
          0.2415773825954036,
          0.4743096925047182,
          0.0513806149905635,
          0.4743096925047182,
          0.0470366446525952,
          0.7542800405500532,
          0.1986833147973516,
          0.7542800405500532,
          0.1986833147973516,
          0.0470366446525952 };

  static double b_save[] = {
          0.9325387028908243,
          0.0337306485545878,
          0.0337306485545878,
          0.5168452348091929,
          0.2415773825954036,
          0.2415773825954036,
          0.0513806149905635,
          0.4743096925047182,
          0.4743096925047182,
          0.7542800405500532,
          0.0470366446525952,
          0.7542800405500532,
          0.1986833147973516,
          0.0470366446525952,
          0.1986833147973516 };

  static double c_save[] = {
          0.0337306485545878,
          0.0337306485545878,
          0.9325387028908243,
          0.2415773825954035,
          0.2415773825954035,
          0.5168452348091929,
          0.4743096925047183,
          0.4743096925047182,
          0.0513806149905636,
          0.1986833147973517,
          0.1986833147973517,
          0.0470366446525953,
          0.0470366446525953,
          0.7542800405500532,
          0.7542800405500532 };

  double w_save[] = {
          0.0165450501107921,
          0.0165450501107921,
          0.0165450501107921,
          0.1279441712301556,
          0.1279441712301556,
          0.1279441712301556,
          0.0770866461859861,
          0.0770866461859861,
          0.0770866461859861,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule08() returns the rule of precision 8.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.4592925882927231,
          0.0814148234145537,
          0.4592925882927231,
          0.1705693077517602,
          0.6588613844964796,
          0.1705693077517602,
          0.0505472283170310,
          0.8989055433659381,
          0.0505472283170310,
          0.0083947774099576,
          0.7284923929554042,
          0.2631128296346381,
          0.7284923929554042,
          0.2631128296346381,
          0.0083947774099576 };

  static double b_save[] = {
          0.3333333333333334,
          0.0814148234145537,
          0.4592925882927231,
          0.4592925882927231,
          0.6588613844964796,
          0.1705693077517602,
          0.1705693077517602,
          0.8989055433659381,
          0.0505472283170310,
          0.0505472283170310,
          0.7284923929554042,
          0.0083947774099576,
          0.7284923929554042,
          0.2631128296346381,
          0.0083947774099576,
          0.2631128296346381 };

  static double c_save[] = {
          0.3333333333333333,
          0.4592925882927232,
          0.4592925882927231,
          0.0814148234145537,
          0.1705693077517603,
          0.1705693077517603,
          0.6588613844964796,
          0.0505472283170310,
          0.0505472283170310,
          0.8989055433659381,
          0.2631128296346381,
          0.2631128296346381,
          0.0083947774099576,
          0.0083947774099576,
          0.7284923929554044,
          0.7284923929554044 };

  double w_save[] = {
          0.1443156076777872,
          0.0950916342672846,
          0.0950916342672846,
          0.0950916342672846,
          0.1032173705347182,
          0.1032173705347182,
          0.1032173705347182,
          0.0324584976231981,
          0.0324584976231981,
          0.0324584976231981,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule09() returns the rule of precision 9.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.4370895914929366,
          0.1258208170141267,
          0.4370895914929366,
          0.1882035356190327,
          0.6235929287619345,
          0.1882035356190327,
          0.4896825191987376,
          0.0206349616025248,
          0.4896825191987376,
          0.0447295133944527,
          0.9105409732110945,
          0.0447295133944527,
          0.0368384120547363,
          0.7411985987844980,
          0.2219629891607657,
          0.7411985987844980,
          0.2219629891607657,
          0.0368384120547363 };

  static double b_save[] = {
          0.3333333333333334,
          0.1258208170141267,
          0.4370895914929366,
          0.4370895914929366,
          0.6235929287619345,
          0.1882035356190327,
          0.1882035356190327,
          0.0206349616025248,
          0.4896825191987376,
          0.4896825191987376,
          0.9105409732110945,
          0.0447295133944527,
          0.0447295133944527,
          0.7411985987844980,
          0.0368384120547363,
          0.7411985987844980,
          0.2219629891607657,
          0.0368384120547363,
          0.2219629891607657 };

  static double c_save[] = {
          0.3333333333333333,
          0.4370895914929367,
          0.4370895914929367,
          0.1258208170141267,
          0.1882035356190327,
          0.1882035356190328,
          0.6235929287619345,
          0.4896825191987376,
          0.4896825191987376,
          0.0206349616025248,
          0.0447295133944527,
          0.0447295133944527,
          0.9105409732110946,
          0.2219629891607657,
          0.2219629891607657,
          0.0368384120547363,
          0.0368384120547363,
          0.7411985987844980,
          0.7411985987844980 };

  double w_save[] = {
          0.0971357962827988,
          0.0778275410047743,
          0.0778275410047743,
          0.0778275410047743,
          0.0796477389272102,
          0.0796477389272102,
          0.0796477389272102,
          0.0313347002271391,
          0.0313347002271391,
          0.0313347002271391,
          0.0255776756586980,
          0.0255776756586980,
          0.0255776756586980,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule10() returns the rule of precision 10.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.0320553732169435,
          0.9358892535661130,
          0.0320553732169435,
          0.1421611010565644,
          0.7156777978868712,
          0.1421611010565644,
          0.3218129952888354,
          0.5300541189273440,
          0.1481328857838206,
          0.5300541189273440,
          0.1481328857838206,
          0.3218129952888354,
          0.0296198894887298,
          0.6012333286834592,
          0.3691467818278110,
          0.6012333286834592,
          0.3691467818278110,
          0.0296198894887298,
          0.0283676653399385,
          0.8079306009228790,
          0.1637017337371825,
          0.8079306009228790,
          0.1637017337371825,
          0.0283676653399385 };

  static double b_save[] = {
          0.3333333333333334,
          0.9358892535661130,
          0.0320553732169435,
          0.0320553732169435,
          0.7156777978868712,
          0.1421611010565644,
          0.1421611010565644,
          0.5300541189273440,
          0.3218129952888354,
          0.5300541189273440,
          0.1481328857838206,
          0.3218129952888354,
          0.1481328857838206,
          0.6012333286834592,
          0.0296198894887298,
          0.6012333286834592,
          0.3691467818278110,
          0.0296198894887298,
          0.3691467818278110,
          0.8079306009228790,
          0.0283676653399385,
          0.8079306009228790,
          0.1637017337371825,
          0.0283676653399385,
          0.1637017337371825 };

  static double c_save[] = {
          0.3333333333333333,
          0.0320553732169435,
          0.0320553732169435,
          0.9358892535661130,
          0.1421611010565644,
          0.1421611010565644,
          0.7156777978868712,
          0.1481328857838204,
          0.1481328857838206,
          0.3218129952888354,
          0.3218129952888354,
          0.5300541189273440,
          0.5300541189273440,
          0.3691467818278109,
          0.3691467818278109,
          0.0296198894887297,
          0.0296198894887297,
          0.6012333286834592,
          0.6012333286834592,
          0.1637017337371824,
          0.1637017337371824,
          0.0283676653399385,
          0.0283676653399385,
          0.8079306009228790,
          0.8079306009228790 };

  double w_save[] = {
          0.0817433291462860,
          0.0133529688131496,
          0.0133529688131496,
          0.0133529688131496,
          0.0459579636047447,
          0.0459579636047447,
          0.0459579636047447,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule11() returns the rule of precision 11.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.0284854176143719,
          0.9430291647712562,
          0.0284854176143719,
          0.2102199567031783,
          0.5795600865936434,
          0.2102199567031783,
          0.1026354827122464,
          0.7947290345755071,
          0.1026354827122464,
          0.4958919009658909,
          0.0082161980682182,
          0.4958919009658909,
          0.4384659267643522,
          0.1230681464712956,
          0.4384659267643522,
          0.1493247886520824,
          0.8433497836618531,
          0.0073254276860645,
          0.8433497836618531,
          0.0073254276860645,
          0.1493247886520824,
          0.0460105001654300,
          0.6644083741968642,
          0.2895811256377059,
          0.6644083741968642,
          0.2895811256377059,
          0.0460105001654300 };

  static double b_save[] = {
          0.3333333333333334,
          0.9430291647712562,
          0.0284854176143719,
          0.0284854176143719,
          0.5795600865936434,
          0.2102199567031783,
          0.2102199567031783,
          0.7947290345755071,
          0.1026354827122464,
          0.1026354827122464,
          0.0082161980682182,
          0.4958919009658909,
          0.4958919009658909,
          0.1230681464712956,
          0.4384659267643522,
          0.4384659267643522,
          0.8433497836618531,
          0.1493247886520824,
          0.8433497836618531,
          0.0073254276860645,
          0.1493247886520824,
          0.0073254276860645,
          0.6644083741968642,
          0.0460105001654300,
          0.6644083741968642,
          0.2895811256377059,
          0.0460105001654300,
          0.2895811256377059 };

  static double c_save[] = {
          0.3333333333333333,
          0.0284854176143718,
          0.0284854176143718,
          0.9430291647712562,
          0.2102199567031782,
          0.2102199567031783,
          0.5795600865936434,
          0.1026354827122464,
          0.1026354827122464,
          0.7947290345755071,
          0.4958919009658909,
          0.4958919009658910,
          0.0082161980682182,
          0.4384659267643523,
          0.4384659267643523,
          0.1230681464712956,
          0.0073254276860645,
          0.0073254276860646,
          0.1493247886520823,
          0.1493247886520824,
          0.8433497836618531,
          0.8433497836618531,
          0.2895811256377059,
          0.2895811256377059,
          0.0460105001654300,
          0.0460105001654300,
          0.6644083741968642,
          0.6644083741968642 };

  double w_save[] = {
          0.0857611797322242,
          0.0104318705128947,
          0.0104318705128947,
          0.0104318705128947,
          0.0705156841117166,
          0.0705156841117166,
          0.0705156841117166,
          0.0386307592370193,
          0.0386307592370193,
          0.0386307592370193,
          0.0166062730545854,
          0.0166062730545854,
          0.0166062730545854,
          0.0673161540794683,
          0.0673161540794683,
          0.0673161540794683,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule12() returns the rule of precision 12.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.4882037509455415,
          0.0235924981089169,
          0.4882037509455415,
          0.1092578276593543,
          0.7814843446812914,
          0.1092578276593543,
          0.2714625070149261,
          0.4570749859701478,
          0.2714625070149261,
          0.0246463634363356,
          0.9507072731273288,
          0.0246463634363356,
          0.4401116486585931,
          0.1197767026828138,
          0.4401116486585931,
          0.2916556797383409,
          0.6853101639063919,
          0.0230341563552671,
          0.6853101639063919,
          0.0230341563552671,
          0.2916556797383409,
          0.1162960196779266,
          0.6282497516835561,
          0.2554542286385174,
          0.6282497516835561,
          0.2554542286385174,
          0.1162960196779266,
          0.8513377925102401,
          0.1272797172335894,
          0.0213824902561706,
          0.1272797172335894,
          0.0213824902561706,
          0.8513377925102401 };

  static double b_save[] = {
          0.0235924981089169,
          0.4882037509455415,
          0.4882037509455415,
          0.7814843446812914,
          0.1092578276593543,
          0.1092578276593543,
          0.4570749859701478,
          0.2714625070149261,
          0.2714625070149261,
          0.9507072731273288,
          0.0246463634363356,
          0.0246463634363356,
          0.1197767026828138,
          0.4401116486585931,
          0.4401116486585931,
          0.6853101639063919,
          0.2916556797383409,
          0.6853101639063919,
          0.0230341563552671,
          0.2916556797383409,
          0.0230341563552671,
          0.6282497516835561,
          0.1162960196779266,
          0.6282497516835561,
          0.2554542286385174,
          0.1162960196779266,
          0.2554542286385174,
          0.1272797172335894,
          0.8513377925102401,
          0.1272797172335894,
          0.0213824902561706,
          0.8513377925102401,
          0.0213824902561706 };

  static double c_save[] = {
          0.4882037509455416,
          0.4882037509455416,
          0.0235924981089169,
          0.1092578276593543,
          0.1092578276593543,
          0.7814843446812915,
          0.2714625070149261,
          0.2714625070149261,
          0.4570749859701478,
          0.0246463634363356,
          0.0246463634363356,
          0.9507072731273288,
          0.4401116486585932,
          0.4401116486585931,
          0.1197767026828138,
          0.0230341563552672,
          0.0230341563552672,
          0.2916556797383409,
          0.2916556797383409,
          0.6853101639063919,
          0.6853101639063919,
          0.2554542286385173,
          0.2554542286385173,
          0.1162960196779266,
          0.1162960196779266,
          0.6282497516835561,
          0.6282497516835561,
          0.0213824902561706,
          0.0213824902561706,
          0.8513377925102400,
          0.8513377925102400,
          0.1272797172335893,
          0.1272797172335893 };

  double w_save[] = {
          0.0242668380814520,
          0.0242668380814520,
          0.0242668380814520,
          0.0284860520688775,
          0.0284860520688775,
          0.0284860520688775,
          0.0625412131959028,
          0.0625412131959028,
          0.0625412131959028,
          0.0079316425099736,
          0.0079316425099736,
          0.0079316425099736,
          0.0499183349280609,
          0.0499183349280609,
          0.0499183349280609,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule13() returns the rule of precision 13.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.4890769464525394,
          0.0218461070949213,
          0.4890769464525394,
          0.2213722862918329,
          0.5572554274163342,
          0.2213722862918329,
          0.4269414142598004,
          0.1461171714803992,
          0.4269414142598004,
          0.0215096811088432,
          0.9569806377823136,
          0.0215096811088432,
          0.0878954830321973,
          0.7485071158999522,
          0.1635974010678505,
          0.7485071158999522,
          0.1635974010678505,
          0.0878954830321973,
          0.1109220428034634,
          0.8647077702954428,
          0.0243701869010938,
          0.8647077702954428,
          0.0243701869010938,
          0.1109220428034634,
          0.3084417608921178,
          0.6235459955536755,
          0.0680122435542067,
          0.6235459955536755,
          0.0680122435542067,
          0.3084417608921178,
          0.0051263891023824,
          0.7223577931241880,
          0.2725158177734297,
          0.7223577931241880,
          0.2725158177734297,
          0.0051263891023824 };

  static double b_save[] = {
          0.3333333333333334,
          0.0218461070949213,
          0.4890769464525394,
          0.4890769464525394,
          0.5572554274163342,
          0.2213722862918329,
          0.2213722862918329,
          0.1461171714803992,
          0.4269414142598004,
          0.4269414142598004,
          0.9569806377823136,
          0.0215096811088432,
          0.0215096811088432,
          0.7485071158999522,
          0.0878954830321973,
          0.7485071158999522,
          0.1635974010678505,
          0.0878954830321973,
          0.1635974010678505,
          0.8647077702954428,
          0.1109220428034634,
          0.8647077702954428,
          0.0243701869010938,
          0.1109220428034634,
          0.0243701869010938,
          0.6235459955536755,
          0.3084417608921178,
          0.6235459955536755,
          0.0680122435542067,
          0.3084417608921178,
          0.0680122435542067,
          0.7223577931241880,
          0.0051263891023824,
          0.7223577931241880,
          0.2725158177734297,
          0.0051263891023824,
          0.2725158177734297 };

  static double c_save[] = {
          0.3333333333333333,
          0.4890769464525393,
          0.4890769464525394,
          0.0218461070949213,
          0.2213722862918328,
          0.2213722862918329,
          0.5572554274163342,
          0.4269414142598004,
          0.4269414142598004,
          0.1461171714803992,
          0.0215096811088432,
          0.0215096811088433,
          0.9569806377823137,
          0.1635974010678505,
          0.1635974010678505,
          0.0878954830321973,
          0.0878954830321973,
          0.7485071158999522,
          0.7485071158999522,
          0.0243701869010938,
          0.0243701869010938,
          0.1109220428034634,
          0.1109220428034634,
          0.8647077702954428,
          0.8647077702954428,
          0.0680122435542067,
          0.0680122435542068,
          0.3084417608921177,
          0.3084417608921178,
          0.6235459955536755,
          0.6235459955536755,
          0.2725158177734297,
          0.2725158177734296,
          0.0051263891023823,
          0.0051263891023823,
          0.7223577931241879,
          0.7223577931241879 };

  double w_save[] = {
          0.0679600365868316,
          0.0239944019288947,
          0.0239944019288947,
          0.0239944019288947,
          0.0582784851192000,
          0.0582784851192000,
          0.0582784851192000,
          0.0556019675304533,
          0.0556019675304533,
          0.0556019675304533,
          0.0060523371035392,
          0.0060523371035392,
          0.0060523371035392,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule14() returns the rule of precision 14.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.1772055324125434,
          0.6455889351749131,
          0.1772055324125434,
          0.4176447193404539,
          0.1647105613190922,
          0.4176447193404539,
          0.0617998830908726,
          0.8764002338182548,
          0.0617998830908726,
          0.4889639103621786,
          0.0220721792756427,
          0.4889639103621786,
          0.2734775283088386,
          0.4530449433823227,
          0.2734775283088386,
          0.0193909612487010,
          0.9612180775025979,
          0.0193909612487010,
          0.2983728821362578,
          0.6869801678080878,
          0.0146469500556544,
          0.6869801678080878,
          0.0146469500556544,
          0.2983728821362578,
          0.0571247574036479,
          0.7706085547749965,
          0.1722666878213556,
          0.7706085547749965,
          0.1722666878213556,
          0.0571247574036479,
          0.3368614597963450,
          0.5702222908466832,
          0.0929162493569718,
          0.5702222908466832,
          0.0929162493569718,
          0.3368614597963450,
          0.0012683309328720,
          0.8797571713701711,
          0.1189744976969568,
          0.8797571713701711,
          0.1189744976969568,
          0.0012683309328720 };

  static double b_save[] = {
          0.6455889351749131,
          0.1772055324125434,
          0.1772055324125434,
          0.1647105613190922,
          0.4176447193404539,
          0.4176447193404539,
          0.8764002338182548,
          0.0617998830908726,
          0.0617998830908726,
          0.0220721792756427,
          0.4889639103621786,
          0.4889639103621786,
          0.4530449433823227,
          0.2734775283088386,
          0.2734775283088386,
          0.9612180775025979,
          0.0193909612487010,
          0.0193909612487010,
          0.6869801678080878,
          0.2983728821362578,
          0.6869801678080878,
          0.0146469500556544,
          0.2983728821362578,
          0.0146469500556544,
          0.7706085547749965,
          0.0571247574036479,
          0.7706085547749965,
          0.1722666878213556,
          0.0571247574036479,
          0.1722666878213556,
          0.5702222908466832,
          0.3368614597963450,
          0.5702222908466832,
          0.0929162493569718,
          0.3368614597963450,
          0.0929162493569718,
          0.8797571713701711,
          0.0012683309328720,
          0.8797571713701711,
          0.1189744976969568,
          0.0012683309328720,
          0.1189744976969568 };

  static double c_save[] = {
          0.1772055324125434,
          0.1772055324125434,
          0.6455889351749131,
          0.4176447193404539,
          0.4176447193404539,
          0.1647105613190921,
          0.0617998830908726,
          0.0617998830908726,
          0.8764002338182548,
          0.4889639103621787,
          0.4889639103621787,
          0.0220721792756428,
          0.2734775283088386,
          0.2734775283088386,
          0.4530449433823227,
          0.0193909612487011,
          0.0193909612487011,
          0.9612180775025979,
          0.0146469500556543,
          0.0146469500556544,
          0.2983728821362577,
          0.2983728821362578,
          0.6869801678080878,
          0.6869801678080878,
          0.1722666878213557,
          0.1722666878213557,
          0.0571247574036480,
          0.0571247574036480,
          0.7706085547749966,
          0.7706085547749966,
          0.0929162493569718,
          0.0929162493569718,
          0.3368614597963450,
          0.3368614597963450,
          0.5702222908466832,
          0.5702222908466832,
          0.1189744976969569,
          0.1189744976969569,
          0.0012683309328720,
          0.0012683309328720,
          0.8797571713701711,
          0.8797571713701711 };

  double w_save[] = {
          0.0421625887369930,
          0.0421625887369930,
          0.0421625887369930,
          0.0327883535441253,
          0.0327883535441253,
          0.0327883535441253,
          0.0144336996697767,
          0.0144336996697767,
          0.0144336996697767,
          0.0218835813694289,
          0.0218835813694289,
          0.0218835813694289,
          0.0517741045072916,
          0.0517741045072916,
          0.0517741045072916,
          0.0049234036024001,
          0.0049234036024001,
          0.0049234036024001,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule15() returns the rule of precision 15.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.4053622141339755,
          0.1892755717320491,
          0.4053622141339755,
          0.0701735528999860,
          0.8596528942000279,
          0.0701735528999860,
          0.4741706814380198,
          0.0516586371239604,
          0.4741706814380198,
          0.2263787134203496,
          0.5472425731593008,
          0.2263787134203496,
          0.4949969567691262,
          0.0100060864617476,
          0.4949969567691262,
          0.0158117262509886,
          0.9683765474980227,
          0.0158117262509886,
          0.0183761123856811,
          0.6669756448018681,
          0.3146482428124509,
          0.6669756448018681,
          0.3146482428124509,
          0.0183761123856811,
          0.0091392370373084,
          0.9199121577262361,
          0.0709486052364555,
          0.9199121577262361,
          0.0709486052364555,
          0.0091392370373084,
          0.1905355894763939,
          0.7152223569314506,
          0.0942420535921554,
          0.7152223569314506,
          0.0942420535921554,
          0.1905355894763939,
          0.1680686452224144,
          0.8132926410494192,
          0.0186387137281664,
          0.8132926410494192,
          0.0186387137281664,
          0.1680686452224144,
          0.3389506114752772,
          0.5652526648771142,
          0.0957967236476086,
          0.5652526648771142,
          0.0957967236476086,
          0.3389506114752772 };

  static double b_save[] = {
          0.3333333333333334,
          0.1892755717320491,
          0.4053622141339755,
          0.4053622141339755,
          0.8596528942000279,
          0.0701735528999860,
          0.0701735528999860,
          0.0516586371239604,
          0.4741706814380198,
          0.4741706814380198,
          0.5472425731593008,
          0.2263787134203496,
          0.2263787134203496,
          0.0100060864617476,
          0.4949969567691262,
          0.4949969567691262,
          0.9683765474980227,
          0.0158117262509886,
          0.0158117262509886,
          0.6669756448018681,
          0.0183761123856811,
          0.6669756448018681,
          0.3146482428124509,
          0.0183761123856811,
          0.3146482428124509,
          0.9199121577262361,
          0.0091392370373084,
          0.9199121577262361,
          0.0709486052364555,
          0.0091392370373084,
          0.0709486052364555,
          0.7152223569314506,
          0.1905355894763939,
          0.7152223569314506,
          0.0942420535921554,
          0.1905355894763939,
          0.0942420535921554,
          0.8132926410494192,
          0.1680686452224144,
          0.8132926410494192,
          0.0186387137281664,
          0.1680686452224144,
          0.0186387137281664,
          0.5652526648771142,
          0.3389506114752772,
          0.5652526648771142,
          0.0957967236476086,
          0.3389506114752772,
          0.0957967236476086 };

  static double c_save[] = {
          0.3333333333333333,
          0.4053622141339754,
          0.4053622141339754,
          0.1892755717320490,
          0.0701735528999861,
          0.0701735528999861,
          0.8596528942000280,
          0.4741706814380198,
          0.4741706814380198,
          0.0516586371239605,
          0.2263787134203497,
          0.2263787134203497,
          0.5472425731593008,
          0.4949969567691261,
          0.4949969567691263,
          0.0100060864617476,
          0.0158117262509887,
          0.0158117262509887,
          0.9683765474980227,
          0.3146482428124509,
          0.3146482428124509,
          0.0183761123856810,
          0.0183761123856810,
          0.6669756448018680,
          0.6669756448018681,
          0.0709486052364555,
          0.0709486052364554,
          0.0091392370373085,
          0.0091392370373083,
          0.9199121577262361,
          0.9199121577262361,
          0.0942420535921553,
          0.0942420535921554,
          0.1905355894763939,
          0.1905355894763940,
          0.7152223569314506,
          0.7152223569314506,
          0.0186387137281663,
          0.0186387137281664,
          0.1680686452224143,
          0.1680686452224144,
          0.8132926410494192,
          0.8132926410494192,
          0.0957967236476086,
          0.0957967236476086,
          0.3389506114752772,
          0.3389506114752772,
          0.5652526648771142,
          0.5652526648771142 };

  double w_save[] = {
          0.0443353873821841,
          0.0427137815714606,
          0.0427137815714606,
          0.0427137815714606,
          0.0164447375626252,
          0.0164447375626252,
          0.0164447375626252,
          0.0173961480007634,
          0.0173961480007634,
          0.0173961480007634,
          0.0467833617287096,
          0.0467833617287096,
          0.0467833617287096,
          0.0095738461824601,
          0.0095738461824601,
          0.0095738461824601,
          0.0029607746379054,
          0.0029607746379054,
          0.0029607746379054,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule16() returns the rule of precision 16.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.2459900704671417,
          0.5080198590657166,
          0.2459900704671417,
          0.4155848968854205,
          0.1688302062291590,
          0.4155848968854205,
          0.0853555665867003,
          0.8292888668265994,
          0.0853555665867003,
          0.1619186441912712,
          0.6761627116174576,
          0.1619186441912712,
          0.5000000000000000,
          0.0000000000000000,
          0.5000000000000000,
          0.4752807275459421,
          0.0494385449081158,
          0.4752807275459421,
          0.0547551749147031,
          0.7541700614447677,
          0.1910747636405291,
          0.7541700614447677,
          0.1910747636405291,
          0.0547551749147031,
          0.0232034277688137,
          0.9682443680309587,
          0.0085522042002276,
          0.9682443680309587,
          0.0085522042002276,
          0.0232034277688137,
          0.0189317782804059,
          0.6493036982454464,
          0.3317645234741476,
          0.6493036982454464,
          0.3317645234741476,
          0.0189317782804059,
          0.0190301297436974,
          0.9002737032704295,
          0.0806961669858730,
          0.9002737032704295,
          0.0806961669858730,
          0.0190301297436974,
          0.1026061902393981,
          0.5891488405642479,
          0.3082449691963540,
          0.5891488405642479,
          0.3082449691963540,
          0.1026061902393981,
          0.0059363500168222,
          0.8066218674993957,
          0.1874417824837821,
          0.8066218674993957,
          0.1874417824837821,
          0.0059363500168222 };

  static double b_save[] = {
          0.3333333333333334,
          0.5080198590657166,
          0.2459900704671417,
          0.2459900704671417,
          0.1688302062291590,
          0.4155848968854205,
          0.4155848968854205,
          0.8292888668265994,
          0.0853555665867003,
          0.0853555665867003,
          0.6761627116174576,
          0.1619186441912712,
          0.1619186441912712,
          0.0000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.0494385449081158,
          0.4752807275459421,
          0.4752807275459421,
          0.7541700614447677,
          0.0547551749147031,
          0.7541700614447677,
          0.1910747636405291,
          0.0547551749147031,
          0.1910747636405291,
          0.9682443680309587,
          0.0232034277688137,
          0.9682443680309587,
          0.0085522042002276,
          0.0232034277688137,
          0.0085522042002276,
          0.6493036982454464,
          0.0189317782804059,
          0.6493036982454464,
          0.3317645234741476,
          0.0189317782804059,
          0.3317645234741476,
          0.9002737032704295,
          0.0190301297436974,
          0.9002737032704295,
          0.0806961669858730,
          0.0190301297436974,
          0.0806961669858730,
          0.5891488405642479,
          0.1026061902393981,
          0.5891488405642479,
          0.3082449691963540,
          0.1026061902393981,
          0.3082449691963540,
          0.8066218674993957,
          0.0059363500168222,
          0.8066218674993957,
          0.1874417824837821,
          0.0059363500168222,
          0.1874417824837821 };

  static double c_save[] = {
          0.3333333333333333,
          0.2459900704671418,
          0.2459900704671417,
          0.5080198590657166,
          0.4155848968854206,
          0.4155848968854206,
          0.1688302062291590,
          0.0853555665867003,
          0.0853555665867003,
          0.8292888668265993,
          0.1619186441912712,
          0.1619186441912712,
          0.6761627116174576,
          0.5000000000000000,
          0.5000000000000000,
          0.0000000000000000,
          0.4752807275459421,
          0.4752807275459421,
          0.0494385449081158,
          0.1910747636405291,
          0.1910747636405292,
          0.0547551749147032,
          0.0547551749147033,
          0.7541700614447678,
          0.7541700614447678,
          0.0085522042002275,
          0.0085522042002275,
          0.0232034277688137,
          0.0232034277688137,
          0.9682443680309587,
          0.9682443680309587,
          0.3317645234741476,
          0.3317645234741476,
          0.0189317782804059,
          0.0189317782804059,
          0.6493036982454464,
          0.6493036982454464,
          0.0806961669858730,
          0.0806961669858730,
          0.0190301297436974,
          0.0190301297436974,
          0.9002737032704295,
          0.9002737032704295,
          0.3082449691963540,
          0.3082449691963540,
          0.1026061902393980,
          0.1026061902393981,
          0.5891488405642479,
          0.5891488405642479,
          0.1874417824837822,
          0.1874417824837821,
          0.0059363500168224,
          0.0059363500168222,
          0.8066218674993957,
          0.8066218674993957 };

  double w_save[] = {
          0.0452645660738188,
          0.0410929231436990,
          0.0410929231436990,
          0.0410929231436990,
          0.0407118333124254,
          0.0407118333124254,
          0.0407118333124254,
          0.0147816346902244,
          0.0147816346902244,
          0.0147816346902244,
          0.0294184096989881,
          0.0294184096989881,
          0.0294184096989881,
          0.0044185463121506,
          0.0044185463121506,
          0.0044185463121506,
          0.0259743332982772,
          0.0259743332982772,
          0.0259743332982772,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule17() returns the rule of precision 17.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.4171034443615992,
          0.1657931112768016,
          0.4171034443615992,
          0.0147554916607540,
          0.9704890166784921,
          0.0147554916607540,
          0.4655978716188903,
          0.0688042567622194,
          0.4655978716188903,
          0.1803581162663706,
          0.6392837674672588,
          0.1803581162663706,
          0.0666540634795969,
          0.8666918730408062,
          0.0666540634795969,
          0.2857065024365866,
          0.4285869951268267,
          0.2857065024365866,
          0.0160176423621193,
          0.8247900701650881,
          0.1591922874727927,
          0.8247900701650881,
          0.1591922874727927,
          0.0160176423621193,
          0.3062815917461865,
          0.6263690303864522,
          0.0673493778673612,
          0.6263690303864522,
          0.0673493778673612,
          0.3062815917461865,
          0.0132296727600869,
          0.5712948679446841,
          0.4154754592952291,
          0.5712948679446841,
          0.4154754592952291,
          0.0132296727600869,
          0.0780423405682824,
          0.7532351459364581,
          0.1687225134952595,
          0.7532351459364581,
          0.1687225134952595,
          0.0780423405682824,
          0.0131358708340027,
          0.7150722591106424,
          0.2717918700553549,
          0.7150722591106424,
          0.2717918700553549,
          0.0131358708340027,
          0.0115751759031806,
          0.9159193532978169,
          0.0725054707990024,
          0.9159193532978169,
          0.0725054707990024,
          0.0115751759031806,
          0.1575054779268699,
          0.5432755795961598,
          0.2992189424769703,
          0.5432755795961598,
          0.2992189424769703,
          0.1575054779268699 };

  static double b_save[] = {
          0.1657931112768016,
          0.4171034443615992,
          0.4171034443615992,
          0.9704890166784921,
          0.0147554916607540,
          0.0147554916607540,
          0.0688042567622194,
          0.4655978716188903,
          0.4655978716188903,
          0.6392837674672588,
          0.1803581162663706,
          0.1803581162663706,
          0.8666918730408062,
          0.0666540634795969,
          0.0666540634795969,
          0.4285869951268267,
          0.2857065024365866,
          0.2857065024365866,
          0.8247900701650881,
          0.0160176423621193,
          0.8247900701650881,
          0.1591922874727927,
          0.0160176423621193,
          0.1591922874727927,
          0.6263690303864522,
          0.3062815917461865,
          0.6263690303864522,
          0.0673493778673612,
          0.3062815917461865,
          0.0673493778673612,
          0.5712948679446841,
          0.0132296727600869,
          0.5712948679446841,
          0.4154754592952291,
          0.0132296727600869,
          0.4154754592952291,
          0.7532351459364581,
          0.0780423405682824,
          0.7532351459364581,
          0.1687225134952595,
          0.0780423405682824,
          0.1687225134952595,
          0.7150722591106424,
          0.0131358708340027,
          0.7150722591106424,
          0.2717918700553549,
          0.0131358708340027,
          0.2717918700553549,
          0.9159193532978169,
          0.0115751759031806,
          0.9159193532978169,
          0.0725054707990024,
          0.0115751759031806,
          0.0725054707990024,
          0.5432755795961598,
          0.1575054779268699,
          0.5432755795961598,
          0.2992189424769703,
          0.1575054779268699,
          0.2992189424769703 };

  static double c_save[] = {
          0.4171034443615993,
          0.4171034443615992,
          0.1657931112768016,
          0.0147554916607540,
          0.0147554916607540,
          0.9704890166784921,
          0.4655978716188902,
          0.4655978716188903,
          0.0688042567622194,
          0.1803581162663707,
          0.1803581162663705,
          0.6392837674672588,
          0.0666540634795969,
          0.0666540634795969,
          0.8666918730408062,
          0.2857065024365866,
          0.2857065024365866,
          0.4285869951268267,
          0.1591922874727927,
          0.1591922874727927,
          0.0160176423621192,
          0.0160176423621192,
          0.8247900701650881,
          0.8247900701650881,
          0.0673493778673613,
          0.0673493778673613,
          0.3062815917461865,
          0.3062815917461866,
          0.6263690303864523,
          0.6263690303864523,
          0.4154754592952290,
          0.4154754592952290,
          0.0132296727600869,
          0.0132296727600869,
          0.5712948679446841,
          0.5712948679446841,
          0.1687225134952595,
          0.1687225134952595,
          0.0780423405682824,
          0.0780423405682824,
          0.7532351459364581,
          0.7532351459364581,
          0.2717918700553549,
          0.2717918700553549,
          0.0131358708340027,
          0.0131358708340027,
          0.7150722591106424,
          0.7150722591106424,
          0.0725054707990024,
          0.0725054707990025,
          0.0115751759031806,
          0.0115751759031806,
          0.9159193532978169,
          0.9159193532978169,
          0.2992189424769703,
          0.2992189424769702,
          0.1575054779268699,
          0.1575054779268699,
          0.5432755795961597,
          0.5432755795961598 };

  double w_save[] = {
          0.0273109265281021,
          0.0273109265281021,
          0.0273109265281021,
          0.0027738875776376,
          0.0027738875776376,
          0.0027738875776376,
          0.0250194509504974,
          0.0250194509504974,
          0.0250194509504974,
          0.0263126305880180,
          0.0263126305880180,
          0.0263126305880180,
          0.0124590008023054,
          0.0124590008023054,
          0.0124590008023054,
          0.0377162371527953,
          0.0377162371527953,
          0.0377162371527953,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule18() returns the rule of precision 18.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.3999556280675762,
          0.2000887438648475,
          0.3999556280675762,
          0.4875803015748695,
          0.0248393968502609,
          0.4875803015748695,
          0.4618095064064492,
          0.0763809871871015,
          0.4618095064064492,
          0.2422647025142720,
          0.5154705949714561,
          0.2422647025142720,
          0.0388302560886856,
          0.9223394878226288,
          0.0388302560886856,
          0.0919477421216432,
          0.8161045157567136,
          0.0919477421216432,
          0.0458049158598608,
          0.7703723762146752,
          0.1838227079254640,
          0.7703723762146752,
          0.1838227079254640,
          0.0458049158598608,
          0.2063492574338379,
          0.6709539851942345,
          0.1226967573719275,
          0.6709539851942345,
          0.1226967573719275,
          0.2063492574338379,
          0.0038976110334734,
          0.6004189546342569,
          0.3956834343322697,
          0.6004189546342569,
          0.3956834343322697,
          0.0038976110334734,
          0.0134620167414450,
          0.8783421894675217,
          0.1081957937910333,
          0.8783421894675217,
          0.1081957937910333,
          0.0134620167414450,
          0.0402602834699081,
          0.6399880920047146,
          0.3197516245253773,
          0.6399880920047146,
          0.3197516245253773,
          0.0402602834699081,
          0.0052983351866098,
          0.7589294798551984,
          0.2357721849581917,
          0.7589294798551984,
          0.2357721849581917,
          0.0052983351866098,
          0.0005483600420423,
          0.9723607289627957,
          0.0270909109951620,
          0.9723607289627957,
          0.0270909109951620,
          0.0005483600420423,
          0.1205876951639246,
          0.5459187753861946,
          0.3334935294498808,
          0.5459187753861946,
          0.3334935294498808,
          0.1205876951639246 };

  static double b_save[] = {
          0.3333333333333334,
          0.2000887438648475,
          0.3999556280675762,
          0.3999556280675762,
          0.0248393968502609,
          0.4875803015748695,
          0.4875803015748695,
          0.0763809871871015,
          0.4618095064064492,
          0.4618095064064492,
          0.5154705949714561,
          0.2422647025142720,
          0.2422647025142720,
          0.9223394878226288,
          0.0388302560886856,
          0.0388302560886856,
          0.8161045157567136,
          0.0919477421216432,
          0.0919477421216432,
          0.7703723762146752,
          0.0458049158598608,
          0.7703723762146752,
          0.1838227079254640,
          0.0458049158598608,
          0.1838227079254640,
          0.6709539851942345,
          0.2063492574338379,
          0.6709539851942345,
          0.1226967573719275,
          0.2063492574338379,
          0.1226967573719275,
          0.6004189546342569,
          0.0038976110334734,
          0.6004189546342569,
          0.3956834343322697,
          0.0038976110334734,
          0.3956834343322697,
          0.8783421894675217,
          0.0134620167414450,
          0.8783421894675217,
          0.1081957937910333,
          0.0134620167414450,
          0.1081957937910333,
          0.6399880920047146,
          0.0402602834699081,
          0.6399880920047146,
          0.3197516245253773,
          0.0402602834699081,
          0.3197516245253773,
          0.7589294798551984,
          0.0052983351866098,
          0.7589294798551984,
          0.2357721849581917,
          0.0052983351866098,
          0.2357721849581917,
          0.9723607289627957,
          0.0005483600420423,
          0.9723607289627957,
          0.0270909109951620,
          0.0005483600420423,
          0.0270909109951620,
          0.5459187753861946,
          0.1205876951639246,
          0.5459187753861946,
          0.3334935294498808,
          0.1205876951639246,
          0.3334935294498808 };

  static double c_save[] = {
          0.3333333333333333,
          0.3999556280675762,
          0.3999556280675762,
          0.2000887438648475,
          0.4875803015748696,
          0.4875803015748695,
          0.0248393968502609,
          0.4618095064064492,
          0.4618095064064492,
          0.0763809871871015,
          0.2422647025142719,
          0.2422647025142719,
          0.5154705949714561,
          0.0388302560886855,
          0.0388302560886856,
          0.9223394878226288,
          0.0919477421216433,
          0.0919477421216433,
          0.8161045157567136,
          0.1838227079254640,
          0.1838227079254640,
          0.0458049158598608,
          0.0458049158598608,
          0.7703723762146752,
          0.7703723762146752,
          0.1226967573719275,
          0.1226967573719275,
          0.2063492574338379,
          0.2063492574338379,
          0.6709539851942345,
          0.6709539851942345,
          0.3956834343322697,
          0.3956834343322697,
          0.0038976110334734,
          0.0038976110334734,
          0.6004189546342569,
          0.6004189546342569,
          0.1081957937910333,
          0.1081957937910333,
          0.0134620167414450,
          0.0134620167414450,
          0.8783421894675217,
          0.8783421894675217,
          0.3197516245253773,
          0.3197516245253773,
          0.0402602834699081,
          0.0402602834699081,
          0.6399880920047146,
          0.6399880920047146,
          0.2357721849581917,
          0.2357721849581917,
          0.0052983351866098,
          0.0052983351866098,
          0.7589294798551984,
          0.7589294798551984,
          0.0270909109951620,
          0.0270909109951620,
          0.0005483600420423,
          0.0005483600420423,
          0.9723607289627956,
          0.9723607289627956,
          0.3334935294498808,
          0.3334935294498808,
          0.1205876951639246,
          0.1205876951639246,
          0.5459187753861946,
          0.5459187753861946 };

  double w_save[] = {
          0.0363557353014267,
          0.0333044700333901,
          0.0333044700333901,
          0.0333044700333901,
          0.0120466476339997,
          0.0120466476339997,
          0.0120466476339997,
          0.0189491715067789,
          0.0189491715067789,
          0.0189491715067789,
          0.0364750894089436,
          0.0364750894089436,
          0.0364750894089436,
          0.0071293260197190,
          0.0071293260197190,
          0.0071293260197190,
          0.0165591599520032,
          0.0165591599520032,
          0.0165591599520032,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule19() returns the rule of precision 19.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.0525238903512090,
          0.8949522192975821,
          0.0525238903512090,
          0.4925126750413369,
          0.0149746499173263,
          0.4925126750413369,
          0.1114488733230214,
          0.7771022533539573,
          0.1114488733230214,
          0.4591942010395437,
          0.0816115979209127,
          0.4591942010395437,
          0.4039697225519012,
          0.1920605548961976,
          0.4039697225519012,
          0.1781701047817643,
          0.6436597904364714,
          0.1781701047817643,
          0.0116394611837894,
          0.9767210776324211,
          0.0116394611837894,
          0.2551616329136077,
          0.4896767341727846,
          0.2551616329136077,
          0.1306976762680324,
          0.8301564644002754,
          0.0391458593316922,
          0.8301564644002754,
          0.0391458593316922,
          0.1306976762680324,
          0.3113176298095413,
          0.5593698057203009,
          0.1293125644701578,
          0.5593698057203009,
          0.1293125644701578,
          0.3113176298095413,
          0.0020689258966048,
          0.6333132931287841,
          0.3646177809746111,
          0.6333132931287841,
          0.3646177809746111,
          0.0020689258966048,
          0.0745602946016267,
          0.7040048199660421,
          0.2214348854323312,
          0.7040048199660421,
          0.2214348854323312,
          0.0745602946016267,
          0.0050072882573545,
          0.8525669543768892,
          0.1424257573657563,
          0.8525669543768892,
          0.1424257573657563,
          0.0050072882573545,
          0.0408880111960169,
          0.6050839790687079,
          0.3540280097352752,
          0.6050839790687079,
          0.3540280097352752,
          0.0408880111960169,
          0.2418945789605796,
          0.7431813689574364,
          0.0149240520819841,
          0.7431813689574364,
          0.0149240520819841,
          0.2418945789605796,
          0.0600862753223067,
          0.9301376988768051,
          0.0097760258008882,
          0.9301376988768051,
          0.0097760258008882,
          0.0600862753223067 };

  static double b_save[] = {
          0.3333333333333334,
          0.8949522192975821,
          0.0525238903512090,
          0.0525238903512090,
          0.0149746499173263,
          0.4925126750413369,
          0.4925126750413369,
          0.7771022533539573,
          0.1114488733230214,
          0.1114488733230214,
          0.0816115979209127,
          0.4591942010395437,
          0.4591942010395437,
          0.1920605548961976,
          0.4039697225519012,
          0.4039697225519012,
          0.6436597904364714,
          0.1781701047817643,
          0.1781701047817643,
          0.9767210776324211,
          0.0116394611837894,
          0.0116394611837894,
          0.4896767341727846,
          0.2551616329136077,
          0.2551616329136077,
          0.8301564644002754,
          0.1306976762680324,
          0.8301564644002754,
          0.0391458593316922,
          0.1306976762680324,
          0.0391458593316922,
          0.5593698057203009,
          0.3113176298095413,
          0.5593698057203009,
          0.1293125644701578,
          0.3113176298095413,
          0.1293125644701578,
          0.6333132931287841,
          0.0020689258966048,
          0.6333132931287841,
          0.3646177809746111,
          0.0020689258966048,
          0.3646177809746111,
          0.7040048199660421,
          0.0745602946016267,
          0.7040048199660421,
          0.2214348854323312,
          0.0745602946016267,
          0.2214348854323312,
          0.8525669543768892,
          0.0050072882573545,
          0.8525669543768892,
          0.1424257573657563,
          0.0050072882573545,
          0.1424257573657563,
          0.6050839790687079,
          0.0408880111960169,
          0.6050839790687079,
          0.3540280097352752,
          0.0408880111960169,
          0.3540280097352752,
          0.7431813689574364,
          0.2418945789605796,
          0.7431813689574364,
          0.0149240520819841,
          0.2418945789605796,
          0.0149240520819841,
          0.9301376988768051,
          0.0600862753223067,
          0.9301376988768051,
          0.0097760258008882,
          0.0600862753223067,
          0.0097760258008882 };

  static double c_save[] = {
          0.3333333333333333,
          0.0525238903512090,
          0.0525238903512090,
          0.8949522192975821,
          0.4925126750413369,
          0.4925126750413369,
          0.0149746499173262,
          0.1114488733230213,
          0.1114488733230212,
          0.7771022533539572,
          0.4591942010395437,
          0.4591942010395437,
          0.0816115979209127,
          0.4039697225519012,
          0.4039697225519012,
          0.1920605548961976,
          0.1781701047817643,
          0.1781701047817643,
          0.6436597904364714,
          0.0116394611837894,
          0.0116394611837894,
          0.9767210776324211,
          0.2551616329136077,
          0.2551616329136077,
          0.4896767341727846,
          0.0391458593316922,
          0.0391458593316922,
          0.1306976762680324,
          0.1306976762680324,
          0.8301564644002754,
          0.8301564644002754,
          0.1293125644701578,
          0.1293125644701578,
          0.3113176298095413,
          0.3113176298095413,
          0.5593698057203009,
          0.5593698057203009,
          0.3646177809746111,
          0.3646177809746111,
          0.0020689258966048,
          0.0020689258966048,
          0.6333132931287841,
          0.6333132931287841,
          0.2214348854323313,
          0.2214348854323311,
          0.0745602946016267,
          0.0745602946016266,
          0.7040048199660421,
          0.7040048199660421,
          0.1424257573657564,
          0.1424257573657564,
          0.0050072882573545,
          0.0050072882573544,
          0.8525669543768892,
          0.8525669543768892,
          0.3540280097352752,
          0.3540280097352753,
          0.0408880111960168,
          0.0408880111960169,
          0.6050839790687079,
          0.6050839790687079,
          0.0149240520819840,
          0.0149240520819840,
          0.2418945789605795,
          0.2418945789605795,
          0.7431813689574364,
          0.7431813689574364,
          0.0097760258008881,
          0.0097760258008882,
          0.0600862753223067,
          0.0600862753223068,
          0.9301376988768051,
          0.9301376988768051 };

  double w_save[] = {
          0.0344693977040123,
          0.0071092565977981,
          0.0071092565977981,
          0.0071092565977981,
          0.0103217551429443,
          0.0103217551429443,
          0.0103217551429443,
          0.0152343510930183,
          0.0152343510930183,
          0.0152343510930183,
          0.0229835900267416,
          0.0229835900267416,
          0.0229835900267416,
          0.0315375348931550,
          0.0315375348931550,
          0.0315375348931550,
          0.0246519148481909,
          0.0246519148481909,
          0.0246519148481909,
          0.0017653227764428,
          0.0017653227764428,
          0.0017653227764428,
          0.0317530193660031,
          0.0317530193660031,
          0.0317530193660031,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//   rule20() returns the rule of precision 20.
//
//   The data is given for the following triangle:
//
//    (0,1) 
//      |  *
//      |    * 
//      |     * 
//      |      * 
//    (0,0)----(1,0)
//
//   We suppose we are given a triangle T with vertices A, B, C.
//   We call a rule with n points, returning barycentric coordinates
//   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
//   approximated by Q as follows:
//
//   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
//   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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]: the barycentric coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double a_save[] = {
          0.3333333333333334,
          0.2545792676733391,
          0.4908414646533218,
          0.2545792676733391,
          0.0109761410283978,
          0.9780477179432046,
          0.0109761410283978,
          0.1093835967117146,
          0.7812328065765708,
          0.1093835967117146,
          0.1862949977445409,
          0.6274100045109181,
          0.1862949977445409,
          0.4455510569559248,
          0.1088978860881504,
          0.4455510569559248,
          0.0373108805988847,
          0.9253782388022307,
          0.0373108805988847,
          0.3934253478170999,
          0.2131493043658003,
          0.3934253478170999,
          0.4762456115404990,
          0.0475087769190020,
          0.4762456115404990,
          0.0075707805046965,
          0.8332955118382362,
          0.1591337076570672,
          0.8332955118382362,
          0.1591337076570672,
          0.0075707805046965,
          0.0465603649076643,
          0.7549215028635474,
          0.1985181322287882,
          0.7549215028635474,
          0.1985181322287882,
          0.0465603649076643,
          0.0640905856084341,
          0.9310544767839422,
          0.0048549376076237,
          0.9310544767839422,
          0.0048549376076237,
          0.0640905856084341,
          0.0549874791429868,
          0.6118777035474257,
          0.3331348173095875,
          0.6118777035474257,
          0.3331348173095875,
          0.0549874791429868,
          0.0999522962881387,
          0.8616840189364867,
          0.0383636847753746,
          0.8616840189364867,
          0.0383636847753746,
          0.0999522962881387,
          0.1062272047202700,
          0.6781657378896355,
          0.2156070573900944,
          0.6781657378896355,
          0.2156070573900944,
          0.1062272047202700,
          0.4200237588162241,
          0.5701446928909734,
          0.0098315482928026,
          0.5701446928909734,
          0.0098315482928026,
          0.4200237588162241,
          0.3178601238357720,
          0.5423318041724281,
          0.1398080719917999,
          0.5423318041724281,
          0.1398080719917999,
          0.3178601238357720,
          0.0107372128560111,
          0.7086813757203236,
          0.2805814114236652,
          0.7086813757203236,
          0.2805814114236652,
          0.0107372128560111 };

  static double b_save[] = {
          0.3333333333333334,
          0.4908414646533218,
          0.2545792676733391,
          0.2545792676733391,
          0.9780477179432046,
          0.0109761410283978,
          0.0109761410283978,
          0.7812328065765708,
          0.1093835967117146,
          0.1093835967117146,
          0.6274100045109181,
          0.1862949977445409,
          0.1862949977445409,
          0.1088978860881504,
          0.4455510569559248,
          0.4455510569559248,
          0.9253782388022307,
          0.0373108805988847,
          0.0373108805988847,
          0.2131493043658003,
          0.3934253478170999,
          0.3934253478170999,
          0.0475087769190020,
          0.4762456115404990,
          0.4762456115404990,
          0.8332955118382362,
          0.0075707805046965,
          0.8332955118382362,
          0.1591337076570672,
          0.0075707805046965,
          0.1591337076570672,
          0.7549215028635474,
          0.0465603649076643,
          0.7549215028635474,
          0.1985181322287882,
          0.0465603649076643,
          0.1985181322287882,
          0.9310544767839422,
          0.0640905856084341,
          0.9310544767839422,
          0.0048549376076237,
          0.0640905856084341,
          0.0048549376076237,
          0.6118777035474257,
          0.0549874791429868,
          0.6118777035474257,
          0.3331348173095875,
          0.0549874791429868,
          0.3331348173095875,
          0.8616840189364867,
          0.0999522962881387,
          0.8616840189364867,
          0.0383636847753746,
          0.0999522962881387,
          0.0383636847753746,
          0.6781657378896355,
          0.1062272047202700,
          0.6781657378896355,
          0.2156070573900944,
          0.1062272047202700,
          0.2156070573900944,
          0.5701446928909734,
          0.4200237588162241,
          0.5701446928909734,
          0.0098315482928026,
          0.4200237588162241,
          0.0098315482928026,
          0.5423318041724281,
          0.3178601238357720,
          0.5423318041724281,
          0.1398080719917999,
          0.3178601238357720,
          0.1398080719917999,
          0.7086813757203236,
          0.0107372128560111,
          0.7086813757203236,
          0.2805814114236652,
          0.0107372128560111,
          0.2805814114236652 };

  static double c_save[] = {
          0.3333333333333333,
          0.2545792676733392,
          0.2545792676733392,
          0.4908414646533218,
          0.0109761410283977,
          0.0109761410283977,
          0.9780477179432044,
          0.1093835967117146,
          0.1093835967117146,
          0.7812328065765708,
          0.1862949977445409,
          0.1862949977445409,
          0.6274100045109181,
          0.4455510569559249,
          0.4455510569559248,
          0.1088978860881504,
          0.0373108805988847,
          0.0373108805988847,
          0.9253782388022306,
          0.3934253478170999,
          0.3934253478170999,
          0.2131493043658003,
          0.4762456115404990,
          0.4762456115404990,
          0.0475087769190020,
          0.1591337076570672,
          0.1591337076570672,
          0.0075707805046965,
          0.0075707805046965,
          0.8332955118382362,
          0.8332955118382362,
          0.1985181322287882,
          0.1985181322287883,
          0.0465603649076644,
          0.0465603649076645,
          0.7549215028635475,
          0.7549215028635475,
          0.0048549376076237,
          0.0048549376076238,
          0.0640905856084341,
          0.0640905856084342,
          0.9310544767839422,
          0.9310544767839422,
          0.3331348173095875,
          0.3331348173095876,
          0.0549874791429867,
          0.0549874791429869,
          0.6118777035474257,
          0.6118777035474257,
          0.0383636847753746,
          0.0383636847753746,
          0.0999522962881386,
          0.0999522962881387,
          0.8616840189364867,
          0.8616840189364867,
          0.2156070573900944,
          0.2156070573900944,
          0.1062272047202700,
          0.1062272047202701,
          0.6781657378896355,
          0.6781657378896355,
          0.0098315482928025,
          0.0098315482928025,
          0.4200237588162241,
          0.4200237588162241,
          0.5701446928909734,
          0.5701446928909732,
          0.1398080719917999,
          0.1398080719917999,
          0.3178601238357720,
          0.3178601238357720,
          0.5423318041724281,
          0.5423318041724281,
          0.2805814114236652,
          0.2805814114236653,
          0.0107372128560111,
          0.0107372128560111,
          0.7086813757203236,
          0.7086813757203236 };

  double w_save[] = {
          0.0278202214029062,
          0.0281664026150405,
          0.0281664026150405,
          0.0281664026150405,
          0.0015976815821332,
          0.0015976815821332,
          0.0015976815821332,
          0.0156604615521491,
          0.0156604615521491,
          0.0156604615521491,
          0.0183469259485058,
          0.0183469259485058,
          0.0183469259485058,
          0.0189047998664649,
          0.0189047998664649,
          0.0189047998664649,
          0.0043225508213312,
          0.0043225508213312,
          0.0043225508213312,
          0.0275761012581409,
          0.0275761012581409,
          0.0275761012581409,
          0.0142036506068169,
          0.0142036506068169,
          0.0142036506068169,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154 };

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

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

double triangle_area ( double vert1[], double vert2[], double vert3[] )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_area() returns the area of a triangle.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    26 June 2014
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    double vert1[2], vert2[2], vert3[2]: the vertices of the triangle.
//
//  Output:
//
//    double triangle_area: the area of the triangle.
//
{
  double value;

  value = 0.5 * 
    ( 
        ( vert2[0] - vert1[0] ) * ( vert3[1] - vert1[1] ) 
      - ( vert3[0] - vert1[0] ) * ( vert2[1] - vert1[1] ) 
    );

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

double triangle_unit_area ( )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_unit_area() returns the area of a unit triangle.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    24 May 2023
//
//  Author:
//
//    John Burkardt.
//
//  Output:
//
//    double triangle_unit_area: the area of the triangle.
//
{
  double value;

  value = 0.5;

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

double triangle_unit_monomial_integral ( int expon[] )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_unit_monomial_integral() integrates a monomial over the unit triangle.
//
//  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 triangle ) x^m y^n dx dy = m! * n! / ( m + n + 2 )!
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    04 July 2007
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int EXPON[3], the exponents.
//
//  Output:
//
//    double triangle_unit_monomial_integral: the value of the integral 
//    of the monomial.
//
{
  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 );
  }

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

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

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

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

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

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

  if ( p == 0 )
  {
    rule00 ( n, a, b, c, w );
  }
  else if ( p == 1 )
  {
    rule01 ( n, a, b, c, w );
  }
  else if ( p == 2 )
  {
    rule02 ( n, a, b, c, w );
  }
  else if ( p == 3 )
  {
    rule03 ( n, a, b, c, w );
  }
  else if ( p == 4 )
  {
    rule04 ( n, a, b, c, w );
  }
  else if ( p == 5 )
  {
    rule05 ( n, a, b, c, w );
  }
  else if ( p == 6 )
  {
    rule06 ( n, a, b, c, w );
  }
  else if ( p == 7 )
  {
    rule07 ( n, a, b, c, w );
  }
  else if ( p == 8 )
  {
    rule08 ( n, a, b, c, w );
  }
  else if ( p == 9 )
  {
    rule09 ( n, a, b, c, w );
  }
  else if ( p == 10 )
  {
    rule10 ( n, a, b, c, w );
  }
  else if ( p == 11 )
  {
    rule11 ( n, a, b, c, w );
  }
  else if ( p == 12 )
  {
    rule12 ( n, a, b, c, w );
  }
  else if ( p == 13 )
  {
    rule13 ( n, a, b, c, w );
  }
  else if ( p == 14 )
  {
    rule14 ( n, a, b, c, w );
  }
  else if ( p == 15 )
  {
    rule15 ( n, a, b, c, w );
  }
  else if ( p == 16 )
  {
    rule16 ( n, a, b, c, w );
  }
  else if ( p == 17 )
  {
    rule17 ( n, a, b, c, w );
  }
  else if ( p == 18 )
  {
    rule18 ( n, a, b, c, w );
  }
  else if ( p == 19 )
  {
    rule19 ( n, a, b, c, w );
  }
  else if ( p == 20 )
  {
    rule20 ( n, a, b, c, w );
  }

  return;
}