## Copyright (C) 2001, 2013 David Billinghurst
## Copyright (C) 2001, 2013 Paul Kienzle
## Copyright (C) 2003, 2013 Jaakko Ruohio
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{k} =} ellipke (@var{m})
## @deftypefnx {Function File} {@var{k} =} ellipke (@var{m}, @var{tol})
## @deftypefnx {Function File} {[@var{k}, @var{e}] =} ellipke (@dots{})
## Compute complete elliptic integrals of the first K(@var{m}) and second
## E(@var{m}) kind.
##
## @var{m} must be a scalar or real array with -Inf @leq{} @var{m} @leq{} 1.
##
## The optional input @var{tol} is currently ignored (@sc{matlab} uses this
## to allow a faster, less accurate approximation).
##
## Called with only one output, elliptic integrals of the first kind are
## returned.
##
## Reference: Milton @nospell{Abramowitz} and Irene A. @nospell{Stegun},
## @cite{Handbook of Mathematical Functions}, Chapter 17, Dover, 1965.
## @seealso{ellipj}
## @end deftypefn

## Author: David Billinghurst <David.Billinghurst@riotinto.com>
## Author: Paul Kienzle <pkienzle@users.sf.net>
## Author: Jaakko Ruohio

function [k, e] = ellipke (m)

  if (nargin < 1 || nargin > 2)
    print_usage ();
  endif

  m = m(:);
  if (! isreal (m))
    error ("ellipke: M must be real");
  elseif (any (m > 1))
    error ("ellipke: M must be <= 1");
  endif

  k = e = zeros (size (m));

  ## Handle extreme values
  idx_1 = (m == 1);
  k(idx_1) = Inf;
  e(idx_1) = 1;

  idx_neginf = (m == -Inf);
  k(idx_neginf) = 0;
  e(idx_neginf) = Inf;

  ## Arithmetic-Geometric Mean (AGM) algorithm
  ## ( Abramowitz and Stegun, Section 17.6 )
  Nmax = 16;
  idx = !idx_1 & !idx_neginf;
  if (any (idx))
    idx_neg = find (m < 0 & !idx_neginf);
    mult_k = 1./sqrt (1 - m(idx_neg));
    mult_e = sqrt (1 - m(idx_neg));
    m(idx_neg) = -m(idx_neg) ./ (1 - m(idx_neg));
    a = ones (sum (idx), 1);
    b = sqrt (1 - m(idx));
    c = sqrt (m(idx));
    f = 0.5;
    sum = f*c.^2;
    n = 2;
    do
      t = (a + b)/2;
      c = (a - b)/2;
      b = sqrt (a.*b);
      a = t;
      f *= 2;
      sum += f*c.^2;
    until (all (c./a < eps) || (++n > Nmax))
    if (n >= Nmax)
      error ("ellipke: algorithm did not converge in %d iterations", Nmax);
    endif
    k(idx) = 0.5*pi ./ a;
    e(idx) = 0.5*pi*(1 - sum) ./ a;
    k(idx_neg) = mult_k .* k(idx_neg);
    e(idx_neg) = mult_e .* e(idx_neg);
  endif

endfunction


## Test complete elliptic functions of first and second kind
## against "exact" solution from Mathematica 3.0
%!test
%! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0];
%! [k,e] = ellipke (m);
%!
%! k_exp = [1.5707963267948966192;
%!          1.5747455615173559527;
%!          1.6124413487202193982;
%!          1.8540746773013719184;
%!          2.5780921133481731882;
%!          3.6956373629898746778;
%!          Inf ];
%! e_exp = [1.5707963267948966192;
%!          1.5668619420216682912;
%!          1.5307576368977632025;
%!          1.3506438810476755025;
%!          1.1047747327040733261;
%!          1.0159935450252239356;
%!          1.0 ];
%! assert (k, k_exp, 8*eps);
%! assert (e, e_exp, 8*eps);

## Test against A&S Table 17.1
%!test
%! m = [0:5:50]'/100;
%! k_exp = [1.570796326794897;
%!          1.591003453790792;
%!          1.612441348720219;
%!          1.635256732264580;
%!          1.659623598610528;
%!          1.685750354812596;
%!          1.713889448178791;
%!          1.744350597225613;
%!          1.777519371491253;
%!          1.813883936816983;
%!          1.854074677301372 ];
%! e_exp = [1.570796327;
%!          1.550973352;
%!          1.530757637;
%!          1.510121831;
%!          1.489035058;
%!          1.467462209;
%!          1.445363064;
%!          1.422691133;
%!          1.399392139;
%!          1.375401972;
%!          1.350643881 ];
%! [k,e] = ellipke (m);
%! assert (k, k_exp, 1e-15);
%! assert (e, e_exp, 1e-8);

## Test negative values against "exact" solution from Mathematica.
%! m = [-0.01; -1; -5; -100; -1000; -Inf];
%! [k,e] = ellipke (m);
%!
%! k_exp = [1.5668912730681963584;
%!          1.3110287771460599052;
%!          0.9555039270640439337;
%!          0.3682192486091410329;
%!          0.1530293349884987857;
%!          0];
%! e_exp = [1.5747159850169884130;
%!          1.9100988945138560089;
%!          2.8301982463458773125;
%!          10.209260919814572009;
%!          31.707204053711259719;
%!          Inf ];
%! assert (k, k_exp, 8*eps);
%! assert (e, e_exp, 8*eps (e_exp));

## Test input validation
%!error ellipke ()
%!error ellipke (1,2,3)