home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ellf / ellik.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.8 KB  |  105 lines
  1. /*                            ellik.c
  2.  *
  3.  *    Incomplete elliptic integral of the first kind
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double phi, m, y, ellik();
  10.  *
  11.  * y = ellik( phi, m );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Approximates the integral
  18.  *
  19.  *
  20.  *
  21.  *                phi
  22.  *                 -
  23.  *                | |
  24.  *                |           dt
  25.  * F(phi_\m)  =    |    ------------------
  26.  *                |                   2
  27.  *              | |    sqrt( 1 - m sin t )
  28.  *               -
  29.  *                0
  30.  *
  31.  * of amplitude phi and modulus m, using the arithmetic -
  32.  * geometric mean algorithm.
  33.  *
  34.  *
  35.  *
  36.  *
  37.  * ACCURACY:
  38.  *
  39.  * Tested at random points with phi in [0, 2] and m in
  40.  * [0, 1].
  41.  *                      Relative error:
  42.  * arithmetic   domain     # trials      peak         rms
  43.  *    DEC       0,2          3700       8.1e-17     2.5e-17
  44.  *    IEEE      0,2         10000       6.0e-16     1.4e-16
  45.  *
  46.  *
  47.  */
  48.  
  49.  
  50. /*
  51. Cephes Math Library Release 2.0:  April, 1987
  52. Copyright 1984, 1987 by Stephen L. Moshier
  53. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54. */
  55.  
  56. /*    Incomplete elliptic integral of first kind    */
  57.  
  58. extern double PI, PIO2, MACHEP;
  59.  
  60. double ellik( phi, m )
  61. double phi, m;
  62. {
  63. double a, b, c, temp;
  64. double t, step;
  65. double sqrt(), fabs(), log(), tan(), atan();
  66. int d, mod, sign;
  67.  
  68. if( m == 0.0 )
  69.     return( phi );
  70. if( phi < 0.0 )
  71.     {
  72.     phi = -phi;
  73.     sign = -1;
  74.     }
  75. else
  76.     sign = 0;
  77. a = 1.0;
  78. b = 1.0 - m;
  79. if( b == 0.0 )
  80.     return(  log(  tan( (PIO2 + phi)/2.0 )  )   );
  81. b = sqrt(b);
  82. c = sqrt(m);
  83. d = 1;
  84. t = tan( phi );
  85. mod = (phi + PIO2)/PI;
  86.  
  87. while( fabs(c/a) > MACHEP )
  88.     {
  89.     temp = b/a;
  90.     phi = phi + atan(t*temp) + mod * PI;
  91.     mod = (phi + PIO2)/PI;
  92.     t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
  93.     c = ( a - b )/2.0;
  94.     temp = sqrt( a * b );
  95.     a = ( a + b )/2.0;
  96.     b = temp;
  97.     d += d;
  98.     }
  99.  
  100. temp = (atan(t) + mod * PI)/(d * a);
  101. if( sign < 0 )
  102.     temp = -temp;
  103. return( temp );
  104. }
  105.