i/elliptic

yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


/*
 * elliptic.i -- $Id$
 * elliptic functions
 *
 * Copyright (c) 1999.  See accompanying LEGAL file for details.
 */

/* Abramowitz and Stegun, sections 16.4, 17.2.17-19, 17.6 */


local elliptic ;

/* DOCUMENT elliptic, ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e


     The elliptic integral of the first kind is:


        u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )


     The functions ell_f and ell_am compute this integral and its
     inverse:

        u   = ell_f(phi, m)
	phi = ell_am(u, m)


     The Jacobian elliptic functions can be computed from the

     "amplitude" ell_am by means of:


        sn(u|m) = sin(ell_am(u,m))

	cn(u|m) = cos(ell_am(u,m))

	dn(u|m) = dn_(ell_am(u,m)) = sqrt(1-m*sn(u|m)^2)

     The other nine functions are sc=sn/cn, cs=cn/sn, nd=1/dn,

     cd=cn/dn, dc=dn/cn, ns=1/sn, sd=sn/dn, nc=1/cn, and ds=dn/sn.
     (The notation u|m does not means yorick's | operator; it is
     the mathematical notation, not valid yorick code!)

     The parameter M is given in three different notations:
       as M, the "parameter",
       as k, the "modulus", or
       as alpha, the "modular angle",

     which are related by: M = k^2 = sin(alpha)^2.  The yorick elliptic
     functions in terms of M may need to be written

       ell_am(u,k^2) or ell_am(u,sin(alpha)^2)
     in order to agree with the definitions in other references.
     Sections 17.2.17-19 of Abramowitz and Stegun explains these notations,
     and chapters 16 and 17 present a compact overview of the subject of

     elliptic functions in general.


     The parameter M must be a scalar; U may be an array.  The

     exceptions are the complete elliptic integrals ellip_k and

     ellip_e which accept an array of M values.


     The ell_am function uses the external variable ell_m if M is
     omitted, otherwise stores M in ell_m.  Hence, you may set ell_m,

     then simply call ell_am(u) if you have a series of calls with

     the same value of M; this also allows the dn_ function to work

     without a second specification of M.


     The elliptic integral of the second kind is:


        u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) )


     The function ell_e computes this integral:

        u   = ell_e(phi, m)


     The special values ell_f(pi/2,m) and ell_e(pi/2,m) are the complete

     elliptic integrals of the first and second kinds; separate functions

     ellip_k and ellip_e are provided to compute them.


     Note that the function ellip_k is infinite for M=1 and for large

     negative M.  The "natural" range for M is 0<=M<=1; all other real

     values can be "reduced" to this range by various transformations;

     the logarithmic singularity of ellip_k is actually very mild, and

     other functions such as ell_am are perfectly well-defined there.


     Here are the sum formulas for elliptic functions:

       sn(u+v) = ( sn(u)*cn(v)*dn(v) + sn(v)*cn(u)*dn(u) ) /
                 ( 1 - m*sn(u)^2*sn(v)^2 )
       cn(u+v) = ( cn(u)*cn(v) - sn(u)*dn(u)*sn(v)*dn(v) ) /
                 ( 1 - m*sn(u)^2*sn(v)^2 )
       dn(u+v) = ( dn(u)*dn(v) - m*sn(u)*cn(u)*sn(v)*cn(v) ) /
                 ( 1 - m*sn(u)^2*sn(v)^2 )

     And the formulas for pure imaginary values:

       sn(1i*u,m) = 1i * sc(u,1-m)
       cn(1i*u,m) = nc(u,1-m)
       dn(1i*u,m) = dc(u,1-m)


   SEE ALSO: ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e
 */


func ell_am (u,m)
/* DOCUMENT ell_am(u)
         or ell_am(u,m)

     returns the "amplitude" (an angle in radians) for the Jacobi

     elliptic functions at U, with parameter M.  That is,
        phi = ell_am(u,m)
     means that

        u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )


     Thus ell_am is the inverse of the incomplete elliptic function

     of the first kind ell_f.  See help,elliptic for more.

   SEE ALSO: elliptic
 */
{
  /* set up the arithmetic-geometric mean scale */
  extern ell_m, _agm_m, _agm_n, _agm_coa, _agm_a, _agm_sn;
  if (is_void(m)) m= ell_m;
  if (m != ell_m) {
    ell_m= m= double(m);
    if (m<0.) {
      _agm_sn= 1./(1.-m);
      m*= -_agm_sn;
      _agm_sn= sqrt(_agm_sn);
    } else if (m>1.) {
      m= 1./m;
      _agm_sn= sqrt(m);
    }
    _agm_m= m;
    _agm_coa()= 0.;
    _agm_n= 0;
    _agm_a= 1.;
    if (m!=1.) {
      b= sqrt(1.-m);
      for (;;) {         /* maximum of 8 passes for 64-bit double */
	c= 0.5*(_agm_a-b);
	if (!c) break;
	am= _agm_a-c;         /* arithmetic mean */
	_agm_coa(++_agm_n)= c/am;
	if (am==_agm_a) break;

	b= sqrt(_agm_a*b);    /* geometric mean */
	_agm_a= am;
      }
      _agm_a*= 2.^_agm_n;
    }
  }

  phi= _agm_a*u;
  if (m<0. || m>1.) phi/= _agm_sn;
  if (m!=1.)

    for (n=_agm_n ; n>0 ; n--) phi= 0.5*(phi+asin(_agm_coa(n)*sin(phi)));
  else

    phi= atan(tanh(phi), sech(phi));
  if (m<0.) {
    cn= cos(phi);
    phi= sin(phi);
    nd= 1./sqrt(1.-m*phi*phi);
    phi= atan(_agm_sn*phi*nd, cn*nd);
  } else if (m>1.) {
    phi= sin(phi);

    phi= atan(_agm_sn*phi, sqrt(1.-_agm_m*phi*phi));
  }
  return phi;
}

_agm_coa= array(0.,16);
_agm_n= 0;
_agm_m= 0.0;
_agm_a= 1.0;
_agm_sn= 1.0;
ell_m= 0.0;


func dn_ (phi, m)

/* DOCUMENT dn_(ell_am(u,m))


     return the Jacobian elliptic function dn(u|m).  The external
     variable ell_m must be set properly before calling dn_.


   SEE ALSO: elliptic, ell_am
 */
{
  if (is_void(m)) m= ell_m;
  phi= sin(phi);
  return sqrt(1.-m*phi*phi);
}


func ell_f (phi,m)
/* DOCUMENT ell_f(phi,m)


     returns the incomplete elliptic integral of the first kind F(phi|M).
     That is,
        u = ell_f(phi,m)
     means that

        u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )


     See help,elliptic for more.


   SEE ALSO: elliptic, ell_e
 */
{
  orig_m= m= double(m);
  if (m>1.) {
    scale= sqrt(m);

    phi= asin(scale*sin(phi));
    m= 1./m;
  } else if (m<0.) {
    scale= sqrt(1.-m);
    phi= 0.5*pi - phi;
    m/= m-1.;
  }

  if (m==1.) {

    phi= atanh(sin(phi));
    a= 0.;
  } else {          /* compute using arithmetic-geometric mean */
    sgn= sign(phi);
    phi= abs(phi);
    a= 1.;
    b= sqrt(1.-m);
    twon= 1.0;
    pi2= 0.5*pi;
    for (;;) {      /* maximum of 8 passes for 64-bit double */
      c= 0.5*(a-b);
      if (!c) break;
      phase= (phi+pi2)/pi;
      cycle= floor(phase);
      phi*= 1. + 1.e-15*(cycle==phase);

      phi+= atan((b/a)*tan(phi)) + pi*cycle;
      twon*= 2.0;
      am= a-c;
      if (am==a) break;
      b= sqrt(a*b);
      a= am;
    }
    phi/= twon*a*sgn;
  }

  if (orig_m>1.) {
    phi/= scale;
  } else if (orig_m<0.) {
    phi= (0.5*pi/a - phi)/scale;
  }

  return phi;
}


func ell_e (phi,m)
/* DOCUMENT ell_e(phi,m)


     returns the incomplete elliptic integral of the second kind E(phi|M).
     That is,
        u = ell_e(phi,m)
     means that

        u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) )


     See help,elliptic for more.


   SEE ALSO: elliptic, ell_f
 */
{
  orig_m= m= double(m);
  if (m>1.) {
    scale= sqrt(m);

    phi= asin(scale*sin(phi));
    m= 1./m;
  } else if (m<0.) {
    scale= sqrt(1.-m);
    phi= 0.5*pi - phi;
    m/= m-1.;
  }

  if (m==1.) {

    per= floor((phi+0.5*pi)/pi);

    phi= sin(phi)*sign(0.5-abs(per%2.)) + 2.*per;
    a= 0.;
  } else {          /* compute using arithmetic-geometric mean */
    sgn= sign(phi);
    phi= abs(phi);
    a= 1.;
    b= sqrt(1.-m);
    eok= 1.-0.5*m;
    cs= 0.;
    twon= 1.0;
    pi2= 0.5*pi;
    for (;;) {      /* maximum of 8 passes for 64-bit double */
      c= 0.5*(a-b);
      if (!c) break;

      phi+= atan((b/a)*tan(phi)) + pi*floor((phi+pi2)/pi);
      cs+= c*sin(phi);
      eok-= twon*c*c;
      twon*= 2.0;
      am= a-c;
      if (am==a) break;
      b= sqrt(a*b);
      a= am;
    }
    f= sgn*phi/(twon*a);
    phi= eok*f + sgn*cs;
  }

  if (orig_m>1.) {
    phi= (phi-(1.-m)*f)/scale;
  } else if (orig_m<0.) {
    phi= (eok*0.5*pi/a - phi)*scale;
  }

  return phi;
}


func ellip_k (m)
/* DOCUMENT ellip_k(m)


     returns the complete elliptic integral of the first kind K(M):

        K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) )


     See help,elliptic for more.


   SEE ALSO: elliptic, ellip_e, ell_f
 */
{

  if (anyof(m>=1.)) error, "ellip_k(m) not computed for m>=1";

  m= double(m);
  mask= (m>=0.);
  list= where(mask);

  if (numberof(list)) scale= array(0.5*pi,numberof(list));
  list= where(!mask);
  if (numberof(list)) {
    sm= 1./(1.-m(list));
    m(list)*= -sm;

    sm= 0.5*pi*sqrt(sm);
  }
  scale= merge(scale,sm,mask);


  a= array(1.,dimsof(m));
  b= sqrt(1.-m);
  for (;;) {
    c= 0.5*(a-b);
    am= a-c;
    if (allof(am==a)) break;
    b= sqrt(a*b);
    a= am;
  }

  return scale/a;
}


func ellip_e (m)
/* DOCUMENT ellip_e(m)


     returns the complete elliptic integral of the second kind E(M):

        E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) )


     See help,elliptic for more.


   SEE ALSO: elliptic, ellip_k, ell_e
 */
{

  if (anyof(m>1.)) error, "ellip_e(m) not computed for m>1";

  m= double(m);
  mask= (m>=0.);
  list= where(mask);

  if (numberof(list)) scale= array(0.5*pi,numberof(list));
  list= where(!mask);
  if (numberof(list)) {
    sm= 1.-m(list);
    m(list)/= -sm;

    sm= 0.5*pi*sqrt(sm);
  }
  scale= merge(scale,sm,mask);

  mask= (m!=1.);
  list= where(mask);
  if (numberof(list)) {
    m= m(list);

    a= array(1.,numberof(m));
    b= sqrt(1.-m);
    e= 1.-0.5*m;
    twon= 1.0;
    for (n=0 ;;) {  /* maximum of 8 passes for 64-bit double */
      c= 0.5*(a-b);
      am= a-c;
      if (allof(am==a)) break;
      e-= twon*c*c;
      twon*= 2.0;
      b= sqrt(a*b);
      a= am;
    }
    e/= a;
  }
  list= where(!mask);

  if (numberof(list)) em= array(2./pi,numberof(list));

  return scale*merge(e,em,mask);
}


func ellip2_k (m)
/* DOCUMENT ellip2_k(m)


     returns the complete elliptic integral of the first kind K(M):

        K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) )
     accurate to 2e-8 for 0<=M<1


   SEE ALSO: elliptic, ellip_e, ell_f
 */
{
  m= 1.-m;

  return poly(m,1.38629436112,0.09666344259,0.03590092383,
	      0.03742563713,0.01451196212) + 

         log(1./m)*poly(m,0.5,0.12498593597,0.06880248576,
			0.03328355346,0.00441787012);
}


func ellip2_e (m)
/* DOCUMENT ellip2_e(m)


     returns the complete elliptic integral of the second kind E(M):

        E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) )
     accurate to 2e-8 for 0<=M<=1


   SEE ALSO: elliptic, ellip_k, ell_e
 */
{
  m= 1.-m;

  return poly(m,1.,0.44325141463,0.06260601220,
	      0.04757383546,0.01736506451) +

         log(1./m)*poly(m,0.,0.24998368310,0.09200180037,
			0.04069697526,0.00526449639);
}