function R = biharmonic_exact_r3 ( X, Y, a, b, c, d, e, f, g )

%*****************************************************************************80
%
%% biharmonic_exact_r3() evaluates exact biharmonic residual for W(X,Y) #3.
%
%  Discussion:
%
%   Note the formula for W:
%
%   R = sqrt ( ( X - e ).^2 + ( Y - f ).^2 );
%
%   W =   a     * R.^2 * log ( R ) ...
%       + b     * R.^2 ...
%       + c * log ( R ) ...
%       + d;
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 August 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, Y: the coordinates of the points.
%
%    real a, b, c, d, e, f, g: parameters.
%
%  Output:
%
%    real R: the residual evaluated at the points (X,Y).
%
  Wxxxx =  2 * ( 8 * a * ( e - X ).^4 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - 12 * a * ( e - X ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) + 3 * a - 24 * c * ( e - X ).^4 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^3 + 24 * c * ( e - X ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - 3 * c ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) ) ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 );

  Wxxyy =  2 * ( 8 * a * ( e - X ).^2 .* ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - 2 * a * ( e - X ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) - 2 * a * ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) + a - 24 * c * ( e - X ).^2 .* ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^3 + 4 * c * ( e - X ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 + 4 * c * ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - c ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) ) ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 );

  Wyyyy =  2 * ( 8 * a * ( f - Y ).^4 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - 12 * a * ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) + 3 * a - 24 * c * ( f - Y ).^4 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^3 + 24 * c * ( f - Y ).^2 ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ).^2 - 3 * c ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 ) ) ...
    ./ ( ( e - X ).^2 + ( f - Y ).^2 );

  R = Wxxxx + 2.0 * Wxxyy + Wyyyy;

  return
end
