gov.nih.mipav.model.algorithms
Class EllipticIntegral

java.lang.Object
  gov.nih.mipav.model.algorithms.EllipticIntegral

public class EllipticIntegral
extends Object

DOCUMENT ME!

Author:

ilb The elliptic integral of the first kind is defined by F(k, phi) = Integral from 0 to phi of d(theta)/sqrt(1 - k*k*sin(theta)*sin(theta)) where k is the modulus Letting t = sin(theta), the elliptic integral of the first kind can be written as F(k, x) = Integral from 0 to x of dt/sqrt((1-t*t)(1-k*k*t*t)) where x = sin(phi)

The elliptic integral of the second kind is defined by E(k, phi) = Integral from 0 to phi of sqrt(1 - k*k*sin(theta)*sin(theta))d(theta) where k is the modulus Letting t = sin(theta), the elliptic integral of the second kind can be written as E(k, x) = Integral from 0 to x of (sqrt(1 - k*k*t*t)/sqrt(1 - t*t))dt where x = sin(phi)

The elliptic integral of the third kind is defined by Pi(k, c, phi) = Integral from 0 to phi of d(theta)/((1 - c*sin(theta)*sin(theta))*sqrt(1 - k*k*sin(theta)*sin(theta))) Letting t = sin(theta), the elliptic integral of the third kind can be written as Pi(k, c, x) = Integral from 0 to x of dt/((1 - c*t*t)*sqrt((1 - t*t)*(1 - k*k*t*t))) where x = sin(phi)

When phi = PI/2 (x = 1), the elliptic integrals defined above are called complete elliptic integrals When phi is not equal to PI/2, the elliptic integrals are often referred to as incomplete elliptic integrals The complete elliptic integral of the first kind is: K(k) = Integral from 0 to PI/2 of d(theta)/sqrt(1 - k*k*sin(theta)*sin(theta)) where k is the modulus = Integral from 0 to 1 of dt/sqrt((1-t*t)(1-k*k*t*t)) The complete elliptic integral of the second kind is: E(k, phi) = Integral from 0 to PI/2 of sqrt(1 - k*k*sin(theta)*sin(theta))d(theta) where k is the modulus = Integral from 0 to 1 of (sqrt(1 - k*k*t*t)/sqrt(1 - t*t))dt The comlplete elliptic integral of the third kind is: Pi(k, c, phi) = Integral from 0 to PI/2 of d(theta)/((1 - c*sin(theta)*sin(theta))*sqrt(1 - k*k*sin(theta)*sin(theta))) = Integral from 0 to 1 of dt/((1 - c*t*t)*sqrt((1 - t*t)*(1 - k*k*t*t)))

The code for EllipticIntegralType = COMPLETE1, COMPLETE2, INCOMPLETE1, and INCOMPLETE2 is ported from a C++ program in Figures 1 and 2 of "Numerical Calculation of the Elliptic Integrals of the First and Second Kinds with Complex Modulus" by Tohru Morita, Interdisciplinary Information Sciences, Vol. 6, No. 1, 2000, pp. 67-74.

How to doublecheck? Use a combination of 2 sources: 1.) Port of FORTRAN code for complete and incomplete elliptic integrals of the first kind with real arguments from Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley & Sons, Inc., 1996, Chapter 18, Elliptic Integrals and Jacobian Elliptic Functions, pp. 654-679. 2.) PORT of FORTRAN code http://wwwasd.web.cern.ch/wwwasd/lhc++/requirements/special_functions/celint64.f revision 1.1.1.1 1996/04/01 by mclareni. This code is a translation of the Algol procedure elco2(x,y,kc,a,b,u,v) in Ronald Bulirsch, Handbook Series Special Functions Numerical Calculation of Elliptic Integrals and Elliptic Functions, Numerische Mathematik, Volume 7, 1965, pp. 78 - 90. elco2 calculates a general elliptic integral of the second kind. It calculates w = u + i*v = Integral from 0 to z = x + i*y of (a + b*e*e)/((1 + e*e)*sqrt((1 + e*e)*(1 + kc*kc*e*e)))de kc is the complementary modulus = sqrt(1 - k*k), where k is the modulus The original algol elco2 function fails if kc = 0 or x < 0. However, the FORTRAN code appears to work for kc = 0. The elco2 may have reduced accuracy at the branchpoints x = 0, y = +-i, +-i/kc.

F(k, arctan(w)) = elco2(w, kc, a = 1, b = 1) Let arctan(w) = z. Then w = tan(z) and for the incomplete elliptic integral of the first kind: F(k, z) = elco2(tan(z), kc, a = 1, b = 1) So elco2 cannot be used to check complex k (modulus), but it can be used as a partial check on complex z. It can check on complex z as long as 2 conditons are met: 1.) Real part of tan(z) >= 0 and 2.) kc > 0

E(k, arctan(w)) = elco2(w, kc, a = 1, b = kc*kc) Let arctan(w) = z. Then w = tan(z) and for the incomplete elliptic integral of the second kind: E(k, z) = elco2(tan(z), kc, a = 1, b = kc*kc) and it can check the incomplete elliptic integral of the second kind for the same conditions that the incomplete elliptic integral of the first kind can be checked.

Test results: For complete elliptic integrals of the first and second kind with real arguments, the Morita arithmetic-geometric mean, Morita Carlson's, ZHANG comelp, and ZHANG elit all give the same answer.

For incomplete elliptic integrals of the first and second kind with real arguemnts, both Morita Carlson's, ZHANG elit and CERN elco2 give the same answer. The problem with the CERN code is that is fails if the real part of z = x + iy is less than zero. However, from the Wolfram EllipticF documentation, the F(z|m) is an odd function with respect to z. F(-z|m) = -F(z|m), so there is no need to use elco2 for values of x < 0. Note that the Wolfram m = k*k, where k is the modulus used in Computation of Special Functions.

For incomplete elliptic integrals of the first and second kind with complex phi and real modulus, the Morita standard Carlson's and the CERN elco2 give the same answer.

Constructor Summary

EllipticIntegral()
~ Constructors ---------------------------------------------------------------------------------------------------

EllipticIntegral(boolean complete, double modr, double modi, boolean complementaryModulusUsed, boolean analyticContinuationUsed, double phir, double phii, double[] firstr, double[] firsti, double[] secondr, double[] secondi, boolean useStandardMethod, int[] errorFlag)
When calculating the complete integrals K(k) and E(K), we can make complete true or false, except when abs(k') <= 1.0E-10, in which case complete must be true and the complementary modulus must be used The only merit of the alternative method over the standard method is that the code is simpler Both methods are based on modifications of Carlson's algortihm.

EllipticIntegral(double mod, double[] first, double[] second)
Used for complete elliptic integrals of first and second kind with real modulus 0 <= mod <= 1, mod is the modulus.

EllipticIntegral(double modr, double modi, boolean complementaryModulusUsed, boolean analyticContinuationUsed, double[] firstr, double[] firsti, double[] secondr, double[] secondi, int[] errorFlag)
Constructor only for complete integrals with phi = PI/2.

EllipticIntegral(double mod, double phi, double[] first, double[] second)
Used for complete and incomplete elliptic integrals of first and second kind with real modulus and argument 0 <= mod <= 1, mod is the modulus phi, argument in radians.

EllipticIntegral(double mod, double phi, double c, double[] third)
Used for complete and incomplete elliptic integrals of the third kind with real modulus, real argument phi, and real parameter c 0 <= mod <= 1 phi, argument in radians c parameter, 0 <= c <=1.

EllipticIntegral(double x, double y, double kc, double a, double b, double[] u, double[] v)
Creates a new EllipticIntegral object.

Method Summary

void	run() DOCUMENT ME!

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Constructor Detail

EllipticIntegral

public EllipticIntegral()

~ Constructors ---------------------------------------------------------------------------------------------------

EllipticIntegral

public EllipticIntegral(double mod,
                        double[] first,
                        double[] second)

Used for complete elliptic integrals of first and second kind with real modulus 0 <= mod <= 1, mod is the modulus.

Parameters:

mod - DOCUMENT ME!

first - DOCUMENT ME!

second - DOCUMENT ME!

EllipticIntegral

public EllipticIntegral(double mod,
                        double phi,
                        double[] first,
                        double[] second)

Used for complete and incomplete elliptic integrals of first and second kind with real modulus and argument 0 <= mod <= 1, mod is the modulus phi, argument in radians.

Parameters:

mod - DOCUMENT ME!

phi - DOCUMENT ME!

first - DOCUMENT ME!

second - DOCUMENT ME!

EllipticIntegral

public EllipticIntegral(double mod,
                        double phi,
                        double c,
                        double[] third)

Used for complete and incomplete elliptic integrals of the third kind with real modulus, real argument phi, and real parameter c 0 <= mod <= 1 phi, argument in radians c parameter, 0 <= c <=1.

Parameters:

mod - DOCUMENT ME!

phi - DOCUMENT ME!

c - DOCUMENT ME!

third - DOCUMENT ME!

EllipticIntegral

public EllipticIntegral(double x,
                        double y,
                        double kc,
                        double a,
                        double b,
                        double[] u,
                        double[] v)

Creates a new EllipticIntegral object.

Parameters:

x - DOCUMENT ME!

y - DOCUMENT ME!

kc - DOCUMENT ME!

a - DOCUMENT ME!

b - DOCUMENT ME!

u - DOCUMENT ME!

v - DOCUMENT ME!

EllipticIntegral

public EllipticIntegral(double modr,
                        double modi,
                        boolean complementaryModulusUsed,
                        boolean analyticContinuationUsed,
                        double[] firstr,
                        double[] firsti,
                        double[] secondr,
                        double[] secondi,
                        int[] errorFlag)

Constructor only for complete integrals with phi = PI/2. The method called is based on the arithmetic-geometric mean procedure If only the complete integrals are calculated, the arithmetic-geometric mean procdure is superior to the 2 methods based on Carlson's algorithm

Parameters:

modr - DOCUMENT ME!

modi - DOCUMENT ME!

complementaryModulusUsed - DOCUMENT ME!

analyticContinuationUsed - DOCUMENT ME!

firstr - DOCUMENT ME!

firsti - DOCUMENT ME!

secondr - DOCUMENT ME!

secondi - DOCUMENT ME!

errorFlag - DOCUMENT ME!

EllipticIntegral

public EllipticIntegral(boolean complete,
                        double modr,
                        double modi,
                        boolean complementaryModulusUsed,
                        boolean analyticContinuationUsed,
                        double phir,
                        double phii,
                        double[] firstr,
                        double[] firsti,
                        double[] secondr,
                        double[] secondi,
                        boolean useStandardMethod,
                        int[] errorFlag)

When calculating the complete integrals K(k) and E(K), we can make complete true or false, except when abs(k') <= 1.0E-10, in which case complete must be true and the complementary modulus must be used The only merit of the alternative method over the standard method is that the code is simpler Both methods are based on modifications of Carlson's algortihm.

Parameters:

complete - DOCUMENT ME!

modr - DOCUMENT ME!

modi - DOCUMENT ME!

complementaryModulusUsed - DOCUMENT ME!

analyticContinuationUsed - DOCUMENT ME!

phir - DOCUMENT ME!

phii - DOCUMENT ME!

firstr - DOCUMENT ME!

firsti - DOCUMENT ME!

secondr - DOCUMENT ME!

secondi - DOCUMENT ME!

useStandardMethod - DOCUMENT ME!

errorFlag - DOCUMENT ME!

Method Detail

run

public void run()

DOCUMENT ME!