X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fspecfun-1.1.0%2Fellipke.m;fp=octave_packages%2Fspecfun-1.1.0%2Fellipke.m;h=2f1e73a09dfa30d9f33080de93540e532d4f77ef;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git
diff --git a/octave_packages/specfun-1.1.0/ellipke.m b/octave_packages/specfun-1.1.0/ellipke.m
new file mode 100644
index 0000000..2f1e73a
--- /dev/null
+++ b/octave_packages/specfun-1.1.0/ellipke.m
@@ -0,0 +1,125 @@
+## Copyright (C) 2001 David Billinghurst
+##
+## This program 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.
+##
+## This program 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 this program; If not, see .
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{k}, @var{e}] =} ellipke (@var{m}[,@var{tol}])
+## Compute complete elliptic integral of first K(@var{m}) and second E(@var{m}).
+##
+## @var{m} is either real array or scalar with 0 <= m <= 1
+##
+## @var{tol} will be ignored (@sc{Matlab} uses this to allow faster, less
+## accurate approximation)
+##
+## Ref: Abramowitz, Milton and Stegun, Irene A. Handbook of Mathematical
+## Functions, Dover, 1965, Chapter 17.
+## @seealso{ellipj}
+## @end deftypefn
+
+## Author: David Billinghurst
+## Created: 31 January 2001
+## 2001-02-01 Paul Kienzle
+## * vectorized
+## * included function name in error messages
+## 2003-1-18 Jaakko Ruohio
+## * extended for m < 0
+
+function [k,e] = ellipke( m )
+
+ if (nargin < 1 || nargin > 2)
+ print_usage;
+ endif
+
+ k = e = zeros(size(m));
+ m = m(:);
+ if any(~isreal(m))
+ error("ellipke must have real m");
+ endif
+ if any(m>1)
+ error("ellipke must have m <= 1");
+ endif
+
+ Nmax = 16;
+ idx = find(m == 1);
+ if (!isempty(idx))
+ k(idx) = Inf;
+ e(idx) = 1.0;
+ endif
+
+ idx = find(m == -Inf);
+ if (!isempty(idx))
+ k(idx) = 0.0;
+ e(idx) = Inf;
+ endif
+
+ ## Arithmetic-Geometric Mean (AGM) algorithm
+ ## ( Abramowitz and Stegun, Section 17.6 )
+ idx = find(m != 1 & m != -Inf);
+ if (!isempty(idx))
+ idx_neg = find(m < 0 & m != -Inf);
+ 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(length(idx),1);
+ b = sqrt(1.0-m(idx));
+ c = sqrt(m(idx));
+ f = 0.5;
+ sum = f*c.*c;
+ for n = 2:Nmax
+ t = (a+b)/2;
+ c = (a-b)/2;
+ b = sqrt(a.*b);
+ a = t;
+ f = f * 2;
+ sum = sum + f*c.*c;
+ if all(c./a < eps), break; endif
+ endfor
+ if n >= Nmax, error("ellipke: not enough workspace"); endif
+ k(idx) = 0.5*pi./a;
+ e(idx) = 0.5*pi.*(1.0-sum)./a;
+ k(idx_neg) = mult_k.*k(idx_neg);
+ e(idx_neg) = mult_e.*e(idx_neg);
+ endif
+
+endfunction
+
+%!test
+%! ## Test complete elliptic functions of first and second kind
+%! ## against "exact" solution from Mathematica 3.0
+%! ##
+%! ## David Billinghurst
+%! ## 1 February 2001
+%! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ];
+%! [k,e] = ellipke(m);
+%!
+%! # K(1.0) is really infinity - see below
+%! K = [
+%! 1.5707963267948966192;
+%! 1.5747455615173559527;
+%! 1.6124413487202193982;
+%! 1.8540746773013719184;
+%! 2.5780921133481731882;
+%! 3.6956373629898746778;
+%! 0.0 ];
+%! E = [
+%! 1.5707963267948966192;
+%! 1.5668619420216682912;
+%! 1.5307576368977632025;
+%! 1.3506438810476755025;
+%! 1.1047747327040733261;
+%! 1.0159935450252239356;
+%! 1.0 ];
+%! if k(7)==Inf, k(7)=0.0; endif;
+%! assert(K,k,8*eps);
+%! assert(E,e,8*eps);