data = load ( '11-90.txt' );
output = fopen ( 'rule11.c', 'w' );
n = size ( data, 1 );

fprintf ( output, 'void rule11 ( int n, double x[], double y[], double z[], double w[] )\n' ); 
fprintf ( output, '\n' );
fprintf ( output, '/******************************************************************************/\n' );
fprintf ( output, '/*\n' );
fprintf ( output, '  Purpose:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    rule11() returns the hexahedron rule of precision 11.\n' );
fprintf ( output, '\n' );
fprintf ( output, '  Discussion:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    We suppose we are given a hexahedron H with vertices\n' );
fprintf ( output, '      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,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\n' );
fprintf ( output, '    H is approximated by Q as follows:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    Q = volume(H) * 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 GNU LGPL license.\n' );
fprintf ( output, '\n' );
fprintf ( output, '  Modified:\n' );
fprintf ( output, '\n' );
fprintf ( output, '    02 May 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, '  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 x(n), y(n), x(n): the coordinates of quadrature points.\n' );
fprintf ( output, '\n' );
fprintf ( output, '    double w[n]: the quadrature weights, which sum to 1.\n' );
fprintf ( output, '*/\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) / 8.0;
fprintf ( output, '  static double x_save[] = {\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, '  static double y_save[] = {\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, '  static double z_save[] = {\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, '  double w_save[] = {\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, '  r8vec_copy ( n, x_save, x );\n' );
fprintf ( output, '  r8vec_copy ( n, y_save, y );\n' );
fprintf ( output, '  r8vec_copy ( n, z_save, z );\n' );
fprintf ( output, '  r8vec_copy ( n, w_save, w );\n' );

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

fclose ( output );
