function rd_test ( )

%*****************************************************************************80
%
%% rd_test() tests rd().
%
%  Discussion:
%
%    This driver tests the real ( kind = 8 ) function subroutine for the
%    integral RD(X,Y,Z), which is symmetric in X and Y.  The first
%    twelve sets of values of X, Y, Z are extreme points of the region of
%    valid arguments defined by the machine-dependent constants LOLIM
%    and UPLIM.  The values of LOLIM, UPLIM, X, Y, Z, and ERRTOL (see
%    comments in subroutine) may be used on most machines but provide a
%    severe test of robustness only on the ibm 360/370 series.  The
%    thirteenth set tests the failure exit.  The fourteenth set is a
%    check value: RD(0,2,1) = 3B = 3(PI)/4A, where A and B are the
%    lemniscate constants.  The remaining sets show the dependence
%    on Z when Y = 1 (no loss of generality because of homogeneity)
%    and X = 0.5 (midway between the complete case X = 0 and the
%    degenerate case X = Y).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 May 2018
%
%  Author:
%
%    Original FORTRAN77 version by Bille Carlson, Elaine Notis.
%    This MATLAB version by John Burkardt.
%
  x = [ ...
   0.00E+00, ...
   0.55E-78, ...
   0.00E+00, ...
   0.55E-78, ...
   0.00E+00, ...
   0.55E-78, ...
   0.00E+00, ...
   0.55E-78, ...
   3.01E-51, ...
   3.01E-51, ...
   0.99E+48, ...
   0.99E+48, ...
   0.00E+00, ...
   0.00E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00 ];

  y = [ ...
   6.01E-51, ...
   6.01E-51, ...
   6.01E-51, ...
   6.01E-51, ...
   0.99E+48, ...
   0.99E+48, ...
   0.99E+48, ...
   0.99E+48, ...
   3.01E-51, ...
   3.01E-51, ...
   0.99E+48, ...
   0.99E+48, ...
   3.01E-51, ...
   2.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00 ];

  z = [ ...
   6.01E-51, ...
   6.01E-51, ...
   0.99E+48, ...
   0.99E+48, ...
   6.01E-51, ...
   6.01E-51, ...
   0.99E+48, ...
   0.99E+48, ...
   6.01E-51, ...
   0.99E+48, ...
   6.01E-51, ...
   0.99E+48, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E-10, ...
   1.00E-05, ...
   1.00E-02, ...
   1.00E-01, ...
   2.00E-01, ...
   5.00E-01, ...
   1.00E+00, ...
   2.00E+00, ...
   5.00E+00, ...
   1.00E+01, ...
   1.00E+02, ...
   1.00E+05, ...
   1.00E+10 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'rd_test():\n' );
  fprintf ( 1, '  rd() evaluates the Carlson elliptic integral\n' );
  fprintf ( 1, '  of the second kind, RD(X,Y,Z)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '               X                          Y          ' );
  fprintf ( 1, '                Z                         RD(X,Y,Z)\n' );
  fprintf ( 1, '\n' );

  errtol = 1.0E-3;

  for i = 1 : 27
    [ eliptc, ierr ] = rd ( x(i), y(i), z(i), errtol );
    if (ierr == 0 )
      fprintf ( 1, '%27.16g%27.16g%27.16g%27.16g\n', x(i), y(i), z(i), eliptc );
    else
      fprintf ( 1, '%27.16g%27.16g%27.16g  ***Error***\n', ...
        x(i), y(i), z(i) );
    end
  end

  return
end
