function rf_test ( )

%*****************************************************************************80
%
%% rf_test() tests rf().
%
%  Discussion:
%
%    This driver tests the function for the
%    integral RF(X,Y,Z), which is symmetric in X, Y, Z.  The first nine
%    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
%    tenth set tests the failure exit.  The eleventh set is a check
%    value: RF(0,1,2) = A, where A is the first lemniscate constant.
%    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 = [ ...
   1.51E-78, ...
   1.51E-78, ...
   0.00E+00, ...
   0.00E+00, ...
   0.00E+00, ...
   0.99E+75, ...
   0.55E-78, ...
   0.55E-78, ...
   0.55E-78, ...
   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, ...
   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, ...
   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, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00, ...
   0.50E+00 ];

  y = [ ...
   1.51E-78, ...
   1.51E-78, ...
   3.01E-78, ...
   3.01E-78, ...
   0.99E+75, ...
   0.99E+75, ...
   3.01E-78, ...
   3.01E-78, ...
   0.99E+75, ...
   2.00E-78, ...
   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, ...
   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, ...
   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, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00 ];

  z = [ ...
   1.51E-78, ...
   0.99E+75, ...
   3.01E-78, ...
   0.99E+75, ...
   0.99E+75, ...
   0.99E+75, ...
   3.01E-78, ...
   0.99E+75, ...
   0.99E+75, ...
   1.00E+00, ...
   2.00E+00, ...
   1.00E+00, ...
   1.10E+00, ...
   1.20E+00, ...
   1.30E+00, ...
   1.40E+00, ...
   1.50E+00, ...
   1.60E+00, ...
   1.70E+00, ...
   1.80E+00, ...
   1.90E+00, ...
   2.00E+00, ...
   2.20E+00, ...
   2.40E+00, ...
   2.60E+00, ...
   2.80E+00, ...
   3.00E+00, ...
   3.50E+00, ...
   4.00E+00, ...
   4.50E+00, ...
   5.00E+00, ...
   6.00E+00, ...
   7.00E+00, ...
   8.00E+00, ...
   9.00E+00, ...
   1.00E+01, ...
   2.00E+01, ...
   3.00E+01, ...
   4.00E+01, ...
   5.00E+01, ...
   1.00E+02, ...
   2.00E+02, ...
   5.00E+02, ...
   1.00E+03, ...
   1.00E+04, ...
   1.00E+05, ...
   1.00E+06, ...
   1.00E+08, ...
   1.00E+10, ...
   1.00E+12, ...
   1.00E+15, ...
   1.00E+20, ...
   1.00E+30, ...
   1.00E+40, ...
   1.00E+50 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'RF_TEST\n' );
  fprintf ( 1, '  RF evaluates the Carlson elliptic integral\n' );
  fprintf ( 1, '  of the first kind, RF(X,Y,Z)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '               X                          Y          ' );
  fprintf ( 1, '                Z                         RF(X,Y,Z)\n' );
  fprintf ( 1, '\n' );

  errtol = 1.0E-3;

  for i = 1 : 55
    [ eliptc, ierr ] = rf ( 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
