function [ value, ierr ] = rj ( x, y, z, p, errtol )

%*****************************************************************************80
%
%% rj() computes an incomplete elliptic integral of the third kind, RJ(X,Y,Z,P).
%
%  Discussion:
%
%    This function computes an incomplete elliptic integral of the third kind.
%
%    RJ(X,Y,Z,P) = Integral ( 0 <= T < oo )
%
%                  -1/2     -1/2     -1/2     -1
%        (3/2)(T+X)    (T+Y)    (T+Z)    (T+P)  DT,
%
%    where X, Y, and Z are nonnegative, at most one of them is
%    zero, and P is positive.
%
%    If X or Y or Z is zero, then the integral is complete.
%
%    The duplication theorem is iterated until the variables are nearly equal,
%    and the function is then expanded in Taylor series to fifth order.
%
%    Check by addition theorem:
%
%      RJ(X,X+Z,X+W,X+P)
%      + RJ(Y,Y+Z,Y+W,Y+P) + (A-B) * RJ(A,B,B,A) + 3 / sqrt ( A)
%      = RJ(0,Z,W,P), where X,Y,Z,W,P are positive and X * Y
%      = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
%      and B - A = P * (P-Z) * (P-W).
%
%    The sum of the third and fourth terms on the left side is 3 * RC(A,B).
%
%  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.
%
%  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.
%
%  Input:
%
%    real X, Y, Z, P, the arguments in the integral.
%
%    real ERRTOL, the error tolerance.
%    Relative error due to truncation of the series for rj
%    is less than 3 * ERRTOL ^ 6 / (1 - ERRTOL) ^ 3/2.
%    An error tolerance (ETOLRC) will be passed to the subroutine
%    for RC to make the truncation error for RC less than for RJ.
%    Sample choices:
%      ERRTOL   Relative truncation error less than
%      1.E-3    4.E-18
%      3.E-3    3.E-15
%      1.E-2    4.E-12
%      3.E-2    3.E-9
%      1.E-1    4.E-6
%
%  Output:
%
%    real VALUE, the value of the function.
%
%    integer IERR, the error flag.
%    0, no error occurred.
%    1, abnormal termination.
%
%  Local:
%
%    real LOLIM, UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
%    LOLIM IS NOT LESS THAN THE CUBE ROOT OF THE VALUE
%    OF LOLIM USED IN THE SUBROUTINE FOR RC.
%    UPLIM IS NOT GREATER THAN 0.3 TIMES THE CUBE ROOT OF
%    THE VALUE OF UPLIM USED IN THE SUBROUTINE FOR RC.
%
  lolim = 2.E-26;
  uplim = 3.E+24;

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

  ierr = 0;
  xn = x;
  yn = y;
  zn = z;
  pn = p;
  sigma = 0.0;
  power4 = 1.0;
  etolrc = 0.5 * errtol;

  while ( true )

    mu = ( xn + yn + zn + pn + pn ) * 0.2;
    xndev = ( mu - xn ) / mu;
    yndev = ( mu - yn ) / mu;
    zndev = ( mu - zn ) / mu;
    pndev = ( mu - pn ) / mu;
    epslon = max ( abs ( xndev ), ...
     max ( abs ( yndev ), ...
     max ( abs ( zndev ), abs ( pndev ) ) ) );

    if ( epslon < errtol )
      c1 = 3.0 / 14.0;
      c2 = 1.0 / 3.0;
      c3 = 3.0 / 22.0;
      c4 = 3.0 / 26.0;
      ea = xndev * ( yndev + zndev ) + yndev * zndev;
      eb = xndev * yndev * zndev;
      ec = pndev * pndev;
      e2 = ea - 3.0 * ec;
      e3 = eb + 2.0 * pndev * ( ea - ec );
      s1 = 1.0 + e2 * ( - c1 + 0.75 * c3 * e2 - 1.5 * c4 * e3 );
      s2 = eb * ( 0.5 * c2 + pndev * ( - c3 - c3 + pndev * c4 ) );
      s3 = pndev * ea * ( c2 - pndev * c3 ) - c2 * pndev * ec;
      value = 3.0 * sigma + power4 * ( s1 + s2 + s3 ) / ( mu * sqrt ( mu ) );
      return
    end

    xnroot = sqrt ( xn );
    ynroot = sqrt ( yn );
    znroot = sqrt ( zn );
    lamda = xnroot * ( ynroot + znroot ) + ynroot * znroot;
    alfa = pn * ( xnroot + ynroot + znroot ) + xnroot * ynroot * znroot;
    alfa = alfa * alfa;
    beta = pn * ( pn + lamda ) * ( pn + lamda );
    [ value, ierr ] = rc ( alfa, beta, etolrc );
    sigma = sigma + power4 * value;

    if ( ierr ~= 0 )
      value = 0.0;
      return
    end

    power4 = power4 * 0.25;
    xn = ( xn + lamda ) * 0.25;
    yn = ( yn + lamda ) * 0.25;
    zn = ( zn + lamda ) * 0.25;
    pn = ( pn + lamda ) * 0.25;

  end

end
