Digital Signal Processing Reference
In-Depth Information
as 100 and 1E-15, respectively, in the ADV_MATH.H include file. Of course
these values could be changed as necessary.
/*====================================================
Ellip_Integral() - calcs complete elliptic integral
using arithmetic-geometric mean method
Prototype: void Ellip_Integral(double k);
Return: complete elliptic integral value
Arguments: k - the modulus of the integral
====================================================*/
double Ellip_Integral(double k)
{ int i; /* Loop counter. */
double A[MAX_TERMS],B[MAX_TERMS],
C[MAX_TERMS]; /* Array storage values. */
/* Square the modulus as required by this method.*/
k = k * k;
/* Initialize the starting values. */
A[0] = 1;
B[0] = sqrt(1-k);
C[0] = sqrt(k);
/* Iterate until error is small enough. */
for(i = 1; i < MAX_TERMS ;i++)
{ A[i] = (A[i-1] + B[i-1])/2;
B[i] = sqrt(A[i-1]*B[i-1]);
C[i] = (A[i-1] - B[i-1])/2;
if(C[i] < ERR_SMALL)
{ break;}
}
return PI / (2 * A[i]);
}
Listing D.5
Ellip_Integral
function.
