elliptic.c File Reference

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.

#include "mdb.h"
#include "mdb_thread.h"

Go to the source code of this file.

Functions
static double	getCeiAccuracy (void)
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.

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 158 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 144 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 72 of file elliptic.c.

                       {
  double a0, b0, c0, a1, b1, c1, K, sum, powerOf2;
  double accuracy = getCeiAccuracy();
 
  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) > accuracy);
 
  K = PI / (2 * a0);
  return K * (1 - sum / 2);
}

getCeiAccuracy()

static double getCeiAccuracy ( void )

static

Definition at line 23 of file elliptic.c.

                                   {
  double accuracy;
  mdb_thread_lock(&cei_accuracy_lock);
  accuracy = ceiAccuracy;
  mdb_thread_unlock(&cei_accuracy_lock);
  return accuracy;
}

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 45 of file elliptic.c.

                       {
  double a0, b0, c0, a1, b1;
  double accuracy = getCeiAccuracy();
 
  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) > accuracy);
  return PI / (2 * a0);
}

KE_cei()

double * KE_cei	(	double	k,
		double *	buffer )

Definition at line 100 of file elliptic.c.

                                         {
  double a0, b0, c0, a1, b1, c1, K, sum, powerOf2;
  double accuracy = getCeiAccuracy();
 
  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) > accuracy);
 
  buffer[0] = K = PI / (2 * a0);
  buffer[1] = K * (1 - sum / 2);
  return buffer;
}

setCeiAccuracy()

void setCeiAccuracy

(

double

newAccuracy

)

Definition at line 31 of file elliptic.c.

                                        {
  mdb_thread_lock(&cei_accuracy_lock);
  ceiAccuracy = newAccuracy;
  mdb_thread_unlock(&cei_accuracy_lock);
}