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

fprintf ( output, 'function [ x, y, z, w ] = rule10 ( )\n' );
fprintf ( output, '\n' );
fprintf ( output, '%%*****************************************************************************80\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%%% rule10() returns the prism rule of precision 10.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Discussion:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    We suppose we are given a triangular prism P with vertices\n' );
fprintf ( output, '%%    (1,0,0), (0,1,0), (0,0,0),\n' );
fprintf ( output, '%%    (1,0,1), (0,1,1), (0,0,1),\n' );
fprintf ( output, '%%    We call a rule with n points, returning coordinates\n' );
fprintf ( output, '%%    x, y, z, and weights w.  Then the integral I of f(x,y,z) over P is \n' );
fprintf ( output, '%%    approximated by Q as follows:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    Q = volume(P) * 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, '%%    27 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, '%%    real x(n), y(n), z(n): the coordinates of quadrature points.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    real w(n): the quadrature weights, which sum to 1.\n' );
fprintf ( output, '%%\n' );
x = 0.5 * ( data(:,1) + 1.0 );
y = 0.5 * ( data(:,2) + 1.0 );
z = 0.5 * ( data(:,3) + 1.0 );
w = data(:,4) / 4.0;
fprintf ( output, '  x = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.16f, ...\n', x(i) );
end
fprintf ( output, '    %24.16f ];\n', x(n) );

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

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

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

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

fclose ( output );
