home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / ELLIPKE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.2 KB  |  52 lines
  1. function [k,e] = ellipke(m)
  2. %ELLIPKE Complete elliptic integrals.
  3. %    [K,E] = ELLIPKE(M) returns the value of the complete elliptic
  4. %    integrals of the first and second kinds, evaluated at M.  As currently
  5. %    implemented, M is limited to 0 < M < 1.
  6. %    The accuracy of ELLIPKE(M) is EPS.
  7. %
  8. %    Be sure you don't confuse the modulus K with the parameter M -
  9. %    they are related in the following way:  M = K^2
  10.  
  11. %    L. Shure 1-9-88
  12. %    Modified to include the second kind by Bjorn Bonnevier
  13. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  14.  
  15. %   ELLIPKE uses the method of the arithmetic-geometric mean
  16. %   described in [1].
  17.  
  18. %   References:
  19. %   [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical
  20. %       Functions" Dover Publications", 1965, 17.6.
  21.  
  22. if any(any(imag(m)))
  23.     error('Input arguments must be real.')
  24. end
  25. a0 = 1;
  26. b0 = sqrt(1-m);
  27. s0 = m;
  28. i1 = 0; mm = 1;
  29. while mm > eps
  30.     a1 = (a0+b0)/2;
  31.     b1 = sqrt(a0.*b0);
  32.     c1 = (a0-b0)/2;
  33.     i1 = i1 + 1;
  34.     w1 = 2^i1*c1.^2;
  35.     mm = max(max(w1));
  36.     s0 = s0 + w1;
  37.     a0 = a1;
  38.     b0 = b1;
  39. end
  40. k = pi./(2*a1);
  41. e = k.*(1-s0/2);
  42. im = find(m ==1);
  43. if ~isempty(im)
  44.     if isieee
  45.         e(im) = ones(max(size(im)),1);
  46.         k(im) = e(im)*inf;
  47.     else
  48.         error('Result is infinite.')
  49.     end
  50. end
  51.  
  52.