data = load ( '15-49.txt' );
output = fopen ( 'rule15.m', 'w' );
n = size ( data, 1 );

fprintf ( output, 'function [ a, b, c, w ] = rule15 ( )\n' );
fprintf ( output, '\n' );
fprintf ( output, '%%*****************************************************************************80\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%%% rule15() returns the rule of precision 15.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Discussion:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    The data is given for the following triangle:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%     (0,1) \n' );
fprintf ( output, '%%       |   \\ \n' );
fprintf ( output, '%%       |    \\ \n' );
fprintf ( output, '%%       |     \\ \n' );
fprintf ( output, '%%       |      \\ \n' );
fprintf ( output, '%%     (0,0)----(1,0)\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    We suppose we are given a triangle T with vertices A, B, C.\n' );
fprintf ( output, '%%    We call a rule with n points, returning barycentric coordinates\n' );
fprintf ( output, '%%    a, b, c, and weights w.  Then the integral I of f(x,y) over T is \n' );
fprintf ( output, '%%    approximated by Q as follows:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C\n' );
fprintf ( output, '%%    Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(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, '%%    23 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 a(n), b(n), c(n): the barycentric coordinates of quadrature points.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    real w(n): the quadrature weights, which sum to 1.\n' );
fprintf ( output, '%%\n' );
a = data(:,1);
a = ( a + 1.0 ) / 2.0;
fprintf ( output, '  a = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %20.16f, ...\n', a(i) );
end
fprintf ( output, '    %20.16f ];\n', a(n) );

b = data(:,2);
b = ( b + 1.0 ) / 2.0;
fprintf ( output, '\n' );
fprintf ( output, '  b = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %20.16f, ...\n', b(i) );
end
fprintf ( output, '    %20.16f ];\n', b(n) );

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

w = data(:,3);
w = w / 2.0;
fprintf ( output, '\n' );
fprintf ( output, '  w = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %20.16f, ...\n', w(i) );
end
fprintf ( output, '    %20.16f ];\n', w(n) );

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

fclose ( output );
