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
|
/* ellikf.c
*
* Incomplete elliptic integral of the first kind
*
*
*
* SYNOPSIS:
*
* float phi, m, y, ellikf();
*
* y = ellikf( 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 phi in [0, 2] and m in
* [0, 1].
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,2 10000 2.9e-7 5.8e-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 first 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);
#else
float sqrtf(), logf(), sinf(), tanf(), atanf();
#endif
float ellikf( float phia, float ma )
{
float phi, m, a, b, c, temp;
float t;
int d, mod, sign;
phi = phia;
m = ma;
if( m == 0.0 )
return( phi );
if( phi < 0.0 )
{
phi = -phi;
sign = -1;
}
else
sign = 0;
a = 1.0;
b = 1.0 - m;
if( b == 0.0 )
return( logf( tanf( 0.5*(PIO2F + phi) ) ) );
b = sqrtf(b);
c = sqrtf(m);
d = 1;
t = tanf( phi );
mod = (phi + PIO2F)/PIF;
while( fabsf(c/a) > MACHEPF )
{
temp = b/a;
phi = phi + atanf(t*temp) + mod * PIF;
mod = (phi + PIO2F)/PIF;
t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
c = ( a - b )/2.0;
temp = sqrtf( a * b );
a = ( a + b )/2.0;
b = temp;
d += d;
}
temp = (atanf(t) + mod * PIF)/(d * a);
if( sign < 0 )
temp = -temp;
return( temp );
}
|
