Re: elliptic integrals?

[Paul Kienzle]

The cephes math library contains C code for complete and incomplete
elliptic integrals of the first and second kind and the
Jacobian elliptic functions.  See <http://netlib.org/cephes/ellf.shar>.

>From these, I've created matlab compatible ellipj.m and ellipke.m
(I really ought to learn more about .oct files).  I will add them
to signalPAK soon, but until I do, here they are.  Note that I have
not translated the incomplete elliptic integral functions since I
won't need them and since they aren't listed in the matlab documentation
that I have access to, but if you do make them callable from octave,
pass them on to me and I will include them in my bundle.

Paul Kienzle
address@hidden

PS, Since I haven't yet written/translated the elliptic filter routines
which use these functions, I haven't extensively tested them.  I did at
one point verify that they returned the same thing as the original C, so
they should be okay.  If you find any problems with them, let me know!

ellipj.m:

## usage [sn cn dn phi] = ellipj(u, m)
##
## Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
## and dn(u|m) of parameter m between 0 and 1, and real
## argument u.  If u and m are vectors, it returns a vector of
## the same shape.
##
## These functions are periodic, with quarter-period on the
## real axis equal to the complete elliptic integral
## ellpk(1.0-m).
##
## Relation to incomplete elliptic integral:
## If u = ellipke(phi,m), then sn(u|m) = sin(phi),
## and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
##
## Computation is by means of the arithmetic-geometric mean
## algorithm, except when m is within 1e-9 of 0 or 1.  In the
## latter case with m close to 1, the approximation applies
## only for phi < pi/2.
##
## Example
##    m = 0:0.005:1; u=[-10:0.1:10]';
##    M = kron(m,ones(size(u))); U = kron(ones(size(m)),u);
##    [sn cn dn phi] = ellipj(U,M);
##    imagesc(sn);

## Cephes Math Library, Release 2.0:  April, 1987
## Copyright 1984, 1987 by Stephen L. Moshier
## Direct inquiries to 30 Frost Street, Cambridge, MA 02140
##
##   Some software in this archive may be from the book _Methods and
## Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
## International, 1989) or from the Cephes Mathematical Library, a
## commercial product. In either event, it is copyrighted by the author.
## What you see here may be used freely but it comes with no support or
## guarantee.
##
##   Stephen L. Moshier
##   address@hidden

## Adapted-By: Paul Kienzle

function [sn cn dn phi]=ellipj(u,m)

  if nargin != 2, usage("[sn cn dn]=ellipj(u,m)"); endif
  if rows(u) != rows(m) || columns(u) != columns(m), 
    error("ellipj must have same shape for u and m");
  endif
  if any(m < 0.0 | m > 1.0), 
    error("ellipj must have m in the range [0,1]");
  endif

  [nr nc] = size(u);
  u = u(:); m = m(:);
  sn=cn=dn=phi=zeros(size(u));

  ## Check for special cases

  idx = find(m<1.0e-9);
  if !isempty(idx)
    t = sin(u(idx));
    b = cos(u(idx));
    ai = 0.25 * m(idx) .* (u(idx) - t.*b);
    sn(idx) = t - ai.*b;
    cn(idx) = b + ai.*t;
    phi(idx) = u(idx) - ai;
    dn(idx) = 1.0 - 0.5*m(idx).*t.*t;
  endif

  idx = find( m >= 0.9999999999 );
  if !isempty(idx)
    ai = 0.25 * (1.0-m(idx));
    b = cosh(u(idx));
    t = tanh(u(idx));
    s = 1.0./b;
    twon = b .* sinh(u(idx));
    sn(idx) = t + ai .* (twon - u(idx))./(b.*b);
    phi(idx) = 2.0*atan(exp(u(idx))) - pi/2 + ai.*(twon - u(idx))./b;
    ai = ai .* t .* s;
    cn(idx) = s - ai .* (twon - u(idx));
    dn(idx) = s + ai .* (twon + u(idx));
  endif

  idx = find( m >= 1e-9 & m < 0.9999999999 );
  if !isempty(idx)
    ## A. G. M. scale
    a = c = zeros(length(idx),9);
    a(:,1) = ones(length(idx),1);
    b = sqrt(1.0 - m(idx));
    c(:,1) = sqrt(m(idx));
    twon = 1;
    i = 1;
    
    while( any(abs(c(:,i)./a(:,i)) > eps ))
      if( i > 8 ), 
        warning("ellipj overflow");
        break;
      endif
      ai = a(:,i);
      i=i+1;
      c(:,i) = ( ai - b )/2.0;
      t = sqrt( ai .* b );
      a(:,i) = ( ai + b )/2.0;
      b = t;
      twon = twon*2;
    endwhile
    
    
    ## backward recurrence
    phi(idx) = twon * a(:,i) .* u(idx);
    for j=i:-1:2
      t = c(:,j) .* sin(phi(idx)) ./ a(:,j);
      b = phi(idx);
      phi(idx) = (asin(t) + phi(idx))/2.0;
    endfor

    sn(idx) = sin(phi(idx));
    cn(idx) = cos(phi(idx));
    dn(idx) = cn(idx)./cos(phi(idx)-b);
  endif
  sn=reshape(sn,nr,nc);


ellipke.m:

## E(0) = K(0) = pi/2.
## E(1) = 1;
##
## If m or x is a matrix, matrices of the same shape are returned.
##
## Example
##    [K E] = ellipke(0:0.01:0.99);
##    plot(K, ";K;", E, ";E;");

## Cephes Math Library, Release 2.0:  April, 1987
## Copyright 1984, 1987 by Stephen L. Moshier
## Direct inquiries to 30 Frost Street, Cambridge, MA 02140
##
##   Some software in this archive may be from the book _Methods and
## Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
## International, 1989) or from the Cephes Mathematical Library, a
## commercial product. In either event, it is copyrighted by the author.
## What you see here may be used freely but it comes with no support or
## guarantee.
##
##   Stephen L. Moshier
##   address@hidden

## Adapted-By: Paul Kienzle

function [K E]=ellipke(m,x)

  if nargin<1 || nargin>2, usage("[k1 k2] = ellipke(m)"); endif
  if any(m < 0.0 | m > 1.0), error("ellipke out of range [0,1]"); endif
  if nargin==1, x = 1-m; endif

  [nr nc]=size(x);
  x = x(:);
  K = zeros(size(x));
  idx = find(x>eps);
  if !isempty(idx)
    P = [1.37982864606273237150E-4,
         2.28025724005875567385E-3,
         7.97404013220415179367E-3,
         9.85821379021226008714E-3,
         6.87489687449949877925E-3,
         6.18901033637687613229E-3,
         8.79078273952743772254E-3,
         1.49380448916805252718E-2,
         3.08851465246711995998E-2,
         9.65735902811690126535E-2,
         1.38629436111989062502E0 ];
    Q = [2.94078955048598507511E-5,
         9.14184723865917226571E-4,
         5.94058303753167793257E-3,
         1.54850516649762399335E-2,
         2.39089602715924892727E-2,
         3.01204715227604046988E-2,
         3.73774314173823228969E-2,
         4.88280347570998239232E-2,
         7.03124996963957469739E-2,
         1.24999999999870820058E-1,
         4.99999999999999999821E-1 ];
    K(idx) = polyval(P,x(idx)) - log(x(idx)) .* polyval(Q,x(idx));
  endif

  idx = find(x==0);
  if !isempty(idx), 
    K(idx) = inf*ones(size(idx)); 
  endif

  idx = find(x>0 & x<=eps);
  if !isempty(idx)
    K(idx) = 1.3862943611198906188E0 - 0.5*log(x(idx)); # log(4) - log(x)/2
  endif

  K = reshape(K,nr,nc);
  if nargout == 1, return; endif

  E = zeros(size(x));

  idx = find(x==0);
  if !isempty(idx), 
    E(idx) = ones(size(idx)); 
  endif

  idx = find(x!=0);
  if !isempty(idx)
    P = [1.53552577301013293365E-4,
         2.50888492163602060990E-3,
         8.68786816565889628429E-3,
         1.07350949056076193403E-2,
         7.77395492516787092951E-3,
         7.58395289413514708519E-3,
         1.15688436810574127319E-2,
         2.18317996015557253103E-2,
         5.68051945617860553470E-2,
         4.43147180560990850618E-1,
         1.00000000000000000299E0];
    Q = [3.27954898576485872656E-5,
         1.00962792679356715133E-3,
         6.50609489976927491433E-3,
         1.68862163993311317300E-2,
         2.61769742454493659583E-2,
         3.34833904888224918614E-2,
         4.27180926518931511717E-2,
         5.85936634471101055642E-2,
         9.37499997197644278445E-2,
         2.49999999999888314361E-1];
    E(idx) = polyval(P,x(idx)) - log(x(idx)) .* (x(idx) .* polyval(Q,x(idx)) );
  endif
  E = reshape(E,nr,nc);
endfunction



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------