function sphere_surface ( )

  t = 117.28612573;
  n = 1;
  for logn = 0 : 6
    q = sphere_surface_quad ( n );
    e = abs ( q - t );
    fprintf ( 1, '  %7d  %14.6g  %14.6g\n', n, q, e );
    n = n * 10;
  end

  return
end
function value = sphere_surface_quad ( n )

  c = [ 1.0; 1.0; 1.0 ];
  r = 2.0;
%
%  Random points on a unit sphere.
%
  x = randn ( 3, n );
  norm2 = x.^2;
  s = sqrt ( sum ( norm2 ) );
  s3 = repmat ( s, 3, 1 );
  x = x ./ s3;
%
%  Shift to a sphere of radius R and center C.
%
  cn = repmat ( c, 1, n );
  x = cn + r * x;

  v = 4 * pi * r^2;
  value = ( v / n ) * sum ( f ( n, x ) );

  return
end
function value = f ( n, x )

  f = zeros ( n, 1 );

  value = ( x(1,:).^2 )';

  return
end
