data = load ( 'cubature_tet_sym_p20_n552_expand_baryc.txt' );
output = fopen ( 'rule20.m', 'w' );
n = size ( data, 1 );

fprintf ( output, 'function [ a, b, c, d, w ] = rule20 ( )\n' );
fprintf ( output, '\n' );
fprintf ( output, '%%*****************************************************************************80\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%%% rule20() returns the tetrahedron rule of precision 20.\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, '%%    02 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, '%%    Jan Jaskowiec, Natarajan Sukumar,\n' );
fprintf ( output, '%%    High order symmetric cubature rules for tetrahedra and pyramids,\n' );
fprintf ( output, '%%    International Journal of Numerical Methods in Engineering,\n' );
fprintf ( output, '%%    Volume 122, Number 1, pages 148-171, 24 August 2020.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%  Output:\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    real a(n), b(n), c(n), d(n): the barycentric coordinates of quadrature points.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '%%    real w(n): the quadrature weights.\n' );
fprintf ( output, '%%\n' );
fprintf ( output, '  a = [ ...\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.16f, ...\n', data(i,1) );
end
fprintf ( output, '    %24.16f ];\n', data(n,1) );

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


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

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

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

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

fclose ( output );
