function [ value, ierr ] = rc ( x, y, errtol )

%*****************************************************************************80
%
%% RC computes the elementary integral RC(X,Y).
%
%  Discussion:
%
%    This function computes the elementary integral
%
%      RC(X,Y) = Integral ( 0 <= T < oo )
%
%                  -1/2     -1
%        (1/2)(T+X)    (T+Y)  DT,
%
%    where X is nonnegative and Y is positive.  The duplication
%    theorem is iterated until the variables are nearly equal,
%    and the function is then expanded in Taylor series to fifth
%    order.
%
%    Logarithmic, inverse circular, and inverse hyperbolic
%    functions can be expressed in terms of RC.
%
%    Check by addition theorem:
%
%      RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z),
%      where X, Y, and Z are positive and X * Y = Z * Z.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2018
%
%  Author:
%
%    Original FORTRAN77 version by Bille Carlson, Elaine Notis.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Bille Carlson,
%    Computing Elliptic Integrals by Duplication,
%    Numerische Mathematik,
%    Volume 33, 1979, pages 1-16.
%
%    Bille Carlson, Elaine Notis,
%    Algorithm 577, Algorithms for Incomplete Elliptic Integrals,
%    ACM Transactions on Mathematical Software,
%    Volume 7, Number 3, pages 398-403, September 1981.
%
%  Parameters:
%
%    Input, real X, Y, the arguments in the integral.
%
%    Input, real ERRTOL, the error tolerance.
%    Relative error due to truncation is less than
%      16 * ERRTOL ^ 6 / (1 - 2 * ERRTOL).
%    Sample choices:
%      ERRTOL   Relative truncation error less than
%      1.0E-3    2.0E-17
%      3.0E-3    2.0E-14
%      1.0E-2    2.0E-11
%      3.0E-2    2.0E-8
%      1.0E-1    2.0E-5
%
%    Output, real VALUE, the value of the function.
%
%    Output, integer IERR, the error flag.
%    0, no error occurred.
%    1, abnormal termination.
%
%  Local Parameters:
%
%    LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
%    LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
%    UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
%
  lolim = 3.E-78;
  uplim = 1.0E+75;

  if ( ...
    x < 0.0 || ...
    y <= 0.0 || ...
    ( x + y ) < lolim || ...
    uplim < x || ...
    uplim < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'RC - Error!\n' );
    fprintf ( 1, '  Invalid input arguments.\n' );
    fprintf ( 1, '  X = %g\n', x );
    fprintf ( 1, '  Y = %g\n', y );
    fprintf ( 1, '\n' );
    ierr = 1;
    value = 0.0;
    return
  end

  ierr = 0;
  xn = x;
  yn = y;

  while ( true )

    mu = ( xn + yn + yn ) / 3.0;
    sn = ( yn + mu ) / mu - 2.0;

    if ( abs ( sn ) < errtol )
      c1 = 1.0 / 7.0;
      c2 = 9.0 / 22.0;
      s = sn * sn * ( 0.3 + sn * ( c1 + sn * ( 0.375 + sn * c2 ) ) );
      value = ( 1.0 + s ) / sqrt ( mu );
      return
    end

    lamda = 2.0 * sqrt ( xn ) * sqrt ( yn ) + yn;
    xn = ( xn + lamda ) * 0.25;
    yn = ( yn + lamda ) * 0.25;

  end

end
