data = load ( 'cubature_tet_sym_p3_n8_expand_baryc.txt' );
output = fopen ( 'rule03.c', 'w' );
n = size ( data, 1 );

fprintf ( output, 'void rule03 ( int n, double a[], double b[], double c[], double d[], double w[] )\n' ); 
fprintf ( output, '\n' );
fprintf ( output, '/******************************************************************************/\n' );
fprintf ( output, '/*\n' );
fprintf ( output, '  Purpose:\n' );
fprintf ( output, '\n' );
fprintf ( output, '   rule03() returns the rule of precision 3.\n' );
fprintf ( output, '\n' );
fprintf ( output, '  Licensing:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    This code is distributed under the GNU LGPL 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, '  Input:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    int n: the number of quadrature points for this rule.\n' );
fprintf ( output, '\n' );
fprintf ( output, '  Output:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    double a[n], b[n], c[n], d[n]: the barycentric coordinates of quadrature points.\n' );
fprintf ( output, '\n' );
fprintf ( output, '    double w[n]: the quadrature weights.\n' );
fprintf ( output, '*/\n' );
fprintf ( output, '{\n' );
fprintf ( output, '  static double a_save[] = {\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f,\n', data(i,1) );
end
fprintf ( output, '    %24.20f };\n', data(n,1) );

fprintf ( output, '\n' );
fprintf ( output, '  static double b_save[] = {\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f,\n', data(i,2) );
end
fprintf ( output, '    %24.20f };\n', data(n,2) );

fprintf ( output, '\n' );
fprintf ( output, '  static double c_save[] = {\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f,\n', data(i,3) );
end
fprintf ( output, '    %24.20f };\n', data(n,3) );

fprintf ( output, '\n' );
fprintf ( output, '  static double d_save[] = {\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f,\n', data(i,4) );
end
fprintf ( output, '    %24.20f };\n', data(n,4) );

fprintf ( output, '\n' );
fprintf ( output, '  double w_save[] = {\n' );
for i = 1 : n - 1
  fprintf ( output, '    %24.20f,\n', data(i,5) );
end
fprintf ( output, '    %24.20f };\n', data(n,5) );
fprintf ( output, '\n' );

fprintf ( output, '  r8vec_copy ( n, a_save, a );\n' );
fprintf ( output, '  r8vec_copy ( n, b_save, b );\n' );
fprintf ( output, '  r8vec_copy ( n, c_save, c );\n' );
fprintf ( output, '  r8vec_copy ( n, d_save, d );\n' );
fprintf ( output, '  r8vec_copy ( n, w_save, w );\n' );

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

fclose ( output );
