octave_packages/specfun-1.1.0/ellipke.m

   1 ## Copyright (C) 2001 David Billinghurst
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{k}, @var{e}] =} ellipke (@var{m}[,@var{tol}])
  18 ## Compute complete elliptic integral of first K(@var{m}) and second E(@var{m}).
  19 ##
  20 ## @var{m} is either real array or scalar with 0 <= m <= 1
  21 ##
  22 ## @var{tol} will be ignored (@sc{Matlab} uses this to allow faster, less
  23 ## accurate approximation)
  24 ##
  25 ## Ref: Abramowitz, Milton and Stegun, Irene A. Handbook of Mathematical
  26 ## Functions, Dover, 1965, Chapter 17.
  27 ## @seealso{ellipj}
  28 ## @end deftypefn
  29
  30 ## Author: David Billinghurst <David.Billinghurst@riotinto.com>
  31 ## Created: 31 January 2001
  32 ## 2001-02-01 Paul Kienzle
  33 ##   * vectorized
  34 ##   * included function name in error messages
  35 ## 2003-1-18 Jaakko Ruohio
  36 ##   * extended for m < 0
  37
  38 function [k,e] = ellipke( m )
  39
  40   if (nargin < 1 || nargin > 2)
  41     print_usage;
  42   endif
  43
  44   k = e = zeros(size(m));
  45   m = m(:);
  46   if any(~isreal(m))
  47     error("ellipke must have real m");
  48   endif
  49   if any(m>1)
  50     error("ellipke must have m <= 1");
  51   endif
  52
  53   Nmax = 16;
  54   idx = find(m == 1);
  55   if (!isempty(idx))
  56     k(idx) = Inf;
  57     e(idx) = 1.0;
  58   endif
  59
  60   idx = find(m == -Inf);
  61   if (!isempty(idx))
  62     k(idx) = 0.0;
  63     e(idx) = Inf;
  64   endif
  65
  66   ## Arithmetic-Geometric Mean (AGM) algorithm
  67   ## ( Abramowitz and Stegun, Section 17.6 )
  68   idx = find(m != 1 & m != -Inf);
  69   if (!isempty(idx))
  70     idx_neg = find(m < 0 & m != -Inf);
  71     mult_k = 1./sqrt(1-m(idx_neg));
  72     mult_e = sqrt(1-m(idx_neg));
  73     m(idx_neg) = -m(idx_neg)./(1-m(idx_neg));
  74     a = ones(length(idx),1);
  75     b = sqrt(1.0-m(idx));
  76     c = sqrt(m(idx));
  77     f = 0.5;
  78     sum = f*c.*c;
  79     for n = 2:Nmax
  80       t = (a+b)/2;
  81       c = (a-b)/2;
  82       b = sqrt(a.*b);
  83       a = t;
  84       f = f * 2;
  85       sum = sum + f*c.*c;
  86       if all(c./a < eps), break; endif
  87     endfor
  88     if n >= Nmax, error("ellipke: not enough workspace"); endif
  89     k(idx) = 0.5*pi./a;
  90     e(idx) = 0.5*pi.*(1.0-sum)./a;
  91     k(idx_neg) = mult_k.*k(idx_neg);
  92     e(idx_neg) = mult_e.*e(idx_neg);
  93   endif
  94
  95 endfunction
  96
  97 %!test
  98 %! ## Test complete elliptic functions of first and second kind
  99 %! ## against "exact" solution from Mathematica 3.0
 100 %! ##
 101 %! ## David Billinghurst <David.Billinghurst@riotinto.com>
 102 %! ## 1 February 2001
 103 %! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ];
 104 %! [k,e] = ellipke(m);
 105 %!
 106 %! # K(1.0) is really infinity - see below
 107 %! K = [
 108 %!  1.5707963267948966192;
 109 %!  1.5747455615173559527;
 110 %!  1.6124413487202193982;
 111 %!  1.8540746773013719184;
 112 %!  2.5780921133481731882;
 113 %!  3.6956373629898746778;
 114 %!  0.0 ];
 115 %! E = [
 116 %!  1.5707963267948966192;
 117 %!  1.5668619420216682912;
 118 %!  1.5307576368977632025;
 119 %!  1.3506438810476755025;
 120 %!  1.1047747327040733261;
 121 %!  1.0159935450252239356;
 122 %!  1.0 ];
 123 %! if k(7)==Inf, k(7)=0.0; endif;
 124 %! assert(K,k,8*eps);
 125 %! assert(E,e,8*eps);