# Category : C Source Code

Archive : CEPHES22.ZIP

Filename : ELLPJ.C

*

* Jacobian Elliptic Functions

*

*

*

* SYNOPSIS:

*

* double u, m, sn, cn, dn, phi, ellpj();

*

* ellpj( u, m, _&sn, _&cn, _&dn, _&phi );

*

*

*

* DESCRIPTION:

*

*

* 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.

*

* 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 = ellik(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.

*

* ACCURACY:

*

* Tested at random points with u between 0 and 10, m between

* 0 and 1.

*

* Absolute error (* = relative error):

* arithmetic function # trials peak rms

* DEC sn 1800 4.5e-16 8.7e-17

* IEEE phi 10000 9.2e-16* 1.4e-16*

* IEEE sn 50000 4.1e-15 4.6e-16

* IEEE cn 40000 3.6e-15 4.4e-16

* IEEE dn 10000 1.3e-12 1.8e-14

*

* Peak error observed in consistency check using addition

* theorem for sn(u+v) was 4e-16 (absolute). Also tested by

* the above relation to the incomplete elliptic integral.

* Accuracy deteriorates when u is large.

*

*/

/* ellpj.c */

/*

Cephes Math Library Release 2.0: April, 1987

Copyright 1984, 1987 by Stephen L. Moshier

Direct inquiries to 30 Frost Street, Cambridge, MA 02140

*/

extern double PIO2, MACHEP;

double ellpj( u, m, sn, cn, dn, ph )

double u, m;

double *sn, *cn, *dn, *ph;

{

double ai, b, phi, t, twon;

double sqrt(), fabs(), sin(), cos(), asin(), tanh();

double sinh(), cosh(), atan(), exp();

double a[9], c[9];

int i;

/* Check for special cases */

if( m < 0.0 || m > 1.0 )

return( puts( "ellpj m out of range" ) );

if( m < 1.0e-9 )

{

t = sin(u);

b = cos(u);

ai = 0.25 * m * (u - t*b);

*sn = t - ai*b;

*cn = b + ai*t;

*ph = u - ai;

*dn = 1.0 - 0.5*m*t*t;

return;

}

if( m >= 0.9999999999 )

{

ai = 0.25 * (1.0-m);

b = cosh(u);

t = tanh(u);

phi = 1.0/b;

twon = b * sinh(u);

*sn = t + ai * (twon - u)/(b*b);

*ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;

ai *= t * phi;

*cn = phi - ai * (twon - u);

*dn = phi + ai * (twon + u);

return;

}

/* A. G. M. scale */

a[0] = 1.0;

b = sqrt(1.0 - m);

c[0] = sqrt(m);

twon = 1.0;

i = 0;

while( fabs(c[i]/a[i]) > MACHEP )

{

if( i > 7 )

{

printf( "ellpj() array full" );

goto done;

}

ai = a[i];

++i;

c[i] = ( ai - b )/2.0;

t = sqrt( ai * b );

a[i] = ( ai + b )/2.0;

b = t;

twon *= 2.0;

}

done:

/* backward recurrence */

phi = twon * a[i] * u;

do

{

t = c[i] * sin(phi) / a[i];

b = phi;

phi = (asin(t) + phi)/2.0;

}

while( --i );

*sn = sin(phi);

t = cos(phi);

*cn = t;

*dn = t/cos(phi-b);

*ph = phi;

return;

}

Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/