function rj_test ( )

%*****************************************************************************80
%
%% rj_test() tests rj().
%
%  Discussion:
%
%    This driver tests the real ( kind = 8 ) function subroutine for the
%    integral Rj(X,Y,Z,P), which is symmetric in X, Y, Z.  The first
%    twenty sets of values of X, Y, Z, P 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, P, 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 twenty-first set tests the failure exit.  The twenty-
%    second set is a check value:
%      RJ(2,3,4,5) = 0.1429757966715675383323308.
%    The remaining sets show the dependence on Z and P
%    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.
%
  p = [ ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   1.00E+00, ...
   5.00E+00, ...
   0.25E+00, ...
   0.75E+00, ...
   1.00E+00, ...
   2.00E+00, ...
   0.25E+00, ...
   0.75E+00, ...
   1.50E+00, ...
   4.00E+00, ...
   0.25E+00, ...
   0.75E+00, ...
   3.00E+00, ...
   1.00E+01, ...
   0.25E+00, ...
   0.75E+00, ...
   5.00E+00, ...
   2.00E+01, ...
   0.25E+00, ...
   0.75E+00, ...
   5.00E+01, ...
   2.00E+02 ];

  x = [ ...
   1.01E-26, ...
   1.01E-26, ...
   0.00E+00, ...
   0.00E+00, ...
   0.00E+00, ...
   2.99E+24, ...
   0.55E-78, ...
   0.55E-78, ...
   0.55E-78, ...
   2.01E-26, ...
   1.01E-26, ...
   1.01E-26, ...
   0.00E+00, ...
   0.00E+00, ...
   0.00E+00, ...
   2.99E+24, ...
   0.55E-78, ...
   0.55E-78, ...
   0.55E-78, ...
   2.01E-26, ...
   0.00E+00, ...
   2.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 ];

  y = [ ...
   1.01E-26, ...
   1.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   2.01E-26, ...
   2.99E+24, ...
   2.01E-26, ...
   1.01E-26, ...
   1.01E-26, ...
   2.01E-26, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   2.01E-26, ...
   2.99E+24, ...
   2.01E-26, ...
   1.90E-26, ...
   3.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.01E-26, ...
   2.99E+24, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   1.01E-26, ...
   2.99E+24, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   2.99E+24, ...
   2.99E+24, ...
   2.01E-26, ...
   1.90E-26, ...
   4.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   1.00E+00, ...
   2.00E+00, ...
   2.00E+00, ...
   2.00E+00, ...
   2.00E+00, ...
   5.00E+00, ...
   5.00E+00, ...
   5.00E+00, ...
   5.00E+00, ...
   1.00E+01, ...
   1.00E+01, ...
   1.00E+01, ...
   1.00E+01, ...
   1.00E+02, ...
   1.00E+02, ...
   1.00E+02, ...
   1.00E+02 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'rj_test():\n' );
  fprintf ( 1, '  rj() evaluates the Carlson elliptic integral\n' );
  fprintf ( 1, '  of the third kind, RJ(X,Y,Z,P)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '               X                          Y          ' );
  fprintf ( 1, '                Z                         P' );
  fprintf ( 1, '                RJ(X,Y,Z,P)\n' );
  fprintf ( 1, '\n' );

  errtol = 1.0E-3;

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

  return
end
