Skip to content

Elliptic integrals

ripplegw.waveforms.elliptic_integrals ¤

JAX implementation of GSL elliptic integral functions.

This module provides JAX-compatible implementations of elliptic integrals, specifically the incomplete elliptic integral of the first kind (F).

ellint_F(phi: Float, k: Float, n_points: int = 1000) -> Float ¤

Compute the incomplete elliptic integral of the first kind.

This function computes F(φ, k) using the Legendre form: F(φ, k) = ∫₀^φ dt / √(1 - k² sin²(t))

This is equivalent to GSL's gsl_sf_ellint_F(phi, k, GSL_PREC_DOUBLE).

Parameters:

Name Type Description Default
phi Float

The amplitude (upper limit of integration) in radians

required
k Float

The modulus, where 0 ≤ k² ≤ 1 (note: this is k, not m = k²)

required
n_points int

Number of integration points (default: 1000)

1000

Returns:

Type Description
Float

The value of F(φ, k)

Notes
  • For k = 0: F(φ, 0) = φ
  • For |k| = 1 and |φ| < π/2: F(φ, ±1) = arctanh(sin(φ))
  • The implementation uses numerical integration via trapezoidal rule
  • JAX-compatible (can be used in jit, grad, vmap, etc.)

Examples:

>>> ellint_F(jnp.pi/2, 0.5)  # Complete elliptic integral K(0.5)
>>> ellint_F(jnp.pi/4, 0.8)  # Incomplete elliptic integral

ellint_Kcomp(k: Float, n_points: int = 1000) -> Float ¤

Compute the complete elliptic integral of the first kind.

This is K(k) = F(π/2, k) = ∫₀^(π/2) dt / √(1 - k² sin²(t))

Equivalent to GSL's gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE).

Parameters:

Name Type Description Default
k Float

The modulus, where 0 ≤ k² ≤ 1

required
n_points int

Number of integration points (default: 1000)

1000

Returns:

Type Description
Float

The value of K(k)

ellint_F_carlson(phi: Float, k: Float) -> Float ¤

Compute the incomplete elliptic integral of the first kind using Carlson's method.

This is an alternative implementation using Carlson symmetric form: F(φ, k) = sin(φ) * R_F(cos²(φ), 1 - k² sin²(φ), 1)

This method can be more accurate for certain parameter ranges but requires implementing Carlson's R_F function.

Parameters:

Name Type Description Default
phi Float

The amplitude in radians

required
k Float

The modulus

required

Returns:

Type Description
Float

The value of F(φ, k)

Note

This is a placeholder for a future implementation using Carlson's symmetric elliptic integrals, which are numerically more stable.

gsl_sf_elljac_e(u: Float, m: Float, max_iter: int = 16) ¤

Compute the Jacobian elliptic functions sn(u|m), cn(u|m), dn(u|m).

This function computes the Jacobian elliptic functions using the descending Landen transformation, which is the same algorithm used by GSL.

The Jacobian elliptic functions are defined by the inverse of the elliptic integral of the first kind: u = F(φ, k) = ∫₀^φ dt / √(1 - k² sin²(t))

Then

sn(u|m) = sin(φ) cn(u|m) = cos(φ) dn(u|m) = √(1 - k² sin²(φ))

where m = k² is the parameter (0 ≤ m ≤ 1).

This is equivalent to GSL's gsl_sf_elljac_e(u, m, &sn, &cn, &dn).

Parameters:

Name Type Description Default
u Float

The argument

required
m Float

The parameter (m = k², where k is the modulus), 0 ≤ m ≤ 1

required
max_iter int

Maximum number of Landen transformations (default: 16)

16

Returns:

Type Description

A tuple (sn, cn, dn) containing the three Jacobian elliptic functions

Notes
  • For m = 0: sn(u|0) = sin(u), cn(u|0) = cos(u), dn(u|0) = 1
  • For m = 1: sn(u|1) = tanh(u), cn(u|1) = sech(u), dn(u|1) = sech(u)
  • The implementation uses the descending Landen transformation for accuracy
  • JAX-compatible (can be used in jit, grad, vmap, etc.)
References
  • Abramowitz & Stegun, Chapter 16
  • GSL library: gsl_sf_elljac.c
  • NIST DLMF: https://dlmf.nist.gov/22.20

Examples:

>>> sn, cn, dn = gsl_sf_elljac_e(1.0, 0.5)
>>> sn, cn, dn = gsl_sf_elljac_e(jnp.array([0.5, 1.0]), 0.8)