elliptic.c File Reference

Functions for evaluating complete elliptic integrals of the first and second kind (K and E), as well as their total derivatives with respect to the modulus k. More...

#include "mdb.h"

Go to the source code of this file.

Functions
void	setCeiAccuracy (double newAccuracy)
double	K_cei (double k)
	Compute the complete elliptic integral of the first kind, K(k).
double	E_cei (double k)
	Compute the complete elliptic integral of the second kind, E(k).
double *	KE_cei (double k, double *buffer)
double	dK_cei (double k)
	Compute the total derivative dK/dk of the complete elliptic integral of the first kind.
double	dE_cei (double k)
	Compute the total derivative dE/dk of the complete elliptic integral of the second kind.

Variables
static double	ceiAccuracy = 1e-14

Detailed Description

Functions for evaluating complete elliptic integrals of the first and second kind (K and E), as well as their total derivatives with respect to the modulus k.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, R. Soliday

Definition in file elliptic.c.

Function Documentation

dE_cei()

double dE_cei

(

double

k

)

Compute the total derivative dE/dk of the complete elliptic integral of the second kind.

Uses K(k) and E(k) to determine the derivative with respect to k.

Parameters

[in]

k

The modulus of the elliptic integral.

Returns

The value of dE/dk.

Definition at line 143 of file elliptic.c.

{
  double buffer[2];
  KE_cei(k, buffer);
  return (buffer[1] - buffer[0]) / k;
}

dK_cei()

double dK_cei

(

double

k

)

Compute the total derivative dK/dk of the complete elliptic integral of the first kind.

Uses K(k) and E(k) to determine the derivative with respect to k.

Parameters

[in]

k

The modulus of the elliptic integral.

Returns

The value of dK/dk.

Definition at line 129 of file elliptic.c.

                        {
  double buffer[2];
  KE_cei(k, buffer);
  return (buffer[1] / (1 - k * k) - buffer[0]) / k;
}

E_cei()

double E_cei

(

double

k

)

Compute the complete elliptic integral of the second kind, E(k).

E(k) = ∫_0^(π/2) √(1 - k² sin² θ) dθ

Parameters

[in]

k

The modulus of the elliptic integral.

Returns

The value of E(k).

Definition at line 59 of file elliptic.c.

                       {
  double a0, b0, c0, a1, b1, c1, K, sum, powerOf2;
 
  a0 = 1;
  b0 = sqrt(1 - sqr(k));
  c0 = k;
  sum = sqr(c0);
  powerOf2 = 1;
 
  do {
    /* do two steps of recurrence per pass in the loop */
    a1 = (a0 + b0) / 2;
    b1 = sqrt(a0 * b0);
    c1 = (a0 - b0) / 2;
    sum += sqr(c1) * (powerOf2 *= 2);
    ;
 
    a0 = (a1 + b1) / 2;
    b0 = sqrt(a1 * b1);
    c0 = (a1 - b1) / 2;
    sum += sqr(c0) * (powerOf2 *= 2);
  } while (fabs(c0) > ceiAccuracy);
 
  K = PI / (2 * a0);
  return K * (1 - sum / 2);
}

K_cei()

double K_cei

(

double

k

)

Compute the complete elliptic integral of the first kind, K(k).

K(k) = ∫_0^(π/2) dθ / √(1 - k² sin² θ)

Parameters

[in]

k

The modulus of the elliptic integral.

Returns

The value of K(k).

Definition at line 33 of file elliptic.c.

                       {
  double a0, b0, c0, a1, b1;
 
  a0 = 1;
  b0 = sqrt(1 - sqr(k));
  c0 = k;
 
  do {
    /* do two steps of recurrence per pass in the loop */
    a1 = (a0 + b0) / 2;
    b1 = sqrt(a0 * b0);
    a0 = (a1 + b1) / 2;
    b0 = sqrt(a1 * b1);
    c0 = (a1 - b1) / 2;
  } while (fabs(c0) > ceiAccuracy);
  return PI / (2 * a0);
}

KE_cei()

double * KE_cei	(	double	k,
		double *	buffer )

Definition at line 86 of file elliptic.c.

                                         {
  double a0, b0, c0, a1, b1, c1, K, sum, powerOf2;
 
  if (!buffer)
    buffer = tmalloc(sizeof(*buffer) * 2);
 
  a0 = 1;
  b0 = sqrt(1 - sqr(k));
  c0 = k;
  sum = sqr(c0);
  powerOf2 = 1;
 
  do {
    /* do two steps of recurrence per pass in the loop */
    a1 = (a0 + b0) / 2;
    b1 = sqrt(a0 * b0);
    c1 = (a0 - b0) / 2;
    sum += sqr(c1) * (powerOf2 *= 2);
    ;
 
    a0 = (a1 + b1) / 2;
    b0 = sqrt(a1 * b1);
    c0 = (a1 - b1) / 2;
    sum += sqr(c0) * (powerOf2 *= 2);
  } while (fabs(c0) > ceiAccuracy);
 
  buffer[0] = K = PI / (2 * a0);
  buffer[1] = K * (1 - sum / 2);
  return buffer;
}

setCeiAccuracy()

void setCeiAccuracy

(

double

newAccuracy

)

Definition at line 21 of file elliptic.c.

                                        {
  ceiAccuracy = newAccuracy;
}

Variable Documentation

ceiAccuracy

double ceiAccuracy = 1e-14

static

Definition at line 19 of file elliptic.c.