function value = elliptic_inc_pik ( phi, n, k )

%*****************************************************************************80
%
%% ELLIPTIC_INC_PIK evaluates the incomplete elliptic integral Pi(PHI,N,K).
%
%  Discussion:
%
%    This is one form of what is sometimes called the incomplete elliptic
%    integral of the third kind.
%
%    The function is defined by the formula:
%
%      Pi(PHI,N,K) = integral ( 0 <= T <= PHI )
%        dT / (1 - N sin^2(T) ) sqrt ( 1 - k^2 * sin ( T )^2 )
%
%    In MATLAB, the function can be evaluated by:
%
%      ellipticPi(n,phi,k^2)
%
%    The value is computed using Carlson elliptic integrals:
%
%      Pi(phi,n,k) = sin(phi) RF ( cos^2(phi), 1 - k^2 sin^2(phi), 1 ) 
%        + 1/3 n sin^3(phi) RJ ( cos^2(phi), 1 - k^2 sin^2(phi), 1, 1 - n sin^2(phi) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PHI, N, K, the arguments.
%    0 <= PHI <= PI/2.
%         N   <= 1 / sin^2 ( PHI )
%    0 <= K^2 <= 1 / sin^2 ( PHI )
%
%    Output, real VALUE, the function value.
%
  sp = sin ( phi );
  cp = cos ( phi );

  x = cp ^ 2;
  y = ( 1.0 - k * sp ) * ( 1.0 + k * sp );
  z = 1.0;
  p = 1.0 - n * sp ^ 2;
  errtol = 1.0E-03;

  [ value1, ierr ] = rf ( x, y, z, errtol );
  if ( ierr ~= 0 )
    fprintf ( 1, 'ELLIPTIC_INC_PIK - Fatal error!\n' );
    fprintf ( 1, '  RF returned IERR = %d\n', ierr );
    error ( 'ELLIPTIC_INC_PIK - Fatal error!' );
    return
  end

  [ value2, ierr ] = rj ( x, y, z, p, errtol );
  if ( ierr ~= 0 )
    fprintf ( 1, 'ELLIPTIC_INC_PIK - Fatal error!\n' );
    fprintf ( 1, '  RJ returned IERR = %d\n', ierr );
    error ( 'ELLIPTIC_INC_PIK - Fatal error!' );
    return
  end

  value = sp * value1 + n * sp ^ 3 * value2 / 3.0;

  return
end
