Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : ELLIE.C

/* ellie.c
*
* Incomplete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* double phi, m, y, ellie();
*
* y = ellie( 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
* DEC 0,2 2000 1.9e-16 3.4e-17
* IEEE 0,2 10000 2.2e-15 2.1e-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 second kind */

extern double PI, PIO2, MACHEP;

double ellie( phi, m )
double phi, m;
{
double a, b, c, e, temp;
double lphi, t, step;
double sqrt(), fabs(), log(), sin(), tan(), atan();
double ellpe(), ellpk();
int d, mod, sign;

if( m == 0.0 )
return( phi );
if( m == 1.0 )
return( sin(phi) );
lphi = phi;
if( lphi < 0.0 )
lphi = -lphi;
a = 1.0;
b = 1.0 - m;
b = sqrt(b);
c = sqrt(m);
d = 1;
e = 0.0;
t = tan( lphi );
mod = (lphi + PIO2)/PI;

while( fabs(c/a) > MACHEP )
{
temp = b/a;
lphi = lphi + atan(t*temp) + mod * PI;
mod = (lphi + PIO2)/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;
e += c * sin(lphi);
}

b = 1.0 - m;
temp = ellpe(b)/ellpk(b);
temp *= (atan(t) + mod * PI)/(d * a);
temp += e;
if( phi < 0.0 )
temp = -temp;
return( temp );
}