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
|
/* ellpe.c
*
* Complete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* double m, y, ellpe();
*
* y = ellpe( m );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
* pi/2
* -
* | | 2
* E(m) = | sqrt( 1 - m sin t ) dt
* | |
* -
* 0
*
* Where m = 1 - m1, using the approximation
*
* P(x) - x log x Q(x).
*
* Though there are no singularities, the argument m1 is used
* internally rather than m for compatibility with ellpk().
*
* E(1) = 1; E(0) = pi/2.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 1 10000 2.1e-16 7.3e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* ellpe domain x<0, x>1 0.0
*
*/
/* ellpe.c */
/* Elliptic integral of second kind */
/*
* Cephes Math Library, Release 2.1: February, 1989
* Copyright 1984, 1987, 1989 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*
* Feb, 2002: altered by Travis Oliphant
* so that it is called with argument m
* (which gets immediately converted to m1 = 1-m)
*/
#include "mconf.h"
static double P[] = {
1.53552577301013293365E-4,
2.50888492163602060990E-3,
8.68786816565889628429E-3,
1.07350949056076193403E-2,
7.77395492516787092951E-3,
7.58395289413514708519E-3,
1.15688436810574127319E-2,
2.18317996015557253103E-2,
5.68051945617860553470E-2,
4.43147180560990850618E-1,
1.00000000000000000299E0
};
static double Q[] = {
3.27954898576485872656E-5,
1.00962792679356715133E-3,
6.50609489976927491433E-3,
1.68862163993311317300E-2,
2.61769742454493659583E-2,
3.34833904888224918614E-2,
4.27180926518931511717E-2,
5.85936634471101055642E-2,
9.37499997197644278445E-2,
2.49999999999888314361E-1
};
double ellpe(x)
double x;
{
x = 1.0 - x;
if (x <= 0.0) {
if (x == 0.0)
return (1.0);
mtherr("ellpe", DOMAIN);
return (NPY_NAN);
}
if (x > 1.0) {
return ellpe(1.0 - 1/x) * sqrt(x);
}
return (polevl(x, P, 10) - log(x) * (x * polevl(x, Q, 9)));
}
|