Quelle  ellik.c

/* ellik.c
 *
 * Incomplete elliptic integral of the first kind
 *
 *
 *
 * SYNOPSIS:
 *
 * double phi, m, y, ellik();
 *
 * y = ellik( phi, m );
 *
 *
 *
 * DESCRIPTION:
 *
 * Approximates the integral
 *
 *
 *
 *                phi
 *                 -
 *                | |
 *                |           dt
 * F(phi_\m)  =    |    ------------------
 *                |                   2
 *              | |    sqrt( 1 - m sin t )
 *               -
 *                0
 *
 * of amplitude phi and modulus m, using the arithmetic -
 * geometric mean algorithm.
 *
 *
 *
 *
 * ACCURACY:
 *
 * Tested at random points with m in [0, 1] and phi as indicated.
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE     -10,10       200000      7.4e-16     1.0e-16
 *
 *
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/

/* Incomplete elliptic integral of first kind */

#include "mconf.h"
#ifdef ANSIPROT
extern double sqrt(double);
extern double fabs(double);
extern double log(double);
extern double tan(double);
extern double atan(double);
extern double floor(double);
extern double ellpk(double);
double ellik(double, double);
#else
double sqrt(), fabs(), log(), tan(), atan(), floor(), ellpk();
double ellik();
#endif
extern double PI, PIO2, MACHEP, MAXNUM;

double ellik(phi, m) double phi, m;
{
  double a, b, c, e, temp, t, K;
  int d, mod, sign, npio2;

  if (m == 0.0)
    return (phi);
  a = 1.0 - m;
  if (a == 0.0) {
    if (fabs(phi) >= PIO2) {
      mtherr("ellik", SING);
      return (MAXNUM);
    }
    return (log(tan((PIO2 + phi) / 2.0)));
  }
  npio2 = floor(phi / PIO2);
  if (npio2 & 1)
    npio2 += 1;
  if (npio2) {
    K = ellpk(a);
    phi = phi - npio2 * PIO2;
  } else
    K = 0.0;
  if (phi < 0.0) {
    phi = -phi;
    sign = -1;
  } else
    sign = 0;
  b = sqrt(a);
  t = tan(phi);
  if (fabs(t) > 10.0) {
    /* Transform the amplitude */
    e = 1.0 / (b * t);
    /* ... but avoid multiple recursions.  */
    if (fabs(e) < 10.0) {
      e = atan(e);
      if (npio2 == 0)
        K = ellpk(a);
      temp = K - ellik(e, m);
      goto done;
    }
  }
  a = 1.0;
  c = sqrt(m);
  d = 1;
  mod = 0;

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

  temp = (atan(t) + mod * PI) / (d * a);

done:
  if (sign < 0)
    temp = -temp;
  temp += npio2 * K;
  return (temp);
}