/* ellpj.c * * Jacobian Elliptic Functions * * * * SYNOPSIS: * * double u, m, sn, cn, dn, phi; * int 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.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier */ #include <math.h> #ifdef ANSIPROT extern double sqrt ( double ); extern double fabs ( double ); extern double sin ( double ); extern double cos ( double ); extern double asin ( double ); extern double tanh ( double ); extern double sinh ( double ); extern double cosh ( double ); extern double atan ( double ); extern double exp ( double ); #else double sqrt(), fabs(), sin(), cos(), asin(), tanh(); double sinh(), cosh(), atan(), exp(); #endif extern double PIO2, MACHEP; int ellpj( u, m, sn, cn, dn, ph ) double u, m; double *sn, *cn, *dn, *ph; { double ai, b, phi, t, twon; double a[9], c[9]; int i; /* Check for special cases */ if( m < 0.0 || m > 1.0 ) { mtherr( "ellpj", DOMAIN ); *sn = 0.0; *cn = 0.0; *ph = 0.0; *dn = 0.0; return(-1); } 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(0); } 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(0); } /* 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 ) { mtherr( "ellpj", OVERFLOW ); 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(0); }