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
|
/* ellief.c
*
* Incomplete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* float phi, m, y, ellief();
*
* y = ellief( phi, m );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
* phi
* -
* | |
* | 2
* E(phi\m) = | sqrt( 1 - m sin t ) dt
* |
* | |
* -
* 0
*
* of amplitude phi and modulus m, using the arithmetic -
* geometric mean algorithm.
*
*
*
* ACCURACY:
*
* Tested at random arguments with phi in [0, 2] and m in
* [0, 1].
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,2 10000 4.5e-7 7.4e-8
*
*
*/
/*
Cephes Math Library Release 2.2: July, 1992
Copyright 1984, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Incomplete elliptic integral of second kind */
#include <math.h>
extern float PIF, PIO2F, MACHEPF;
#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
#ifdef ANSIC
float sqrtf(float), logf(float), sinf(float), tanf(float), atanf(float);
float ellpef(float), ellpkf(float);
#else
float sqrtf(), logf(), sinf(), tanf(), atanf();
float ellpef(), ellpkf();
#endif
float ellief( float phia, float ma )
{
float phi, m, a, b, c, e, temp;
float lphi, t;
int d, mod;
phi = phia;
m = ma;
if( m == 0.0 )
return( phi );
if( m == 1.0 )
return( sinf(phi) );
lphi = phi;
if( lphi < 0.0 )
lphi = -lphi;
a = 1.0;
b = 1.0 - m;
b = sqrtf(b);
c = sqrtf(m);
d = 1;
e = 0.0;
t = tanf( lphi );
mod = (lphi + PIO2F)/PIF;
while( fabsf(c/a) > MACHEPF )
{
temp = b/a;
lphi = lphi + atanf(t*temp) + mod * PIF;
mod = (lphi + PIO2F)/PIF;
t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
c = 0.5 * ( a - b );
temp = sqrtf( a * b );
a = 0.5 * ( a + b );
b = temp;
d += d;
e += c * sinf(lphi);
}
b = 1.0 - m;
temp = ellpef(b)/ellpkf(b);
temp *= (atanf(t) + mod * PIF)/(d * a);
temp += e;
if( phi < 0.0 )
temp = -temp;
return( temp );
}
|
