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
|
/* ellik.c
*
* Incomplete elliptic integral of the first kind
*
*
*
* SYNOPSIS:
*
* double phi, m, y, ellik();
*
* y = ellik( phi, m );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
*
* phi
* -
* | |
* | dt
* F(phi | m) = | ------------------
* | 2
* | | sqrt( 1 - m sin t )
* -
* 0
*
* of amplitude phi and modulus m, using the arithmetic -
* geometric mean algorithm.
*
*
*
*
* ACCURACY:
*
* Tested at random points with m in [0, 1] and phi as indicated.
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,10 200000 7.4e-16 1.0e-16
*
*
*/
/*
* Cephes Math Library Release 2.0: April, 1987
* Copyright 1984, 1987 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Incomplete elliptic integral of first kind */
#include "mconf.h"
extern double MACHEP;
double ellik(phi, m)
double phi, m;
{
double a, b, c, e, temp, t, K;
int d, mod, sign, npio2;
if (m == 0.0)
return (phi);
a = 1.0 - m;
if (a == 0.0) {
if (fabs(phi) >= NPY_PI_2) {
mtherr("ellik", SING);
return (NPY_INFINITY);
}
return (log(tan((NPY_PI_2 + phi) / 2.0)));
}
npio2 = floor(phi / NPY_PI_2);
if (npio2 & 1)
npio2 += 1;
if (npio2) {
K = ellpk(a);
phi = phi - npio2 * NPY_PI_2;
}
else
K = 0.0;
if (phi < 0.0) {
phi = -phi;
sign = -1;
}
else
sign = 0;
b = sqrt(a);
t = tan(phi);
if (fabs(t) > 10.0) {
/* Transform the amplitude */
e = 1.0 / (b * t);
/* ... but avoid multiple recursions. */
if (fabs(e) < 10.0) {
e = atan(e);
if (npio2 == 0)
K = ellpk(a);
temp = K - ellik(e, m);
goto done;
}
}
a = 1.0;
c = sqrt(m);
d = 1;
mod = 0;
while (fabs(c / a) > MACHEP) {
temp = b / a;
phi = phi + atan(t * temp) + mod * NPY_PI;
mod = (phi + NPY_PI_2) / NPY_PI;
t = t * (1.0 + temp) / (1.0 - temp * t * t);
c = (a - b) / 2.0;
temp = sqrt(a * b);
a = (a + b) / 2.0;
b = temp;
d += d;
}
temp = (atan(t) + mod * NPY_PI) / (d * a);
done:
if (sign < 0)
temp = -temp;
temp += npio2 * K;
return (temp);
}
|
