1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
|
/* 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(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
* 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
*/
#include "mconf.h"
extern double 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 || cephes_isnan(m)) {
mtherr("ellpj", DOMAIN);
*sn = NPY_NAN;
*cn = NPY_NAN;
*ph = NPY_NAN;
*dn = NPY_NAN;
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)) - NPY_PI_2 + 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);
}
|