home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / ELLIPJ.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.3 KB  |  86 lines
  1. function [sn,cn,dn] = ellipj(u,m)
  2. %ELLIPJ    Jacobi elliptic functions SN, CN and DN.
  3. %    [Sn,Cn,Dn] = ELLIPJ(U,M) returns the values of the Jacobi
  4. %    elliptic functions SN, CN and DN, evaluated at argument U
  5. %    and parameter M.  As currently implemented, M is limited
  6. %    to 0 < M < 1.
  7. %    ELLIPJ(U,M) are accurate to EPS.
  8. %    U and M must either be matrices of the same size or either
  9. %    one may be a scalar.
  10. %
  11. %    Be sure you don't confuse the modulus K with the parameter M -
  12. %    they are related in the following way:  M = K^2
  13.  
  14. %    L. Shure 1-9-88
  15. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  16.  
  17. %    ELLIPJ uses the method of the arithmetic-geometric mean
  18. %    described in [1].
  19.  
  20. % References:
  21. %    [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical
  22. %        Functions" Dover Publications", 1965, Ch. 16-17.6.
  23.  
  24. tol = eps(1);
  25. if any(any(imag(u))) | any(any(imag(m)))
  26.     error('Input arguments must be real.')
  27. end
  28. [mm,nm] = size(m);
  29. [mu,nu] = size(u);
  30. if mu ~= mm | nu ~= nm
  31.     if mu*nu == 1 
  32.       u = u*ones(mm,nm);
  33.     elseif mm*nm == 1 
  34.       m = m*ones(mu,nu);
  35.     else
  36.       error('Matrix dimensions must agree.')
  37.     end
  38. end
  39. [ms,ns] = size(u);
  40. mmax = ms*ns;
  41. cn = zeros(ms,ns);
  42. sn = cn;
  43. dn = sn;
  44. m = m(:).';    % make a row vector
  45. u = u(:).';
  46. % pre-allocate space and augment if needed
  47. chunk = 10;
  48. a = zeros(chunk,mmax);
  49. c = a;
  50. b = a;
  51. a(1,:) = ones(1,mmax);
  52. c(1,:) = sqrt(m);
  53. b(1,:) = sqrt(1-m);
  54. n = zeros(1,mmax);
  55. i = 1;
  56. while any(abs(c(i,:)) > tol)
  57.     i = i + 1;
  58.     if i > size(a,1)
  59.       a = [a; zeros(chunk,mmax)];
  60.       b = [b; zeros(chunk,mmax)];
  61.       c = [c; zeros(chunk,mmax)];
  62.     end
  63.     a(i,:) = 0.5 * (a(i-1,:) + b(i-1,:));
  64.     b(i,:) = sqrt(a(i-1,:) .* b(i-1,:));
  65.     c(i,:) = 0.5 * (a(i-1,:) - b(i-1,:));
  66.     in = find((abs(c(i,:)) <= tol) & (abs(c(i-1,:)) > tol));
  67.     if ~isempty(in)
  68.       [mi,ni] = size(in);
  69.       n(in) = ones(mi,ni)*(i-1);
  70.     end
  71. end
  72. phin = zeros(i,mmax);
  73. phin(i,:) = (2 .^ n).*a(i,:).*u;
  74. while i > 1
  75.     i = i - 1;
  76.     in = find(n >= i);
  77.     phin(i,:) = phin(i+1,:);
  78.     if ~isempty(in)
  79.       phin(i,in) = 0.5 * ...
  80.       (asin(c(i+1,in).*sin(rem(phin(i+1,in),2*pi))./a(i+1,in)) + phin(i+1,in));
  81.     end
  82. end
  83. sn(:) = sin(rem(phin(1,:),2*pi));
  84. cn(:) = cos(rem(phin(1,:),2*pi));
  85. dn(:) = sqrt(1 - m .* (sn(:).').^2);
  86.