subroutine comp_next ( n, k, a, more, h, t )

!*****************************************************************************80
!
!! 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.
!
!    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.
!    This 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, integer N, the integer whose compositions are desired.
!
!    Input, integer K, the number of parts in the composition.
!
!    Input/output, integer A(K), the parts of the composition.
!
!    Input/output, logical MORE, set by the user to start the computation,
!    and by the routine to terminate it.
!
!    Input/output, integer 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!
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer k

  integer a(k)
  integer h
  logical more
  integer n
  integer t
!
!  The first computation.
!
  if ( .not. more ) then

    t = n
    h = 0
    a(1) = n
    a(2:k) = 0
!
!  The next computation.
!
  else
!
!  If the first entry A(1) is positive, then set H to zero,
!  so that when we increment H, it points to A(1); we will decrement A(1) by 1
!  and increment A(2).
!
    if ( 1 < t ) then
      h = 0
    end if
!
!  Otherwise, A(1) is 0.  Then by H + 1 is the entry we incremented last time.
!  Set H = H + 1, zero A(H), adding all but one of its value to A(1),
!  and incrementing A(H+1) by 1.
!
    h = h + 1
    t = a(h)
    a(h) = 0
    a(1) = t - 1
    a(h+1) = a(h+1) + 1

  end if
!
!  This is the last element of the sequence if all the
!  items are in the last slot.
!
  more = ( a(k) /= n )

  return
end
subroutine monomial_value ( d, n, e, x, v )

!*****************************************************************************80
!
!! monomial_value() evaluates a monomial.
!
!  Discussion:
!
!    This routine evaluates a monomial of the form
!
!      product ( 1 <= i <= d ) x(i)^e(i)
!
!    The combination 0.0^0, if encountered, is treated as 1.0.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    16 May 2023
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer D, the spatial dimension.
!
!    integer N, the number of evaluation points.
!
!    integer E(D), the exponents.
!
!    real ( kind = rk ) X(N,D), the point coordinates.
!
!  Output:
!
!    real ( kind = rk ) V(N), the monomial values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer d
  integer n

  integer e(d)
  integer j
  real ( kind = rk ) v(n)
  real ( kind = rk ) x(n,d)

  v(1:n) = 1.0D+00

  do j = 1, d
    if ( 0 /= e(j) ) then
      v(1:n) = v(1:n) * x(1:n,j) ** e(j)
    end if
  end do

  return
end
subroutine prism_unit_monomial_integral ( expon, value )

!*****************************************************************************80
!
!! prism_unit_monomial_integral(): monomial integral in a unit prism.
!
!  Discussion:
!
!    This routine returns the integral of
!
!      product ( 1 <= I <= 3 ) X(I)^EXPON(I)
!
!    over the unit triangular prism.
!
!    A unit triangular prism has a triangular base in the (x,y) plane,
!    projected vertically one unit in the z direction.  
!
!    The vertices are 
!      (1,0,0), (0,1,0), (0,0,0), 
!      (1,0,1), (0,1,1), (0,0,1), 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 May 2023
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer EXPON(3), the exponents.
!
!  Output:
!
!    real ( kind = rk ) VALUE, the integral of the monomial.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer expon(3)
  real ( kind = rk ) r8_factorial
  real ( kind = rk ) value

  value = r8_factorial ( expon(1) ) &
        * r8_factorial ( expon(2) ) &
        / ( expon(3) + 1 ) &
        / r8_factorial ( expon(1) + expon(2) + 2 )

  return
end
function prism_unit_volume ( )

!*****************************************************************************80
!
!! prism_unit_volume() returns the volume of a unit prism.
!
!  Discussion:
!
!    A unit triangular prism has a triangular base in the (x,y) plane,
!    projected vertically one unit in the z direction.  
!
!    The vertices are 
!      (1,0,0), (0,1,0), (0,0,0), 
!      (1,0,1), (0,1,1), (0,0,1), 
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    28 April 2023
!
!  Author:
!
!    John Burkardt
!
!  Output:
!
!    real ( kind = rk ) value: the volume.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) prism_unit_volume

  prism_unit_volume = 0.5D+00

  return
end
subroutine prism_witherden_rule ( p, n, x, y, z, w )

!*****************************************************************************80
!
!! prism_witherden_rule() returns a prism quadrature rule of given precision.
!
!  Discussion:
!
!    A unit triangular prism has a triangular base in the (x,y) plane,
!    projected vertically one unit in the z direction.  
!
!    The vertices are 
!      (1,0,0), (0,1,0), (0,0,0), 
!      (1,0,1), (0,1,1), (0,0,1), 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 May 2023
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Freddie Witherden, Peter Vincent,
!    On the identification of symmetric quadrature rules for finite element methods,
!    Computers and Mathematics with Applications,
!    Volume 69, pages 1232-1241, 2015.
!
!  Input:
!
!    int p: the precision, 0 <= p <= 10.
!
!    int n: the order of the rule.
!
!  Output:
!
!    double x(n), y(n), z(n): the coordinates of quadrature points.
!
!    double w(n): the quadrature weights.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  integer p
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)
  real ( kind = rk ) y(n)
  real ( kind = rk ) z(n)

  if ( p < 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'prism_witherden_rule(): Fatal error!'
    write ( *, '(a)' ) '  Input p < 0.'
    stop ( 1 )
  end if

  if ( 10 < p ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'prism_witherden_rule(): Fatal error!'
    write ( *, '(a)' ) '  Input 10 < p.'
    stop ( 1 )
  end if

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

  return
end
function r8_factorial ( n )

!*****************************************************************************80
!
!! r8_factorial() computes the factorial of N.
!
!  Discussion:
!
!    factorial ( N ) = product ( 1 <= I <= N ) I
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the argument of the factorial function.
!    If N is less than 1, the function value is returned as 1.
!
!  Output:
!
!    real ( kind = rk ) R8_FACTORIAL: the factorial of N.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) r8_factorial
  integer i
  integer n
  real ( kind = rk ) value

  value = 1.0D+00

  do i = 1, n
    value = value * real ( i, kind = rk )
  end do

  r8_factorial = value

  return
end
subroutine rule_order ( p, order )

!*****************************************************************************80
!
!! 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:
!
!    27 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:
!
!    integer p: the precision, 0 <= p <= 10.
!
!  Output:
!
!    integer order: the order of the rule.
!
  implicit none

  integer order
  integer, dimension ( 0 : 10 ) :: order_vec = (/ &
      1, &
      1,   5,   8,  11,  16, &
     28,  35,  46,  60,  85 /)
  integer p

  if ( p < 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'rule_order(): Fatal error.'
    write ( *, '(a)' ) '  Input p < 0.'
  end if

  if ( 10 < p ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'rule_order(): Fatal error.'
    write ( *, '(a)' ) '  Input 10 < p.'
  end if

  order = order_vec(p)

  return
end
subroutine rule00 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule00() returns the prism rule of precision 0.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 1

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.5000000000000000D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          1.0000000000000000D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule01 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule01() returns the prism rule of precision 1.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 1

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.5000000000000000D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          1.0000000000000000D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule02 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule02() returns the prism rule of precision 2.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 5

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.1292091881014018D+00, &
          0.7415816237971964D+00, &
          0.1292091881014018D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.7415816237971964D+00, &
          0.1292091881014018D+00, &
          0.1292091881014018D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.0000000000000001D+00, &
          1.0000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.1666666666666667D+00, &
          0.1666666666666667D+00, &
          0.2222222222222222D+00, &
          0.2222222222222222D+00, &
          0.2222222222222222D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule03 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule03() returns the prism rule of precision 3.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 8

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.2388834359565783D+00, &
          0.7611141478748917D+00, &
          0.0000024161685299D+00, &
          0.7611141478748917D+00, &
          0.0000024161685299D+00, &
          0.2388834359565783D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.7611141478748917D+00, &
          0.2388834359565783D+00, &
          0.7611141478748917D+00, &
          0.0000024161685299D+00, &
          0.2388834359565783D+00, &
          0.0000024161685299D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.0696653977676110D+00, &
          0.9303346022323890D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.2249967381449005D+00, &
          0.2249967381449005D+00, &
          0.0916677539516998D+00, &
          0.0916677539516998D+00, &
          0.0916677539516998D+00, &
          0.0916677539516998D+00, &
          0.0916677539516998D+00, &
          0.0916677539516998D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule04 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule04() returns the prism rule of precision 4.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 11

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.4686558098619952D+00, &
          0.0626883802760096D+00, &
          0.4686558098619952D+00, &
          0.1007404057989106D+00, &
          0.7985191884021787D+00, &
          0.1007404057989106D+00, &
          0.1007404057989106D+00, &
          0.7985191884021787D+00, &
          0.1007404057989106D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.0626883802760096D+00, &
          0.4686558098619952D+00, &
          0.4686558098619952D+00, &
          0.7985191884021787D+00, &
          0.1007404057989106D+00, &
          0.1007404057989106D+00, &
          0.7985191884021787D+00, &
          0.1007404057989106D+00, &
          0.1007404057989106D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.0665690129954826D+00, &
          0.9334309870045174D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.1621800881588701D+00, &
          0.1621800881588701D+00, &
          0.1621800881588701D+00, &
          0.8378199118411298D+00, &
          0.8378199118411298D+00, &
          0.8378199118411298D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.1079119748155355D+00, &
          0.1079119748155355D+00, &
          0.1364146126054776D+00, &
          0.1364146126054776D+00, &
          0.1364146126054776D+00, &
          0.0624887020920827D+00, &
          0.0624887020920827D+00, &
          0.0624887020920827D+00, &
          0.0624887020920827D+00, &
          0.0624887020920827D+00, &
          0.0624887020920827D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule05 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule05() returns the prism rule of precision 5.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 16

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.4872999645502457D+00, &
          0.0254000708995087D+00, &
          0.4872999645502457D+00, &
          0.4455981046703720D+00, &
          0.1088037906592560D+00, &
          0.4455981046703720D+00, &
          0.4455981046703720D+00, &
          0.1088037906592560D+00, &
          0.4455981046703720D+00, &
          0.1008589459827085D+00, &
          0.7982821080345830D+00, &
          0.1008589459827085D+00, &
          0.1008589459827085D+00, &
          0.7982821080345830D+00, &
          0.1008589459827085D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.0254000708995087D+00, &
          0.4872999645502457D+00, &
          0.4872999645502457D+00, &
          0.1088037906592560D+00, &
          0.4455981046703720D+00, &
          0.4455981046703720D+00, &
          0.1088037906592560D+00, &
          0.4455981046703720D+00, &
          0.4455981046703720D+00, &
          0.7982821080345830D+00, &
          0.1008589459827085D+00, &
          0.1008589459827085D+00, &
          0.7982821080345830D+00, &
          0.1008589459827085D+00, &
          0.1008589459827085D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.0644985325672780D+00, &
          0.0644985325672780D+00, &
          0.0644985325672780D+00, &
          0.9355014674327220D+00, &
          0.9355014674327220D+00, &
          0.9355014674327220D+00, &
          0.2147865096474204D+00, &
          0.2147865096474204D+00, &
          0.2147865096474204D+00, &
          0.7852134903525796D+00, &
          0.7852134903525796D+00, &
          0.7852134903525796D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.1778638889828717D+00, &
          0.0561775168070667D+00, &
          0.0561775168070667D+00, &
          0.0561775168070667D+00, &
          0.0464153553290395D+00, &
          0.0464153553290395D+00, &
          0.0464153553290395D+00, &
          0.0464153553290395D+00, &
          0.0464153553290395D+00, &
          0.0464153553290395D+00, &
          0.0625185714369485D+00, &
          0.0625185714369485D+00, &
          0.0625185714369485D+00, &
          0.0625185714369485D+00, &
          0.0625185714369485D+00, &
          0.0625185714369485D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule06 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule06() returns the prism rule of precision 6.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 28

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.1652032650280925D+00, &
          0.6695934699438151D+00, &
          0.1652032650280925D+00, &
          0.0204732696122569D+00, &
          0.9590534607754861D+00, &
          0.0204732696122569D+00, &
          0.4661116655901595D+00, &
          0.0677766688196810D+00, &
          0.4661116655901595D+00, &
          0.4661116655901595D+00, &
          0.0677766688196810D+00, &
          0.4661116655901595D+00, &
          0.2025738967252308D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.2025738967252308D+00, &
          0.2025738967252308D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.2025738967252308D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.6695934699438151D+00, &
          0.1652032650280925D+00, &
          0.1652032650280925D+00, &
          0.9590534607754861D+00, &
          0.0204732696122569D+00, &
          0.0204732696122569D+00, &
          0.0677766688196810D+00, &
          0.4661116655901595D+00, &
          0.4661116655901595D+00, &
          0.0677766688196810D+00, &
          0.4661116655901595D+00, &
          0.4661116655901595D+00, &
          0.7627467408337223D+00, &
          0.2025738967252308D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.2025738967252308D+00, &
          0.0346793624410470D+00, &
          0.7627467408337223D+00, &
          0.2025738967252308D+00, &
          0.7627467408337223D+00, &
          0.0346793624410470D+00, &
          0.2025738967252308D+00, &
          0.0346793624410470D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.2366886475975126D+00, &
          0.7633113524024874D+00, &
          0.0045818495903776D+00, &
          0.9954181504096224D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.2591629325578062D+00, &
          0.2591629325578062D+00, &
          0.2591629325578062D+00, &
          0.7408370674421938D+00, &
          0.7408370674421938D+00, &
          0.7408370674421938D+00, &
          0.0953571033720836D+00, &
          0.0953571033720836D+00, &
          0.0953571033720836D+00, &
          0.0953571033720836D+00, &
          0.0953571033720836D+00, &
          0.0953571033720836D+00, &
          0.9046428966279163D+00, &
          0.9046428966279163D+00, &
          0.9046428966279163D+00, &
          0.9046428966279163D+00, &
          0.9046428966279163D+00, &
          0.9046428966279163D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.0556732835230563D+00, &
          0.0556732835230563D+00, &
          0.0345990273435311D+00, &
          0.0345990273435311D+00, &
          0.0779570037033338D+00, &
          0.0779570037033338D+00, &
          0.0779570037033338D+00, &
          0.0125366338096741D+00, &
          0.0125366338096741D+00, &
          0.0125366338096741D+00, &
          0.0490126400515205D+00, &
          0.0490126400515205D+00, &
          0.0490126400515205D+00, &
          0.0490126400515205D+00, &
          0.0490126400515205D+00, &
          0.0490126400515205D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00, &
          0.0211582187848898D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule07 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule07() returns the prism rule of precision 7.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 35

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.0083375954660215D+00, &
          0.9833248090679571D+00, &
          0.0083375954660215D+00, &
          0.4815219753291366D+00, &
          0.0369560493417268D+00, &
          0.4815219753291366D+00, &
          0.4815219753291366D+00, &
          0.0369560493417268D+00, &
          0.4815219753291366D+00, &
          0.0954832483714894D+00, &
          0.8090335032570213D+00, &
          0.0954832483714894D+00, &
          0.0954832483714894D+00, &
          0.8090335032570213D+00, &
          0.0954832483714894D+00, &
          0.0121491315983783D+00, &
          0.7429966820728956D+00, &
          0.2448541863287261D+00, &
          0.7429966820728956D+00, &
          0.2448541863287261D+00, &
          0.0121491315983783D+00, &
          0.3051562164322261D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.3051562164322261D+00, &
          0.3051562164322261D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.3051562164322261D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.9833248090679571D+00, &
          0.0083375954660215D+00, &
          0.0083375954660215D+00, &
          0.0369560493417268D+00, &
          0.4815219753291366D+00, &
          0.4815219753291366D+00, &
          0.0369560493417268D+00, &
          0.4815219753291366D+00, &
          0.4815219753291366D+00, &
          0.8090335032570213D+00, &
          0.0954832483714894D+00, &
          0.0954832483714894D+00, &
          0.8090335032570213D+00, &
          0.0954832483714894D+00, &
          0.0954832483714894D+00, &
          0.7429966820728956D+00, &
          0.0121491315983783D+00, &
          0.7429966820728956D+00, &
          0.2448541863287261D+00, &
          0.0121491315983783D+00, &
          0.2448541863287261D+00, &
          0.1529845984247976D+00, &
          0.3051562164322261D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.3051562164322261D+00, &
          0.5418591851429763D+00, &
          0.1529845984247976D+00, &
          0.3051562164322261D+00, &
          0.1529845984247976D+00, &
          0.5418591851429763D+00, &
          0.3051562164322261D+00, &
          0.5418591851429763D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.0098859652045542D+00, &
          0.9901140347954458D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.0793063228967370D+00, &
          0.0793063228967370D+00, &
          0.0793063228967370D+00, &
          0.9206936771032630D+00, &
          0.9206936771032630D+00, &
          0.9206936771032630D+00, &
          0.1020754547065085D+00, &
          0.1020754547065085D+00, &
          0.1020754547065085D+00, &
          0.8979245452934915D+00, &
          0.8979245452934915D+00, &
          0.8979245452934915D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.2980327197339451D+00, &
          0.2980327197339451D+00, &
          0.2980327197339451D+00, &
          0.2980327197339451D+00, &
          0.2980327197339451D+00, &
          0.2980327197339451D+00, &
          0.7019672802660549D+00, &
          0.7019672802660549D+00, &
          0.7019672802660549D+00, &
          0.7019672802660549D+00, &
          0.7019672802660549D+00, &
          0.7019672802660549D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.0284456811126893D+00, &
          0.0284456811126893D+00, &
          0.0061182421730950D+00, &
          0.0061182421730950D+00, &
          0.0061182421730950D+00, &
          0.0215508640862265D+00, &
          0.0215508640862265D+00, &
          0.0215508640862265D+00, &
          0.0215508640862265D+00, &
          0.0215508640862265D+00, &
          0.0215508640862265D+00, &
          0.0291785249020985D+00, &
          0.0291785249020985D+00, &
          0.0291785249020985D+00, &
          0.0291785249020985D+00, &
          0.0291785249020985D+00, &
          0.0291785249020985D+00, &
          0.0255148563351493D+00, &
          0.0255148563351493D+00, &
          0.0255148563351493D+00, &
          0.0255148563351493D+00, &
          0.0255148563351493D+00, &
          0.0255148563351493D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00, &
          0.0389407032762076D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule08 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule08() returns the prism rule of precision 8.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 46

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.0534123509369072D+00, &
          0.8931752981261858D+00, &
          0.0534123509369072D+00, &
          0.4600889628137106D+00, &
          0.0798220743725788D+00, &
          0.4600889628137106D+00, &
          0.4585690687909513D+00, &
          0.0828618624180973D+00, &
          0.4585690687909513D+00, &
          0.4585690687909513D+00, &
          0.0828618624180973D+00, &
          0.4585690687909513D+00, &
          0.1740616079243704D+00, &
          0.6518767841512592D+00, &
          0.1740616079243704D+00, &
          0.1740616079243704D+00, &
          0.6518767841512592D+00, &
          0.1740616079243704D+00, &
          0.0472387858397694D+00, &
          0.9055224283204613D+00, &
          0.0472387858397694D+00, &
          0.0472387858397694D+00, &
          0.9055224283204613D+00, &
          0.0472387858397694D+00, &
          0.1597492639425890D+00, &
          0.6805014721148220D+00, &
          0.1597492639425890D+00, &
          0.1597492639425890D+00, &
          0.6805014721148220D+00, &
          0.1597492639425890D+00, &
          0.2628138006912410D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.2628138006912410D+00, &
          0.2628138006912410D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.2628138006912410D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.8931752981261858D+00, &
          0.0534123509369072D+00, &
          0.0534123509369072D+00, &
          0.0798220743725788D+00, &
          0.4600889628137106D+00, &
          0.4600889628137106D+00, &
          0.0828618624180973D+00, &
          0.4585690687909513D+00, &
          0.4585690687909513D+00, &
          0.0828618624180973D+00, &
          0.4585690687909513D+00, &
          0.4585690687909513D+00, &
          0.6518767841512592D+00, &
          0.1740616079243704D+00, &
          0.1740616079243704D+00, &
          0.6518767841512592D+00, &
          0.1740616079243704D+00, &
          0.1740616079243704D+00, &
          0.9055224283204613D+00, &
          0.0472387858397694D+00, &
          0.0472387858397694D+00, &
          0.9055224283204613D+00, &
          0.0472387858397694D+00, &
          0.0472387858397694D+00, &
          0.6805014721148220D+00, &
          0.1597492639425890D+00, &
          0.1597492639425890D+00, &
          0.6805014721148220D+00, &
          0.1597492639425890D+00, &
          0.1597492639425890D+00, &
          0.7285980718010000D+00, &
          0.2628138006912410D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.2628138006912410D+00, &
          0.0085881275077590D+00, &
          0.7285980718010000D+00, &
          0.2628138006912410D+00, &
          0.7285980718010000D+00, &
          0.0085881275077590D+00, &
          0.2628138006912410D+00, &
          0.0085881275077590D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.3877830669470047D+00, &
          0.6122169330529953D+00, &
          0.1590872792145671D+00, &
          0.8409127207854329D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.0619020449998729D+00, &
          0.0619020449998729D+00, &
          0.0619020449998729D+00, &
          0.9380979550001272D+00, &
          0.9380979550001272D+00, &
          0.9380979550001272D+00, &
          0.2969578853227048D+00, &
          0.2969578853227048D+00, &
          0.2969578853227048D+00, &
          0.7030421146772952D+00, &
          0.7030421146772952D+00, &
          0.7030421146772952D+00, &
          0.0855922784706599D+00, &
          0.0855922784706599D+00, &
          0.0855922784706599D+00, &
          0.9144077215293400D+00, &
          0.9144077215293400D+00, &
          0.9144077215293400D+00, &
          0.0170674167296728D+00, &
          0.0170674167296728D+00, &
          0.0170674167296728D+00, &
          0.9829325832703272D+00, &
          0.9829325832703272D+00, &
          0.9829325832703272D+00, &
          0.2113248654051871D+00, &
          0.2113248654051871D+00, &
          0.2113248654051871D+00, &
          0.2113248654051871D+00, &
          0.2113248654051871D+00, &
          0.2113248654051871D+00, &
          0.7886751345948129D+00, &
          0.7886751345948129D+00, &
          0.7886751345948129D+00, &
          0.7886751345948129D+00, &
          0.7886751345948129D+00, &
          0.7886751345948129D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.0207349366428555D+00, &
          0.0207349366428555D+00, &
          0.0510894271121815D+00, &
          0.0510894271121815D+00, &
          0.0181522487841497D+00, &
          0.0181522487841497D+00, &
          0.0181522487841497D+00, &
          0.0526997449065072D+00, &
          0.0526997449065072D+00, &
          0.0526997449065072D+00, &
          0.0210865092935865D+00, &
          0.0210865092935865D+00, &
          0.0210865092935865D+00, &
          0.0210865092935865D+00, &
          0.0210865092935865D+00, &
          0.0210865092935865D+00, &
          0.0406292626589893D+00, &
          0.0406292626589893D+00, &
          0.0406292626589893D+00, &
          0.0406292626589893D+00, &
          0.0406292626589893D+00, &
          0.0406292626589893D+00, &
          0.0071453577535386D+00, &
          0.0071453577535386D+00, &
          0.0071453577535386D+00, &
          0.0071453577535386D+00, &
          0.0071453577535386D+00, &
          0.0071453577535386D+00, &
          0.0111430273484653D+00, &
          0.0111430273484653D+00, &
          0.0111430273484653D+00, &
          0.0111430273484653D+00, &
          0.0111430273484653D+00, &
          0.0111430273484653D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00, &
          0.0136475290908731D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule09 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule09() returns the prism rule of precision 9.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 60

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.4994650351706859D+00, &
          0.0010699296586282D+00, &
          0.4994650351706859D+00, &
          0.4450675000425037D+00, &
          0.1098649999149927D+00, &
          0.4450675000425037D+00, &
          0.4450675000425037D+00, &
          0.1098649999149927D+00, &
          0.4450675000425037D+00, &
          0.2111544337105105D+00, &
          0.5776911325789790D+00, &
          0.2111544337105105D+00, &
          0.2111544337105105D+00, &
          0.5776911325789790D+00, &
          0.2111544337105105D+00, &
          0.0438424814414928D+00, &
          0.9123150371170143D+00, &
          0.0438424814414928D+00, &
          0.0438424814414928D+00, &
          0.9123150371170143D+00, &
          0.0438424814414928D+00, &
          0.4827499316708669D+00, &
          0.0345001366582662D+00, &
          0.4827499316708669D+00, &
          0.4827499316708669D+00, &
          0.0345001366582662D+00, &
          0.4827499316708669D+00, &
          0.1771782441116773D+00, &
          0.6456435117766455D+00, &
          0.1771782441116773D+00, &
          0.1771782441116773D+00, &
          0.6456435117766455D+00, &
          0.1771782441116773D+00, &
          0.0462311804109438D+00, &
          0.9075376391781124D+00, &
          0.0462311804109438D+00, &
          0.0462311804109438D+00, &
          0.9075376391781124D+00, &
          0.0462311804109438D+00, &
          0.0593670097433656D+00, &
          0.7246952695529995D+00, &
          0.2159377207036349D+00, &
          0.7246952695529995D+00, &
          0.2159377207036349D+00, &
          0.0593670097433656D+00, &
          0.2241275457533561D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.2241275457533561D+00, &
          0.2241275457533561D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.2241275457533561D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.0010699296586282D+00, &
          0.4994650351706859D+00, &
          0.4994650351706859D+00, &
          0.1098649999149927D+00, &
          0.4450675000425037D+00, &
          0.4450675000425037D+00, &
          0.1098649999149927D+00, &
          0.4450675000425037D+00, &
          0.4450675000425037D+00, &
          0.5776911325789790D+00, &
          0.2111544337105105D+00, &
          0.2111544337105105D+00, &
          0.5776911325789790D+00, &
          0.2111544337105105D+00, &
          0.2111544337105105D+00, &
          0.9123150371170143D+00, &
          0.0438424814414928D+00, &
          0.0438424814414928D+00, &
          0.9123150371170143D+00, &
          0.0438424814414928D+00, &
          0.0438424814414928D+00, &
          0.0345001366582662D+00, &
          0.4827499316708669D+00, &
          0.4827499316708669D+00, &
          0.0345001366582662D+00, &
          0.4827499316708669D+00, &
          0.4827499316708669D+00, &
          0.6456435117766455D+00, &
          0.1771782441116773D+00, &
          0.1771782441116773D+00, &
          0.6456435117766455D+00, &
          0.1771782441116773D+00, &
          0.1771782441116773D+00, &
          0.9075376391781124D+00, &
          0.0462311804109438D+00, &
          0.0462311804109438D+00, &
          0.9075376391781124D+00, &
          0.0462311804109438D+00, &
          0.0462311804109438D+00, &
          0.7246952695529995D+00, &
          0.0593670097433656D+00, &
          0.7246952695529995D+00, &
          0.2159377207036349D+00, &
          0.0593670097433656D+00, &
          0.2159377207036349D+00, &
          0.7484054595660388D+00, &
          0.2241275457533561D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.2241275457533561D+00, &
          0.0274669946806053D+00, &
          0.7484054595660388D+00, &
          0.2241275457533561D+00, &
          0.7484054595660388D+00, &
          0.0274669946806053D+00, &
          0.2241275457533561D+00, &
          0.0274669946806053D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.5000000000000000D+00, &
          0.0620804206000612D+00, &
          0.9379195793999389D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.2351832223969616D+00, &
          0.2351832223969616D+00, &
          0.2351832223969616D+00, &
          0.7648167776030383D+00, &
          0.7648167776030383D+00, &
          0.7648167776030383D+00, &
          0.3380115589156150D+00, &
          0.3380115589156150D+00, &
          0.3380115589156150D+00, &
          0.6619884410843850D+00, &
          0.6619884410843850D+00, &
          0.6619884410843850D+00, &
          0.2991741657146056D+00, &
          0.2991741657146056D+00, &
          0.2991741657146056D+00, &
          0.7008258342853944D+00, &
          0.7008258342853944D+00, &
          0.7008258342853944D+00, &
          0.0161039372678948D+00, &
          0.0161039372678948D+00, &
          0.0161039372678948D+00, &
          0.9838960627321052D+00, &
          0.9838960627321052D+00, &
          0.9838960627321052D+00, &
          0.0546091297715147D+00, &
          0.0546091297715147D+00, &
          0.0546091297715147D+00, &
          0.9453908702284853D+00, &
          0.9453908702284853D+00, &
          0.9453908702284853D+00, &
          0.0211468980535365D+00, &
          0.0211468980535365D+00, &
          0.0211468980535365D+00, &
          0.9788531019464635D+00, &
          0.9788531019464635D+00, &
          0.9788531019464635D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.1523897230641826D+00, &
          0.1523897230641826D+00, &
          0.1523897230641826D+00, &
          0.1523897230641826D+00, &
          0.1523897230641826D+00, &
          0.1523897230641826D+00, &
          0.8476102769358174D+00, &
          0.8476102769358174D+00, &
          0.8476102769358174D+00, &
          0.8476102769358174D+00, &
          0.8476102769358174D+00, &
          0.8476102769358174D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.0477869971115503D+00, &
          0.0265242298973318D+00, &
          0.0265242298973318D+00, &
          0.0123656013604738D+00, &
          0.0123656013604738D+00, &
          0.0123656013604738D+00, &
          0.0377019831531896D+00, &
          0.0377019831531896D+00, &
          0.0377019831531896D+00, &
          0.0377019831531896D+00, &
          0.0377019831531896D+00, &
          0.0377019831531896D+00, &
          0.0233456174363846D+00, &
          0.0233456174363846D+00, &
          0.0233456174363846D+00, &
          0.0233456174363846D+00, &
          0.0233456174363846D+00, &
          0.0233456174363846D+00, &
          0.0095970715963016D+00, &
          0.0095970715963016D+00, &
          0.0095970715963016D+00, &
          0.0095970715963016D+00, &
          0.0095970715963016D+00, &
          0.0095970715963016D+00, &
          0.0066321078005801D+00, &
          0.0066321078005801D+00, &
          0.0066321078005801D+00, &
          0.0066321078005801D+00, &
          0.0066321078005801D+00, &
          0.0066321078005801D+00, &
          0.0166799856433621D+00, &
          0.0166799856433621D+00, &
          0.0166799856433621D+00, &
          0.0166799856433621D+00, &
          0.0166799856433621D+00, &
          0.0166799856433621D+00, &
          0.0030327934591018D+00, &
          0.0030327934591018D+00, &
          0.0030327934591018D+00, &
          0.0030327934591018D+00, &
          0.0030327934591018D+00, &
          0.0030327934591018D+00, &
          0.0219227903282098D+00, &
          0.0219227903282098D+00, &
          0.0219227903282098D+00, &
          0.0219227903282098D+00, &
          0.0219227903282098D+00, &
          0.0219227903282098D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00, &
          0.0123828035424656D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end
subroutine rule10 ( n, x, y, z, w )

!*****************************************************************************80
!
!! rule10() returns the prism rule of precision 10.
!
!  Discussion:
!
!    We suppose we are given a triangular prism P with vertices
!    (1,0,0), (0,1,0), (0,0,0),
!    (1,0,1), (0,1,1), (0,0,1),
!    We call a rule with n points, returning coordinates
!    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is 
!    approximated by Q as follows:
!
!    Q = volume(P) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 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.
!
!  Output:
!
!    real x(n), y(n), z(n): the coordinates of quadrature points.
!
!    real w(n): the quadrature weights, which sum to 1.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer, parameter :: n_save = 85

  real ( kind = rk ) x(n)
  real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/  &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.2074503749892171D+00, &
          0.5850992500215657D+00, &
          0.2074503749892171D+00, &
          0.4561273628680818D+00, &
          0.0877452742638365D+00, &
          0.4561273628680818D+00, &
          0.0833192521500528D+00, &
          0.8333614956998945D+00, &
          0.0833192521500528D+00, &
          0.2071663521398311D+00, &
          0.5856672957203378D+00, &
          0.2071663521398311D+00, &
          0.2071663521398311D+00, &
          0.5856672957203378D+00, &
          0.2071663521398311D+00, &
          0.1429911965428408D+00, &
          0.7140176069143184D+00, &
          0.1429911965428408D+00, &
          0.1429911965428408D+00, &
          0.7140176069143184D+00, &
          0.1429911965428408D+00, &
          0.4970408872397459D+00, &
          0.0059182255205082D+00, &
          0.4970408872397459D+00, &
          0.4970408872397459D+00, &
          0.0059182255205082D+00, &
          0.4970408872397459D+00, &
          0.4277824598534921D+00, &
          0.1444350802930158D+00, &
          0.4277824598534921D+00, &
          0.4277824598534921D+00, &
          0.1444350802930158D+00, &
          0.4277824598534921D+00, &
          0.0367246894695538D+00, &
          0.9265506210608924D+00, &
          0.0367246894695538D+00, &
          0.0367246894695538D+00, &
          0.9265506210608924D+00, &
          0.0367246894695538D+00, &
          0.0018926195990494D+00, &
          0.9979100409523263D+00, &
          0.0001973394486243D+00, &
          0.9979100409523263D+00, &
          0.0001973394486243D+00, &
          0.0018926195990494D+00, &
          0.0432896837794919D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.0432896837794919D+00, &
          0.0432896837794919D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.0432896837794919D+00, &
          0.0382915191722036D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.0382915191722036D+00, &
          0.0382915191722036D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.0382915191722036D+00, &
          0.0160393936272334D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.0160393936272334D+00, &
          0.0160393936272334D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.0160393936272334D+00 /)

  real ( kind = rk ) y(n)
  real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.3333333333333334D+00, &
          0.5850992500215657D+00, &
          0.2074503749892171D+00, &
          0.2074503749892171D+00, &
          0.0877452742638365D+00, &
          0.4561273628680818D+00, &
          0.4561273628680818D+00, &
          0.8333614956998945D+00, &
          0.0833192521500528D+00, &
          0.0833192521500528D+00, &
          0.5856672957203378D+00, &
          0.2071663521398311D+00, &
          0.2071663521398311D+00, &
          0.5856672957203378D+00, &
          0.2071663521398311D+00, &
          0.2071663521398311D+00, &
          0.7140176069143184D+00, &
          0.1429911965428408D+00, &
          0.1429911965428408D+00, &
          0.7140176069143184D+00, &
          0.1429911965428408D+00, &
          0.1429911965428408D+00, &
          0.0059182255205082D+00, &
          0.4970408872397459D+00, &
          0.4970408872397459D+00, &
          0.0059182255205082D+00, &
          0.4970408872397459D+00, &
          0.4970408872397459D+00, &
          0.1444350802930158D+00, &
          0.4277824598534921D+00, &
          0.4277824598534921D+00, &
          0.1444350802930158D+00, &
          0.4277824598534921D+00, &
          0.4277824598534921D+00, &
          0.9265506210608924D+00, &
          0.0367246894695538D+00, &
          0.0367246894695538D+00, &
          0.9265506210608924D+00, &
          0.0367246894695538D+00, &
          0.0367246894695538D+00, &
          0.9979100409523263D+00, &
          0.0018926195990494D+00, &
          0.9979100409523263D+00, &
          0.0001973394486243D+00, &
          0.0018926195990494D+00, &
          0.0001973394486243D+00, &
          0.6627532432395955D+00, &
          0.0432896837794919D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.0432896837794919D+00, &
          0.2939570729809126D+00, &
          0.6627532432395955D+00, &
          0.0432896837794919D+00, &
          0.6627532432395955D+00, &
          0.2939570729809126D+00, &
          0.0432896837794919D+00, &
          0.2939570729809126D+00, &
          0.6875657341783573D+00, &
          0.0382915191722036D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.0382915191722036D+00, &
          0.2741427466494392D+00, &
          0.6875657341783573D+00, &
          0.0382915191722036D+00, &
          0.6875657341783573D+00, &
          0.2741427466494392D+00, &
          0.0382915191722036D+00, &
          0.2741427466494392D+00, &
          0.8702636086393689D+00, &
          0.0160393936272334D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.0160393936272334D+00, &
          0.1136969977333976D+00, &
          0.8702636086393689D+00, &
          0.0160393936272334D+00, &
          0.8702636086393689D+00, &
          0.1136969977333976D+00, &
          0.0160393936272334D+00, &
          0.1136969977333976D+00 /)

  real ( kind = rk ) z(n)
  real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ &
          0.0000000000000000D+00, &
          1.0000000000000000D+00, &
          0.3580817993007142D+00, &
          0.6419182006992858D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.2608548695498160D+00, &
          0.2608548695498160D+00, &
          0.2608548695498160D+00, &
          0.7391451304501840D+00, &
          0.7391451304501840D+00, &
          0.7391451304501840D+00, &
          0.0739666680004302D+00, &
          0.0739666680004302D+00, &
          0.0739666680004302D+00, &
          0.9260333319995698D+00, &
          0.9260333319995698D+00, &
          0.9260333319995698D+00, &
          0.1862549441760233D+00, &
          0.1862549441760233D+00, &
          0.1862549441760233D+00, &
          0.8137450558239767D+00, &
          0.8137450558239767D+00, &
          0.8137450558239767D+00, &
          0.1155140587235295D+00, &
          0.1155140587235295D+00, &
          0.1155140587235295D+00, &
          0.8844859412764705D+00, &
          0.8844859412764705D+00, &
          0.8844859412764705D+00, &
          0.0273190720391215D+00, &
          0.0273190720391215D+00, &
          0.0273190720391215D+00, &
          0.9726809279608786D+00, &
          0.9726809279608786D+00, &
          0.9726809279608786D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.5000000000000000D+00, &
          0.0288182662745812D+00, &
          0.0288182662745812D+00, &
          0.0288182662745812D+00, &
          0.0288182662745812D+00, &
          0.0288182662745812D+00, &
          0.0288182662745812D+00, &
          0.9711817337254187D+00, &
          0.9711817337254187D+00, &
          0.9711817337254187D+00, &
          0.9711817337254187D+00, &
          0.9711817337254187D+00, &
          0.9711817337254187D+00, &
          0.3323319456663029D+00, &
          0.3323319456663029D+00, &
          0.3323319456663029D+00, &
          0.3323319456663029D+00, &
          0.3323319456663029D+00, &
          0.3323319456663029D+00, &
          0.6676680543336971D+00, &
          0.6676680543336971D+00, &
          0.6676680543336971D+00, &
          0.6676680543336971D+00, &
          0.6676680543336971D+00, &
          0.6676680543336971D+00, &
          0.1928382490614613D+00, &
          0.1928382490614613D+00, &
          0.1928382490614613D+00, &
          0.1928382490614613D+00, &
          0.1928382490614613D+00, &
          0.1928382490614613D+00, &
          0.8071617509385387D+00, &
          0.8071617509385387D+00, &
          0.8071617509385387D+00, &
          0.8071617509385387D+00, &
          0.8071617509385387D+00, &
          0.8071617509385387D+00 /)

  real ( kind = rk ) w(n)
  real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ &
          0.0095107770210602D+00, &
          0.0095107770210602D+00, &
          0.0320123016081395D+00, &
          0.0320123016081395D+00, &
          0.0113557896200342D+00, &
          0.0113557896200342D+00, &
          0.0113557896200342D+00, &
          0.0301966545250866D+00, &
          0.0301966545250866D+00, &
          0.0301966545250866D+00, &
          0.0176309095004028D+00, &
          0.0176309095004028D+00, &
          0.0176309095004028D+00, &
          0.0260882953260851D+00, &
          0.0260882953260851D+00, &
          0.0260882953260851D+00, &
          0.0260882953260851D+00, &
          0.0260882953260851D+00, &
          0.0260882953260851D+00, &
          0.0113382624425785D+00, &
          0.0113382624425785D+00, &
          0.0113382624425785D+00, &
          0.0113382624425785D+00, &
          0.0113382624425785D+00, &
          0.0113382624425785D+00, &
          0.0070316434137385D+00, &
          0.0070316434137385D+00, &
          0.0070316434137385D+00, &
          0.0070316434137385D+00, &
          0.0070316434137385D+00, &
          0.0070316434137385D+00, &
          0.0243250168786856D+00, &
          0.0243250168786856D+00, &
          0.0243250168786856D+00, &
          0.0243250168786856D+00, &
          0.0243250168786856D+00, &
          0.0243250168786856D+00, &
          0.0018924357100839D+00, &
          0.0018924357100839D+00, &
          0.0018924357100839D+00, &
          0.0018924357100839D+00, &
          0.0018924357100839D+00, &
          0.0018924357100839D+00, &
          0.0007802159100457D+00, &
          0.0007802159100457D+00, &
          0.0007802159100457D+00, &
          0.0007802159100457D+00, &
          0.0007802159100457D+00, &
          0.0007802159100457D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0060458851624517D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0136507771841729D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00, &
          0.0061923846298526D+00 /)

  x = x_save
  y = y_save
  z = z_save
  w = w_save

  return
end