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

using namespace std;

# include "quadrilateral_witherden_rule.hpp"

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

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

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

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

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

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

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

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

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

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

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

double quadrilateral_unit_area ( )

//****************************************************************************80
//
//  Purpose:
//
//    quadrilateral_unit_area() returns the area of a unit quadrilateral.
//
//  Discussion:
//
//    The unit quadrilateral has vertices (0,0), (1,0), (1,1), (0,1).
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    real value: the area.
//
{
  double value;

  value = 1.0;

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

double quadrilateral_unit_monomial_integral ( int expon[3] )

//****************************************************************************80
//
//  Purpose:
//
//    quadrilateral_unit_monomial_integral(): monomial integral in a unit quadrilateral.
//
//  Discussion:
//
//    This function returns the value of the integral of X^ALPHA Y^BETA
//    over the unit quadrilateral.
//
//    The unit quadrilateral has vertices 
//    (0,0), (1,0), (1,1), (0,1).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int EXPON[2], the exponents.
//
//  Output:
//
//    double quadrilateral_unit_monomial_integral, the integral of the monomial.
//
{
  double value;

  value = 1.0 / ( expon[0] + 1 ) / ( expon[1] + 1 );

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

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

//****************************************************************************80
//
//  Purpose:
//
//    quadrilateral_witherden_rule() returns a quadrilateral quadrature rule of given precision.
//
//  Discussion:
//
//    The unit quadrilateral has vertices 
//    (0,0), (1,0), (1,1), (0,1).
//
//  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:
//
//    int p: the precision, 0 <= p <= 21.
//
//    int n: the order of the rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w(n): the quadrature weights.
//
{
  if ( p < 0 ) 
  {
    cout << "\n";
    cout << "quadrilateral_witherden_rule(): Fatal error!\n";
    cout << "  Input p < 0.\n";
    exit ( 1 );
  }

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

  if ( p <= 1 )
  {
    rule01 ( n, x, y, w );
  }
  else if ( p <= 3 )
  {
    rule03 ( n, x, y, w );
  }
  else if ( p <= 5 )
  {
    rule05 ( n, x, y, w );
  }
  else if ( p <= 7 )
  {
    rule07 ( n, x, y, w );
  }
  else if ( p <= 9 )
  {
    rule09 ( n, x, y, w );
  }
  else if ( p <= 11 )
  {
    rule11 ( n, x, y, w );
  }
  else if ( p <= 13 )
  {
    rule13 ( n, x, y, w );
  }
  else if ( p <= 15 )
  {
    rule15 ( n, x, y, w );
  }
  else if ( p <= 17 )
  {
    rule17 ( n, x, y, w );
  }
  else if ( p <= 19 )
  {
    rule19 ( n, x, y, w );
  }
  else if ( p <= 21 )
  {
    rule21 ( n, x, y, w );
  }

  return;
}
//****************************************************************************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 prism 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:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int p: the precision, 0 <= p <= 21.
//
//  Output:
//
//    int rule_order: the order of the rule.
//
{
  int order;
  static int order_save[] = {
       1,
       1,  4,  4,  8,  8, 12, 12, 20, 20, 28,
      28, 37, 37, 48, 48, 60, 60, 72, 72, 85,
      85 };

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

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

  order = order_save[p];

  return order;
}

void rule01 ( int n, double x[], double y[], double w[] )

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

  static double y_save[] = {
          0.5000000000000000 };

  double w_save[] = {
          1.0000000000000000 };

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

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

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

  static double y_save[] = {
          0.7886751345948129,
          0.2113248654051871,
          0.7886751345948129,
          0.2113248654051871 };

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

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

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

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

  static double y_save[] = {
          0.5000000000000000,
          0.8415650255319866,
          0.5000000000000000,
          0.1584349744680134,
          0.9409585518440984,
          0.0590414481559016,
          0.9409585518440984,
          0.0590414481559016 };

  double w_save[] = {
          0.2040816326530612,
          0.2040816326530612,
          0.2040816326530612,
          0.2040816326530612,
          0.0459183673469388,
          0.0459183673469388,
          0.0459183673469388,
          0.0459183673469388 };

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

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

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

  static double y_save[] = {
          0.5000000000000000,
          0.9629100498862757,
          0.5000000000000000,
          0.0370899501137243,
          0.9029898914592993,
          0.0970101085407006,
          0.9029898914592993,
          0.0970101085407006,
          0.6902772166041579,
          0.3097227833958421,
          0.6902772166041579,
          0.3097227833958421 };

  double w_save[] = {
          0.0604938271604938,
          0.0604938271604938,
          0.0604938271604938,
          0.0604938271604938,
          0.0593579436726576,
          0.0593579436726576,
          0.0593579436726576,
          0.0593579436726576,
          0.1301482291668486,
          0.1301482291668486,
          0.1301482291668486,
          0.1301482291668486 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule09() returns the quadrilateral rule of precision 9.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.7444634284871845,
          0.5000000000000000,
          0.2555365715128155,
          0.5000000000000000,
          0.9698276290484189,
          0.9698276290484189,
          0.0301723709515812,
          0.0301723709515812,
          0.8454402752431720,
          0.8454402752431720,
          0.1545597247568281,
          0.1545597247568281,
          0.9593102205283611,
          0.6724360126822018,
          0.9593102205283611,
          0.3275639873177982,
          0.0406897794716389,
          0.6724360126822018,
          0.0406897794716389,
          0.3275639873177982 };

  static double y_save[] = {
          0.5000000000000000,
          0.7444634284871845,
          0.5000000000000000,
          0.2555365715128155,
          0.9698276290484189,
          0.0301723709515812,
          0.9698276290484189,
          0.0301723709515812,
          0.8454402752431720,
          0.1545597247568281,
          0.8454402752431720,
          0.1545597247568281,
          0.6724360126822018,
          0.9593102205283611,
          0.3275639873177982,
          0.9593102205283611,
          0.6724360126822018,
          0.0406897794716389,
          0.3275639873177982,
          0.0406897794716389 };

  double w_save[] = {
          0.1135409901716873,
          0.1135409901716873,
          0.1135409901716873,
          0.1135409901716873,
          0.0106828079664439,
          0.0106828079664439,
          0.0106828079664439,
          0.0106828079664439,
          0.0535500902317154,
          0.0535500902317154,
          0.0535500902317154,
          0.0535500902317154,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767,
          0.0361130558150767 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule11() returns the quadrilateral rule of precision 11.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.8573089148323030,
          0.5000000000000000,
          0.1426910851676970,
          0.5000000000000000,
          0.6368286050857298,
          0.6368286050857298,
          0.3631713949142702,
          0.3631713949142702,
          0.8183019661061506,
          0.8183019661061506,
          0.1816980338938495,
          0.1816980338938495,
          0.9758151943920168,
          0.9077827168448191,
          0.9758151943920168,
          0.0922172831551808,
          0.0241848056079833,
          0.9077827168448191,
          0.0241848056079833,
          0.0922172831551808,
          0.6731036000238227,
          0.9677839357437956,
          0.6731036000238227,
          0.0322160642562044,
          0.3268963999761773,
          0.9677839357437956,
          0.3268963999761773,
          0.0322160642562044 };

  static double y_save[] = {
          0.5000000000000000,
          0.8573089148323030,
          0.5000000000000000,
          0.1426910851676970,
          0.6368286050857298,
          0.3631713949142702,
          0.6368286050857298,
          0.3631713949142702,
          0.8183019661061506,
          0.1816980338938495,
          0.8183019661061506,
          0.1816980338938495,
          0.9077827168448191,
          0.9758151943920168,
          0.0922172831551808,
          0.9758151943920168,
          0.9077827168448191,
          0.0241848056079833,
          0.0922172831551808,
          0.0241848056079833,
          0.9677839357437956,
          0.6731036000238227,
          0.0322160642562044,
          0.6731036000238227,
          0.9677839357437956,
          0.3268963999761773,
          0.0322160642562044,
          0.3268963999761773 };

  double w_save[] = {
          0.0543501099671780,
          0.0543501099671780,
          0.0543501099671780,
          0.0543501099671780,
          0.0693185257459628,
          0.0693185257459628,
          0.0693185257459628,
          0.0693185257459628,
          0.0534834094695620,
          0.0534834094695620,
          0.0534834094695620,
          0.0534834094695620,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0110186422787458,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028,
          0.0254053351299028 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule13() returns the quadrilateral rule of precision 13.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.5000000000000000,
          0.9917306663066228,
          0.5000000000000000,
          0.0082693336933772,
          0.5000000000000000,
          0.8199307091835548,
          0.5000000000000000,
          0.1800692908164451,
          0.5000000000000000,
          0.9593892440398557,
          0.9593892440398557,
          0.0406107559601443,
          0.0406107559601443,
          0.6898142525933719,
          0.6898142525933719,
          0.3101857474066281,
          0.3101857474066281,
          0.8497771082566756,
          0.8497771082566756,
          0.1502228917433244,
          0.1502228917433244,
          0.8217986874983090,
          0.9875344419529918,
          0.8217986874983090,
          0.0124655580470082,
          0.1782013125016910,
          0.9875344419529918,
          0.1782013125016910,
          0.0124655580470082,
          0.6667699405823916,
          0.9322138046335309,
          0.6667699405823916,
          0.0677861953664691,
          0.3332300594176084,
          0.9322138046335309,
          0.3332300594176084,
          0.0677861953664691 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.9917306663066228,
          0.5000000000000000,
          0.0082693336933772,
          0.5000000000000000,
          0.8199307091835548,
          0.5000000000000000,
          0.1800692908164451,
          0.9593892440398557,
          0.0406107559601443,
          0.9593892440398557,
          0.0406107559601443,
          0.6898142525933719,
          0.3101857474066281,
          0.6898142525933719,
          0.3101857474066281,
          0.8497771082566756,
          0.1502228917433244,
          0.8497771082566756,
          0.1502228917433244,
          0.9875344419529918,
          0.8217986874983090,
          0.0124655580470082,
          0.8217986874983090,
          0.9875344419529918,
          0.1782013125016910,
          0.0124655580470082,
          0.1782013125016910,
          0.9322138046335309,
          0.6667699405823916,
          0.0677861953664691,
          0.6667699405823916,
          0.9322138046335309,
          0.3332300594176084,
          0.0677861953664691,
          0.3332300594176084 };

  double w_save[] = {
          0.0749983360897274,
          0.0095321563374632,
          0.0095321563374632,
          0.0095321563374632,
          0.0095321563374632,
          0.0461338672424518,
          0.0461338672424518,
          0.0461338672424518,
          0.0461338672424518,
          0.0098767851618631,
          0.0098767851618631,
          0.0098767851618631,
          0.0098767851618631,
          0.0578496298501713,
          0.0578496298501713,
          0.0578496298501713,
          0.0578496298501713,
          0.0343105241782509,
          0.0343105241782509,
          0.0343105241782509,
          0.0343105241782509,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0083799450125954,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886,
          0.0283937815910886 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule15() returns the quadrilateral rule of precision 15.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.8990494939485035,
          0.5000000000000000,
          0.1009505060514965,
          0.5000000000000000,
          0.6519265131992299,
          0.5000000000000000,
          0.3480734868007701,
          0.5000000000000000,
          0.9412111350534427,
          0.9412111350534427,
          0.0587888649465573,
          0.0587888649465573,
          0.9888948997699514,
          0.9888948997699514,
          0.0111051002300486,
          0.0111051002300486,
          0.9043560679097518,
          0.7836067866982954,
          0.9043560679097518,
          0.2163932133017046,
          0.0956439320902482,
          0.7836067866982954,
          0.0956439320902482,
          0.2163932133017046,
          0.6523273995185263,
          0.7893680970179033,
          0.6523273995185263,
          0.2106319029820967,
          0.3476726004814737,
          0.7893680970179033,
          0.3476726004814737,
          0.2106319029820967,
          0.9902524218622659,
          0.8487318159954835,
          0.9902524218622659,
          0.1512681840045165,
          0.0097475781377340,
          0.8487318159954835,
          0.0097475781377340,
          0.1512681840045165,
          0.6324220779361581,
          0.9778712599047559,
          0.6324220779361581,
          0.0221287400952441,
          0.3675779220638419,
          0.9778712599047559,
          0.3675779220638419,
          0.0221287400952441 };

  static double y_save[] = {
          0.5000000000000000,
          0.8990494939485035,
          0.5000000000000000,
          0.1009505060514965,
          0.5000000000000000,
          0.6519265131992299,
          0.5000000000000000,
          0.3480734868007701,
          0.9412111350534427,
          0.0587888649465573,
          0.9412111350534427,
          0.0587888649465573,
          0.9888948997699514,
          0.0111051002300486,
          0.9888948997699514,
          0.0111051002300486,
          0.7836067866982954,
          0.9043560679097518,
          0.2163932133017046,
          0.9043560679097518,
          0.7836067866982954,
          0.0956439320902482,
          0.2163932133017046,
          0.0956439320902482,
          0.7893680970179033,
          0.6523273995185263,
          0.2106319029820967,
          0.6523273995185263,
          0.7893680970179033,
          0.3476726004814737,
          0.2106319029820967,
          0.3476726004814737,
          0.8487318159954835,
          0.9902524218622659,
          0.1512681840045165,
          0.9902524218622659,
          0.8487318159954835,
          0.0097475781377340,
          0.1512681840045165,
          0.0097475781377340,
          0.9778712599047559,
          0.6324220779361581,
          0.0221287400952441,
          0.6324220779361581,
          0.9778712599047559,
          0.3675779220638419,
          0.0221287400952441,
          0.3675779220638419 };

  double w_save[] = {
          0.0287332318158813,
          0.0287332318158813,
          0.0287332318158813,
          0.0287332318158813,
          0.0454214397418180,
          0.0454214397418180,
          0.0454214397418180,
          0.0454214397418180,
          0.0103096209471915,
          0.0103096209471915,
          0.0103096209471915,
          0.0103096209471915,
          0.0014834366214598,
          0.0014834366214598,
          0.0014834366214598,
          0.0014834366214598,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0252872858693581,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0369591804082203,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0058182985361683,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780,
          0.0139613706230780 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule17() returns the quadrilateral rule of precision 17.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.7725214936045741,
          0.5000000000000000,
          0.2274785063954258,
          0.5000000000000000,
          0.9587485941095766,
          0.5000000000000000,
          0.0412514058904234,
          0.5000000000000000,
          0.5962771213070235,
          0.5962771213070235,
          0.4037228786929766,
          0.4037228786929766,
          0.9356643423477216,
          0.9356643423477216,
          0.0643356576522784,
          0.0643356576522784,
          0.8471940450447754,
          0.8471940450447754,
          0.1528059549552246,
          0.1528059549552246,
          0.7288414964415777,
          0.7288414964415777,
          0.2711585035584222,
          0.2711585035584222,
          0.9846021280457544,
          0.9846021280457544,
          0.0153978719542456,
          0.0153978719542456,
          0.8810441380400991,
          0.6397388833575961,
          0.8810441380400991,
          0.3602611166424039,
          0.1189558619599009,
          0.6397388833575961,
          0.1189558619599009,
          0.3602611166424039,
          0.9536971760991118,
          0.7734493352071550,
          0.9536971760991118,
          0.2265506647928450,
          0.0463028239008882,
          0.7734493352071550,
          0.0463028239008882,
          0.2265506647928450,
          0.8806243332195414,
          0.9918939857507747,
          0.8806243332195414,
          0.0081060142492252,
          0.1193756667804586,
          0.9918939857507747,
          0.1193756667804586,
          0.0081060142492252,
          0.9935473535223840,
          0.6514037408944306,
          0.9935473535223840,
          0.3485962591055693,
          0.0064526464776161,
          0.6514037408944306,
          0.0064526464776161,
          0.3485962591055693 };

  static double y_save[] = {
          0.5000000000000000,
          0.7725214936045741,
          0.5000000000000000,
          0.2274785063954258,
          0.5000000000000000,
          0.9587485941095766,
          0.5000000000000000,
          0.0412514058904234,
          0.5962771213070235,
          0.4037228786929766,
          0.5962771213070235,
          0.4037228786929766,
          0.9356643423477216,
          0.0643356576522784,
          0.9356643423477216,
          0.0643356576522784,
          0.8471940450447754,
          0.1528059549552246,
          0.8471940450447754,
          0.1528059549552246,
          0.7288414964415777,
          0.2711585035584222,
          0.7288414964415777,
          0.2711585035584222,
          0.9846021280457544,
          0.0153978719542456,
          0.9846021280457544,
          0.0153978719542456,
          0.6397388833575961,
          0.8810441380400991,
          0.3602611166424039,
          0.8810441380400991,
          0.6397388833575961,
          0.1189558619599009,
          0.3602611166424039,
          0.1189558619599009,
          0.7734493352071550,
          0.9536971760991118,
          0.2265506647928450,
          0.9536971760991118,
          0.7734493352071550,
          0.0463028239008882,
          0.2265506647928450,
          0.0463028239008882,
          0.9918939857507747,
          0.8806243332195414,
          0.0081060142492252,
          0.8806243332195414,
          0.9918939857507747,
          0.1193756667804586,
          0.0081060142492252,
          0.1193756667804586,
          0.6514037408944306,
          0.9935473535223840,
          0.3485962591055693,
          0.9935473535223840,
          0.6514037408944306,
          0.0064526464776161,
          0.3485962591055693,
          0.0064526464776161 };

  double w_save[] = {
          0.0347010944218050,
          0.0347010944218050,
          0.0347010944218050,
          0.0347010944218050,
          0.0163369559518662,
          0.0163369559518662,
          0.0163369559518662,
          0.0163369559518662,
          0.0357141049055499,
          0.0357141049055499,
          0.0357141049055499,
          0.0357141049055499,
          0.0091789771213129,
          0.0091789771213129,
          0.0091789771213129,
          0.0091789771213129,
          0.0215639249095360,
          0.0215639249095360,
          0.0215639249095360,
          0.0215639249095360,
          0.0336390683029792,
          0.0336390683029792,
          0.0336390683029792,
          0.0336390683029792,
          0.0019139083537275,
          0.0019139083537275,
          0.0019139083537275,
          0.0019139083537275,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0257941811962338,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0137066000593229,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0037799525970642,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908,
          0.0051952491639908 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule19() returns the quadrilateral rule of precision 19.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.8571721297863971,
          0.5000000000000000,
          0.1428278702136029,
          0.5000000000000000,
          0.6328360260604818,
          0.5000000000000000,
          0.3671639739395182,
          0.5000000000000000,
          0.9822171346289836,
          0.5000000000000000,
          0.0177828653710164,
          0.5000000000000000,
          0.9016898647338425,
          0.9016898647338425,
          0.0983101352661575,
          0.0983101352661575,
          0.9960827029485674,
          0.9960827029485674,
          0.0039172970514326,
          0.0039172970514326,
          0.9647248013998047,
          0.9647248013998047,
          0.0352751986001953,
          0.0352751986001953,
          0.7551391286846717,
          0.6333201572972811,
          0.7551391286846717,
          0.3666798427027189,
          0.2448608713153283,
          0.6333201572972811,
          0.2448608713153283,
          0.3666798427027189,
          0.6953671028876249,
          0.9939464665708677,
          0.6953671028876249,
          0.0060535334291323,
          0.3046328971123751,
          0.9939464665708677,
          0.3046328971123751,
          0.0060535334291323,
          0.8585839606548726,
          0.7562459386080489,
          0.8585839606548726,
          0.2437540613919512,
          0.1414160393451274,
          0.7562459386080489,
          0.1414160393451274,
          0.2437540613919512,
          0.6327200039056480,
          0.9344512170772521,
          0.6327200039056480,
          0.0655487829227479,
          0.3672799960943520,
          0.9344512170772521,
          0.3672799960943520,
          0.0655487829227479,
          0.8100176993466282,
          0.9631514779035646,
          0.8100176993466282,
          0.0368485220964354,
          0.1899823006533718,
          0.9631514779035646,
          0.1899823006533718,
          0.0368485220964354,
          0.9008357923592984,
          0.9942232653419869,
          0.9008357923592984,
          0.0057767346580131,
          0.0991642076407016,
          0.9942232653419869,
          0.0991642076407016,
          0.0057767346580131 };

  static double y_save[] = {
          0.5000000000000000,
          0.8571721297863971,
          0.5000000000000000,
          0.1428278702136029,
          0.5000000000000000,
          0.6328360260604818,
          0.5000000000000000,
          0.3671639739395182,
          0.5000000000000000,
          0.9822171346289836,
          0.5000000000000000,
          0.0177828653710164,
          0.9016898647338425,
          0.0983101352661575,
          0.9016898647338425,
          0.0983101352661575,
          0.9960827029485674,
          0.0039172970514326,
          0.9960827029485674,
          0.0039172970514326,
          0.9647248013998047,
          0.0352751986001953,
          0.9647248013998047,
          0.0352751986001953,
          0.6333201572972811,
          0.7551391286846717,
          0.3666798427027189,
          0.7551391286846717,
          0.6333201572972811,
          0.2448608713153283,
          0.3666798427027189,
          0.2448608713153283,
          0.9939464665708677,
          0.6953671028876249,
          0.0060535334291323,
          0.6953671028876249,
          0.9939464665708677,
          0.3046328971123751,
          0.0060535334291323,
          0.3046328971123751,
          0.7562459386080489,
          0.8585839606548726,
          0.2437540613919512,
          0.8585839606548726,
          0.7562459386080489,
          0.1414160393451274,
          0.2437540613919512,
          0.1414160393451274,
          0.9344512170772521,
          0.6327200039056480,
          0.0655487829227479,
          0.6327200039056480,
          0.9344512170772521,
          0.3672799960943520,
          0.0655487829227479,
          0.3672799960943520,
          0.9631514779035646,
          0.8100176993466282,
          0.0368485220964354,
          0.8100176993466282,
          0.9631514779035646,
          0.1899823006533718,
          0.0368485220964354,
          0.1899823006533718,
          0.9942232653419869,
          0.9008357923592984,
          0.0057767346580131,
          0.9008357923592984,
          0.9942232653419869,
          0.0991642076407016,
          0.0057767346580131,
          0.0991642076407016 };

  double w_save[] = {
          0.0244343411732072,
          0.0244343411732072,
          0.0244343411732072,
          0.0244343411732072,
          0.0348269032306206,
          0.0348269032306206,
          0.0348269032306206,
          0.0348269032306206,
          0.0087173958729720,
          0.0087173958729720,
          0.0087173958729720,
          0.0087173958729720,
          0.0119511367906989,
          0.0119511367906989,
          0.0119511367906989,
          0.0119511367906989,
          0.0004044287944779,
          0.0004044287944779,
          0.0004044287944779,
          0.0004044287944779,
          0.0043602620085889,
          0.0043602620085889,
          0.0043602620085889,
          0.0043602620085889,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0294314582100140,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0040448944029145,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0207116547458512,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0162472703731490,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0097451380991026,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859,
          0.0024723502336859 };

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rule21() returns the quadrilateral rule of precision 21.
//
//  Discussion:
//
//    We suppose we are given a quadrilateral Quad with vertices
//      (0,0), (1,0), (1,1), (0,1),
//    We call a rule with n points, returning coordinates
//    x, y, and weights w.  Then the integral I of f(x,y) over
//    Quad is approximated by Q as follows:
//
//    Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Freddie Witherden, Peter Vincent,
//    On the identification of symmetric quadrature rules for finite element methods,
//    Computers and Mathematics with Applications,
//    Volume 69, pages 1232-1241, 2015.
//
//  Input:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.5000000000000000,
          0.7366755371791120,
          0.5000000000000000,
          0.2633244628208880,
          0.5000000000000000,
          0.9176392148779342,
          0.5000000000000000,
          0.0823607851220658,
          0.5000000000000000,
          0.6286859503036145,
          0.6286859503036145,
          0.3713140496963855,
          0.3713140496963855,
          0.9816689160578117,
          0.9816689160578117,
          0.0183310839421883,
          0.0183310839421883,
          0.9312259626898258,
          0.9312259626898258,
          0.0687740373101742,
          0.0687740373101742,
          0.7484489812596729,
          0.7484489812596729,
          0.2515510187403271,
          0.2515510187403271,
          0.8521660875977003,
          0.8521660875977003,
          0.1478339124022997,
          0.1478339124022997,
          0.6209390927383510,
          0.8370731099776589,
          0.6209390927383510,
          0.1629268900223411,
          0.3790609072616490,
          0.8370731099776589,
          0.3790609072616490,
          0.1629268900223411,
          0.7400784831563976,
          0.9123236876354603,
          0.7400784831563976,
          0.0876763123645397,
          0.2599215168436024,
          0.9123236876354603,
          0.2599215168436024,
          0.0876763123645397,
          0.9661209679608770,
          0.6353495920508325,
          0.9661209679608770,
          0.3646504079491675,
          0.0338790320391230,
          0.6353495920508325,
          0.0338790320391230,
          0.3646504079491675,
          0.9674511629120053,
          0.8375197837185376,
          0.9674511629120053,
          0.1624802162814624,
          0.0325488370879947,
          0.8375197837185376,
          0.0325488370879947,
          0.1624802162814624,
          0.9952277337647621,
          0.7444581255885834,
          0.9952277337647621,
          0.2555418744114166,
          0.0047722662352379,
          0.7444581255885834,
          0.0047722662352379,
          0.2555418744114166,
          0.9155676229879507,
          0.9944803056932863,
          0.9155676229879507,
          0.0055196943067138,
          0.0844323770120493,
          0.9944803056932863,
          0.0844323770120493,
          0.0055196943067138,
          0.5457348004611456,
          0.9920085742447995,
          0.5457348004611456,
          0.0079914257552005,
          0.4542651995388544,
          0.9920085742447995,
          0.4542651995388544,
          0.0079914257552005 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.7366755371791120,
          0.5000000000000000,
          0.2633244628208880,
          0.5000000000000000,
          0.9176392148779342,
          0.5000000000000000,
          0.0823607851220658,
          0.6286859503036145,
          0.3713140496963855,
          0.6286859503036145,
          0.3713140496963855,
          0.9816689160578117,
          0.0183310839421883,
          0.9816689160578117,
          0.0183310839421883,
          0.9312259626898258,
          0.0687740373101742,
          0.9312259626898258,
          0.0687740373101742,
          0.7484489812596729,
          0.2515510187403271,
          0.7484489812596729,
          0.2515510187403271,
          0.8521660875977003,
          0.1478339124022997,
          0.8521660875977003,
          0.1478339124022997,
          0.8370731099776589,
          0.6209390927383510,
          0.1629268900223411,
          0.6209390927383510,
          0.8370731099776589,
          0.3790609072616490,
          0.1629268900223411,
          0.3790609072616490,
          0.9123236876354603,
          0.7400784831563976,
          0.0876763123645397,
          0.7400784831563976,
          0.9123236876354603,
          0.2599215168436024,
          0.0876763123645397,
          0.2599215168436024,
          0.6353495920508325,
          0.9661209679608770,
          0.3646504079491675,
          0.9661209679608770,
          0.6353495920508325,
          0.0338790320391230,
          0.3646504079491675,
          0.0338790320391230,
          0.8375197837185376,
          0.9674511629120053,
          0.1624802162814624,
          0.9674511629120053,
          0.8375197837185376,
          0.0325488370879947,
          0.1624802162814624,
          0.0325488370879947,
          0.7444581255885834,
          0.9952277337647621,
          0.2555418744114166,
          0.9952277337647621,
          0.7444581255885834,
          0.0047722662352379,
          0.2555418744114166,
          0.0047722662352379,
          0.9944803056932863,
          0.9155676229879507,
          0.0055196943067138,
          0.9155676229879507,
          0.9944803056932863,
          0.0844323770120493,
          0.0055196943067138,
          0.0844323770120493,
          0.9920085742447995,
          0.5457348004611456,
          0.0079914257552005,
          0.5457348004611456,
          0.9920085742447995,
          0.4542651995388544,
          0.0079914257552005,
          0.4542651995388544 };

  double w_save[] = {
          0.0337982100890357,
          0.0266985396964851,
          0.0266985396964851,
          0.0266985396964851,
          0.0266985396964851,
          0.0165643020012361,
          0.0165643020012361,
          0.0165643020012361,
          0.0165643020012361,
          0.0299711209115851,
          0.0299711209115851,
          0.0299711209115851,
          0.0299711209115851,
          0.0020573653223050,
          0.0020573653223050,
          0.0020573653223050,
          0.0020573653223050,
          0.0076509102321939,
          0.0076509102321939,
          0.0076509102321939,
          0.0076509102321939,
          0.0241979479733988,
          0.0241979479733988,
          0.0241979479733988,
          0.0241979479733988,
          0.0153790853703291,
          0.0153790853703291,
          0.0153790853703291,
          0.0153790853703291,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0216821248783021,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0139697793510318,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0092822068602285,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0071939002127551,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0028773801106233,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0018821205326004,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628,
          0.0026280760395628 };

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

  return;
}