lib/libm/ellie.c

   1 /*                                                      ellie.c
   2  *
   3  *      Incomplete elliptic integral of the second kind
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double phi, m, y, ellie();
  10  *
  11  * y = ellie( phi, m );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Approximates the integral
  18  *
  19  *
  20  *                phi
  21  *                 -
  22  *                | |
  23  *                |                   2
  24  * E(phi_\m)  =    |    sqrt( 1 - m sin t ) dt
  25  *                |
  26  *              | |
  27  *               -
  28  *                0
  29  *
  30  * of amplitude phi and modulus m, using the arithmetic -
  31  * geometric mean algorithm.
  32  *
  33  *
  34  *
  35  * ACCURACY:
  36  *
  37  * Tested at random arguments with phi in [-10, 10] and m in
  38  * [0, 1].
  39  *                      Relative error:
  40  * arithmetic   domain     # trials      peak         rms
  41  *    DEC        0,2         2000       1.9e-16     3.4e-17
  42  *    IEEE     -10,10      150000       3.3e-15     1.4e-16
  43  *
  44  *
  45  */
  46 \f
  47
  48 /*
  49 Cephes Math Library Release 2.8:  June, 2000
  50 Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
  51 */
  52
  53 /*      Incomplete elliptic integral of second kind     */
  54 #include "mconf.h"
  55 extern double PI, PIO2, MACHEP;
  56 #ifdef ANSIPROT
  57 extern double sqrt ( double );
  58 extern double fabs ( double );
  59 extern double log ( double );
  60 extern double sin ( double x );
  61 extern double tan ( double x );
  62 extern double atan ( double );
  63 extern double floor ( double );
  64 extern double ellpe ( double );
  65 extern double ellpk ( double );
  66 double ellie ( double, double );
  67 #else
  68 double sqrt(), fabs(), log(), sin(), tan(), atan(), floor();
  69 double ellpe(), ellpk(), ellie();
  70 #endif
  71
  72 double ellie( phi, m )
  73 double phi, m;
  74 {
  75 double a, b, c, e, temp;
  76 double lphi, t, E;
  77 int d, mod, npio2, sign;
  78
  79 if( m == 0.0 )
  80         return( phi );
  81 lphi = phi;
  82 npio2 = floor( lphi/PIO2 );
  83 if( npio2 & 1 )
  84         npio2 += 1;
  85 lphi = lphi - npio2 * PIO2;
  86 if( lphi < 0.0 )
  87         {
  88         lphi = -lphi;
  89         sign = -1;
  90         }
  91 else
  92         {
  93         sign = 1;
  94         }
  95 a = 1.0 - m;
  96 E = ellpe( a );
  97 if( a == 0.0 )
  98         {
  99         temp = sin( lphi );
 100         goto done;
 101         }
 102 t = tan( lphi );
 103 b = sqrt(a);
 104 /* Thanks to Brian Fitzgerald <fitzgb@mml0.meche.rpi.edu>
 105    for pointing out an instability near odd multiples of pi/2.  */
 106 if( fabs(t) > 10.0 )
 107         {
 108         /* Transform the amplitude */
 109         e = 1.0/(b*t);
 110         /* ... but avoid multiple recursions.  */
 111         if( fabs(e) < 10.0 )
 112                 {
 113                 e = atan(e);
 114                 temp = E + m * sin( lphi ) * sin( e ) - ellie( e, m );
 115                 goto done;
 116                 }
 117         }
 118 c = sqrt(m);
 119 a = 1.0;
 120 d = 1;
 121 e = 0.0;
 122 mod = 0;
 123
 124 while( fabs(c/a) > MACHEP )
 125         {
 126         temp = b/a;
 127         lphi = lphi + atan(t*temp) + mod * PI;
 128         mod = (lphi + PIO2)/PI;
 129         t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
 130         c = ( a - b )/2.0;
 131         temp = sqrt( a * b );
 132         a = ( a + b )/2.0;
 133         b = temp;
 134         d += d;
 135         e += c * sin(lphi);
 136         }
 137
 138 temp = E / ellpk( 1.0 - m );
 139 temp *= (atan(t) + mod * PI)/(d * a);
 140 temp += e;
 141
 142 done:
 143
 144 if( sign < 0 )
 145         temp = -temp;
 146 temp += npio2 * E;
 147 return( temp );
 148 }