data = load ( '10-81.txt' );
output = fopen ( 'rule10.m', 'w' );
n = size ( data, 1 );

fprintf ( output, 'function [ a, b, c, d, w ] = rule10 ( )\n' );
fprintf ( output, '\n' );
fprintf ( output, '%%*****************************************************************************80\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%%% rule10() returns the rule of precision 10.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    We suppose we are given a tetrahedron T with vertices A, B, C, D.\n' );
fprintf ( output, '%%    We call a rule with n points, returning barycentric coordinates\n' );
fprintf ( output, '%%    a, b, c, d, and weights w.  Then the integral I of f(x,y,z) over T is \n' );
fprintf ( output, '%%    approximated by Q as follows:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    (x,y,z) = a(1:n) * A + b(1:n) * B + c(1:n) * C + d(1:n * D\n' );
fprintf ( output, '%%    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) \n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Licensing:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    This code is distributed under the MIT license.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Modified:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    24 April 2023\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Author:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    John Burkardt\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Reference:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    Freddie Witherden, Peter Vincent,\n' );
fprintf ( output, '%%    On the identification of symmetric quadrature rules for finite element methods,\n' );
fprintf ( output, '%%    Computers and Mathematics with Applications,\n' );
fprintf ( output, '%%    Volume 69, pages 1232-1241, 2015.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Output:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    a(n), b(n), c(n), d(n): the barycentric coordinates of quadrature points.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    w(n): the quadrature weights, which sum to 1.\n' );
fprintf ( output, '%%\n' );
a = 0.5 * ( data(:,1) + 1.0 );
b = 0.5 * ( data(:,2) + 1.0 );
c = 0.5 * ( data(:,3) + 1.0 );
d = 1.0 - a - b - c;
w = 3.0 * data(:,4) / 4.0;
fprintf ( output, '  a = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f, ...\n', a(i) );
end
fprintf ( output, '    %24.20f ];\n', a(n) );

fprintf ( output, '\n' );
fprintf ( output, '  b = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f, ...\n', b(i) );
end
fprintf ( output, '    %24.20f ];\n', b(n) );


fprintf ( output, '\n' );
fprintf ( output, '  c = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f, ...\n', c(i) );
end
fprintf ( output, '    %24.20f ];\n', c(n) );

fprintf ( output, '\n' );
fprintf ( output, '  d = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f, ...\n', d(i) );
end
fprintf ( output, '    %24.20f ];\n', d(n) );

fprintf ( output, '\n' );
fprintf ( output, '  w = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f, ...\n', w(i) );
end
fprintf ( output, '    %24.20f ];\n', w(n) );

fprintf ( output, '\n' );
fprintf ( output, '  return\n' );
fprintf ( output, 'end\n' );

fclose ( output );
