Package 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: 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 invalid input: '&' Sons, Inc., 1996, Chapter 18, Elliptic Integrals and Jacobian Elliptic Functions, pp. 654-679.

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 and ZHANG elit give the same answer.

On 2016/04/09 0:05, Gandler, William (NIH/CIT) [E] wrote: > Dear Interdisciplinary Information Sciences: > I am a computer programmer at the National Institutes of Health in Bethesda, Maryland, USA. I would like to incorporate elliptic integral calculation routines into the MIPAV software package. These routines are in Figures 1 and 2 of "Numerical Calculation of the Elliptic Integrals of the First and Second Kind with Complex Modulus" by Tohru Morita, Interdisciplinary Information Sciences, Vol. 6, No. 1, 2000, pp. 67-74. Who should I ask for permission? > > Sincerely, > > William Gandler Dear Dr William Gandler, Thank you very much for your inquiry. You may use the calculation routines published in our journal with clear statement of attribution. Sincerely, Nobuaki Obata Editor-in-Chief IIS

Field Summary

Fields

Modifier and Type	Field	Description
private double[]	aki	DOCUMENT ME!
private double[]	akr	DOCUMENT ME!
private boolean	analyticContinuationUsed	If false, function used If true, analytic continuation of function used.
private double	BIGd	DOCUMENT ME!
private double	c	DOCUMENT ME!
private double	C1	private double BIG = 3.0E137; // 3.0E37.
private double	C1d	DOCUMENT ME!
private double	C2	DOCUMENT ME!
private double	C2d	DOCUMENT ME!
private double	C3	DOCUMENT ME!
private double	C3d	DOCUMENT ME!
private double	C4	DOCUMENT ME!
private double	C4d	DOCUMENT ME!
private double	C5d	DOCUMENT ME!
private double	C6d	DOCUMENT ME!
private double[]	cki	DOCUMENT ME!
private double[]	ckr	DOCUMENT ME!
private boolean	complementaryModulusUsed	k is the modulus k', the complementary modulus = sqrt(1 - k*k).
private boolean	complete	DOCUMENT ME!
private int[]	errorFlag	DOCUMENT ME!
private double	ERRTOL	DOCUMENT ME!
private double	ERRTOLb	DOCUMENT ME!
private double[]	first	DOCUMENT ME!
private double[]	firsti	DOCUMENT ME!
private double[]	firstr	DOCUMENT ME!
private double	mod	modulus, 0 invalid input: '<'= mod invalid input: '<'= 1.
private double	modi	Imaginary part of modulus or complementary modulus.
private double	modr	Real part of modulus or complementary modulus.
private static final int	MORITA_SOURCE	F(k, phi) INCOMPLETE1 E(k, phi) INCOMPLETE2 K(k) COMPLETE1 E(k) COMPLETE2 Complex arguments for first and second complete and incomplete elliptic integrals.
private int	nrep	DOCUMENT ME!
private double	phi	DOCUMENT ME!
private double	phii	DOCUMENT ME!
private double	phir	For complete integrals phir = PI/2, phii = 0.
private double	pihf	DOCUMENT ME!
private double[]	rdi	DOCUMENT ME!
private double[]	rdr	DOCUMENT ME!
private double	RERR	DOCUMENT ME!
private double[]	rfi	DOCUMENT ME!
private double[]	rfr	DOCUMENT ME!
private boolean	runcel12	DOCUMENT ME!
private boolean	runcomelp	DOCUMENT ME!
private boolean	runelit	DOCUMENT ME!
private double[]	second	DOCUMENT ME!
private double[]	secondi	DOCUMENT ME!
private double[]	secondr	DOCUMENT ME!
private int	sgnck	DOCUMENT ME!
private int	source	DOCUMENT ME!
private static final int	TEST_SOURCE	Compare outputs of different modules for self-testing.
private double[]	third	DOCUMENT ME!
private double	TINY	DOCUMENT ME!
private double	TINYd	DOCUMENT ME!
private ViewUserInterface	UI	DOCUMENT ME!
private boolean	useStandardMethod	DOCUMENT ME!
private static final int	ZHANG_SOURCE	Real arguments for first and second complete and incomplete elliptic integrals.

Constructor Summary

Constructors

Constructor	Description
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') invalid input: '<'= 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 invalid input: '<'= mod invalid input: '<'= 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 invalid input: '<'= mod invalid input: '<'= 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 invalid input: '<'= mod invalid input: '<'= 1 phi, argument in radians c parameter, 0 invalid input: '<'= c invalid input: '<'=1.

Method Summary

Modifier and Type	Method	Description
private void	cel12(double ckr, double cki)	Routine cel12 only calculates complete integrals of the first and second kind with complex arguments Based on the arithmetic-geometric mean procedure Ported from Tohru Morita article.
private void	comelp(double mod)	Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete first and second elliptic integrals with a real modulus 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus.
private double	cosh(double x)	DOCUMENT ME!
private void	elit(double mod, double phi)	Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete and incomplete first and second elliptic integrals with a real modulus and real argument 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus arg is the argument in radians.
private void	elit3(double mod, double phi, double c)	Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete and incomplete elliptic integrals of the third kind with a real modulus real argument, and real parameter 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus arg is the argument in radians c is the parameter, 0 invalid input: '<'= c invalid input: '<'= 1.
private void	ell12(double phir, double phii, double akr, double aki, double ckr, double cki, int sgnck)	Calculates complete and incomplete elliptic integrals of the first and second kind with complex arguments Ported from the Tohru Morita article.
private void	rfrd(double xrtr, double xrti, double yrtr, double yrti)	Ported from Tohru Morita article.
private void	rfrdb(double xrtr, double xrti, double yrtr, double yrti)	Ported from Tohru Morita article.
void	run()	DOCUMENT ME!
private void	runTest()	DOCUMENT ME!
private double	sinh(double x)	DOCUMENT ME!
private void	sqrtc(double zinr, double zini, double[] zsqr, double[] zsqi)	Ported from Tohru Morita article Performs square root of complex argument.
private double	zabs(double zr, double zi)	zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.
private void	zcos(double inr, double ini, double[] outr, double[] outi)	DOCUMENT ME!
private void	zdiv(double ar, double ai, double br, double bi, double[] cr, double[] ci)	complex divide c = a/b.
private void	zmlt(double ar, double ai, double br, double bi, double[] cr, double[] ci)	complex multiply c = a * b.
private void	zsin(double inr, double ini, double[] outr, double[] outi)	DOCUMENT ME!

Methods inherited from class java.lang.Object

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

Field Details

MORITA_SOURCE

private static final int MORITA_SOURCE

F(k, phi) INCOMPLETE1 E(k, phi) INCOMPLETE2 K(k) COMPLETE1 E(k) COMPLETE2 Complex arguments for first and second complete and incomplete elliptic integrals.

See Also:

• Constant Field Values

ZHANG_SOURCE

private static final int ZHANG_SOURCE

Real arguments for first and second complete and incomplete elliptic integrals.

See Also:

• Constant Field Values

TEST_SOURCE

private static final int TEST_SOURCE

Compare outputs of different modules for self-testing.

See Also:

• Constant Field Values

aki

private double[] aki

DOCUMENT ME!

akr

private double[] akr

DOCUMENT ME!

analyticContinuationUsed

private boolean analyticContinuationUsed

If false, function used If true, analytic continuation of function used.

BIGd

private double BIGd

DOCUMENT ME!

c

private double c

DOCUMENT ME!

C1

private double C1

private double BIG = 3.0E137; // 3.0E37.

C1d

private double C1d

DOCUMENT ME!

C2

private double C2

DOCUMENT ME!

C2d

private double C2d

DOCUMENT ME!

C3

private double C3

DOCUMENT ME!

C3d

private double C3d

DOCUMENT ME!

C4

private double C4

DOCUMENT ME!

C4d

private double C4d

DOCUMENT ME!

C5d

private double C5d

DOCUMENT ME!

C6d

private double C6d

DOCUMENT ME!

cki

private double[] cki

DOCUMENT ME!

ckr

private double[] ckr

DOCUMENT ME!

complementaryModulusUsed

private boolean complementaryModulusUsed

k is the modulus k', the complementary modulus = sqrt(1 - k*k).

complete

private boolean complete

DOCUMENT ME!

errorFlag

private int[] errorFlag

DOCUMENT ME!

ERRTOL

private double ERRTOL

DOCUMENT ME!

ERRTOLb

private double ERRTOLb

DOCUMENT ME!

first

private double[] first

DOCUMENT ME!

firsti

private double[] firsti

DOCUMENT ME!

firstr

private double[] firstr

DOCUMENT ME!

mod

private double mod

modulus, 0 invalid input: '<'= mod invalid input: '<'= 1.

modi

private double modi

Imaginary part of modulus or complementary modulus.

modr

private double modr

Real part of modulus or complementary modulus.

nrep

private int nrep

DOCUMENT ME!

phi

private double phi

DOCUMENT ME!

phii

private double phii

DOCUMENT ME!

phir

private double phir

For complete integrals phir = PI/2, phii = 0.

pihf

private double pihf

DOCUMENT ME!

rdi

private double[] rdi

DOCUMENT ME!

rdr

private double[] rdr

DOCUMENT ME!

RERR

private double RERR

DOCUMENT ME!

rfi

private double[] rfi

DOCUMENT ME!

rfr

private double[] rfr

DOCUMENT ME!

runcel12

private boolean runcel12

DOCUMENT ME!

runcomelp

private boolean runcomelp

DOCUMENT ME!

runelit

private boolean runelit

DOCUMENT ME!

second

private double[] second

DOCUMENT ME!

secondi

private double[] secondi

DOCUMENT ME!

secondr

private double[] secondr

DOCUMENT ME!

sgnck

private int sgnck

DOCUMENT ME!

source

private int source

DOCUMENT ME!

third

private double[] third

DOCUMENT ME!

TINY

private double TINY

DOCUMENT ME!

TINYd

private double TINYd

DOCUMENT ME!

UI

private ViewUserInterface UI

DOCUMENT ME!

useStandardMethod

private boolean useStandardMethod

DOCUMENT ME!

Constructor Details

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 invalid input: '<'= mod invalid input: '<'= 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 invalid input: '<'= mod invalid input: '<'= 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 invalid input: '<'= mod invalid input: '<'= 1 phi, argument in radians c parameter, 0 invalid input: '<'= c invalid input: '<'=1.

Parameters:

mod - DOCUMENT ME!

phi - DOCUMENT ME!

c - DOCUMENT ME!

third - 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') invalid input: '<'= 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 Details

run

public void run()

DOCUMENT ME!

cel12

private void cel12(double ckr,
 double cki)

Routine cel12 only calculates complete integrals of the first and second kind with complex arguments Based on the arithmetic-geometric mean procedure Ported from Tohru Morita article.

Parameters:

ckr - DOCUMENT ME!

cki - DOCUMENT ME!

comelp

private void comelp(double mod)

Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete first and second elliptic integrals with a real modulus 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus.

Parameters:

mod - DOCUMENT ME!

cosh

private double cosh(double x)

DOCUMENT ME!

Parameters:

x - DOCUMENT ME!

Returns:

DOCUMENT ME!

elit

private void elit(double mod,
 double phi)

Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete and incomplete first and second elliptic integrals with a real modulus and real argument 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus arg is the argument in radians.

Parameters:

mod - DOCUMENT ME!

phi - DOCUMENT ME!

elit3

private void elit3(double mod,
 double phi,
 double c)

Ported from Computation of Special Functions by Shanjie Zhang and Jianming Jin Computes complete and incomplete elliptic integrals of the third kind with a real modulus real argument, and real parameter 0 invalid input: '<'= mod invalid input: '<'= 1, mod is the modulus arg is the argument in radians c is the parameter, 0 invalid input: '<'= c invalid input: '<'= 1.

Parameters:

mod - DOCUMENT ME!

phi - DOCUMENT ME!

c - DOCUMENT ME!

ell12

private void ell12(double phir,
 double phii,
 double akr,
 double aki,
 double ckr,
 double cki,
 int sgnck)

Calculates complete and incomplete elliptic integrals of the first and second kind with complex arguments Ported from the Tohru Morita article.

Parameters:

phir - DOCUMENT ME!

phii - DOCUMENT ME!

akr - DOCUMENT ME!

aki - DOCUMENT ME!

ckr - DOCUMENT ME!

cki - DOCUMENT ME!

sgnck - DOCUMENT ME!

rfrd

private void rfrd(double xrtr,
 double xrti,
 double yrtr,
 double yrti)

Ported from Tohru Morita article.

Parameters:

xrtr - DOCUMENT ME!

xrti - DOCUMENT ME!

yrtr - DOCUMENT ME!

yrti - DOCUMENT ME!

rfrdb

private void rfrdb(double xrtr,
 double xrti,
 double yrtr,
 double yrti)

Ported from Tohru Morita article.

Parameters:

xrtr - DOCUMENT ME!

xrti - DOCUMENT ME!

yrtr - DOCUMENT ME!

yrti - DOCUMENT ME!

runTest

private void runTest()

DOCUMENT ME!

sinh

private double sinh(double x)

DOCUMENT ME!

Parameters:

x - DOCUMENT ME!

Returns:

DOCUMENT ME!

sqrtc

private void sqrtc(double zinr,
 double zini,
 double[] zsqr,
 double[] zsqi)

Ported from Tohru Morita article Performs square root of complex argument.

Parameters:

zinr - DOCUMENT ME!

zini - DOCUMENT ME!

zsqr - DOCUMENT ME!

zsqi - DOCUMENT ME!

zabs

private double zabs(double zr,
 double zi)

zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.

Parameters:

zr - double

zi - double

Returns:

double

zcos

private void zcos(double inr,
 double ini,
 double[] outr,
 double[] outi)

DOCUMENT ME!

Parameters:

inr - DOCUMENT ME!

ini - DOCUMENT ME!

outr - DOCUMENT ME!

outi - DOCUMENT ME!

zdiv

private void zdiv(double ar,
 double ai,
 double br,
 double bi,
 double[] cr,
 double[] ci)

complex divide c = a/b.

Parameters:

ar - double

ai - double

br - double

bi - double

cr - double[]

ci - double[]

zmlt

private void zmlt(double ar,
 double ai,
 double br,
 double bi,
 double[] cr,
 double[] ci)

complex multiply c = a * b.

Parameters:

ar - double

ai - double

br - double

bi - double

cr - double[]

ci - double[]

zsin

private void zsin(double inr,
 double ini,
 double[] outr,
 double[] outi)

DOCUMENT ME!

Parameters:

inr - DOCUMENT ME!

ini - DOCUMENT ME!

outr - DOCUMENT ME!

outi - DOCUMENT ME!