function tetrahedron_witherden_rule_test04 ( )

%*****************************************************************************80
%
%% tetrahedron_witherden_rule_test04() integrates 1/sqrt(r).
%
%  Discussion:
%
%    The integral is taken over the reference tetrahedron with vertices
%    (0,0,0), (1,0,0), (0,1,0), (0,0,1).
%
%    The test integral is taken from Jaskowiec and Sukumar.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    23 May 2023
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Jan Jaskowiec, Natarajan Sukumar,
%    High order symmetric cubature rules for tetrahedra and pyramids,
%    International Journal of Numerical Methods in Engineering,
%    Volume 122, Number 1, pages 148-171, 24 August 2020.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'tetrahedron_witherden_rule_test04():\n' );
  fprintf ( 1, '  Integrate 1/sqrt(r) over the reference tetrahedron.\n' );
  fprintf ( 1, '  Witherden rule #9 fails because a quadrature point is\n' );
  fprintf ( 1, '  very near the singularity at the origin.\n' );

  exact = 163.0 / 679.0;
  fprintf ( 1, '  Exact integral value is %g.\n', exact );

  tetra = [ 0.0, 0.0, 0.0; ...
            1.0, 0.0, 0.0; ...
            0.0, 1.0, 0.0; ...
            0.0, 0.0, 1.0 ];
  volume = tetrahedron_volume ( tetra );
  fprintf ( 1, '  Volume of tetrahedron is %g.\n', volume );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   P    N     Q     |Q-Exact]\n' );
  fprintf ( 1, '\n' );
  for p = 0 : 10

    [ n, a, b, c, d, w ] = tetrahedron_witherden_rule ( p );
    r = sqrt ( a.^2 + b.^2 + c.^2 );
    f = 1.0 ./ sqrt ( r );
    q = volume * sum ( w .* f );
    e = abs ( q - exact );
    fprintf ( 1, '  %2d  %3d  %20.16g  %20.16g\n', p, n, q, e )

  end

  return
end
