home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ellf / ellie.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.9 KB  |  105 lines
  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 [0, 2] 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       0,2        10000       2.2e-15     2.1e-16
  43.  *
  44.  *
  45.  */
  46.  
  47.  
  48. /*
  49. Cephes Math Library Release 2.0:  April, 1987
  50. Copyright 1984, 1987 by Stephen L. Moshier
  51. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  52. */
  53.  
  54. /*    Incomplete elliptic integral of second kind    */
  55.  
  56. extern double PI, PIO2, MACHEP;
  57.  
  58. double ellie( phi, m )
  59. double phi, m;
  60. {
  61. double a, b, c, e, temp;
  62. double lphi, t, step;
  63. double sqrt(), fabs(), log(), sin(), tan(), atan();
  64. double ellpe(), ellpk();
  65. int d, mod, sign;
  66.  
  67. if( m == 0.0 )
  68.     return( phi );
  69. if( m == 1.0 )
  70.     return( sin(phi) );
  71. lphi = phi;
  72. if( lphi < 0.0 )
  73.     lphi = -lphi;
  74. a = 1.0;
  75. b = 1.0 - m;
  76. b = sqrt(b);
  77. c = sqrt(m);
  78. d = 1;
  79. e = 0.0;
  80. t = tan( lphi );
  81. mod = (lphi + PIO2)/PI;
  82.  
  83. while( fabs(c/a) > MACHEP )
  84.     {
  85.     temp = b/a;
  86.     lphi = lphi + atan(t*temp) + mod * PI;
  87.     mod = (lphi + PIO2)/PI;
  88.     t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
  89.     c = ( a - b )/2.0;
  90.     temp = sqrt( a * b );
  91.     a = ( a + b )/2.0;
  92.     b = temp;
  93.     d += d;
  94.     e += c * sin(lphi);
  95.     }
  96.  
  97. b = 1.0 - m;
  98. temp = ellpe(b)/ellpk(b);
  99. temp *= (atan(t) + mod * PI)/(d * a);
  100. temp += e;
  101. if( phi < 0.0 )
  102.     temp = -temp;
  103. return( temp );
  104. }
  105.