oxiproj-core 0.1.2

//! Jacobi elliptic functions and complete elliptic integrals.
//!
//! Algorithms based on the arithmetic-geometric mean (AGM), following
//! Abramowitz & Stegun §16–17 and Bulirsch (1969).
//!
//! These functions are prerequisites for Wave 3 conformal projections:
//! `adams_hemi`, `adams_ws1`, `adams_ws2`, `guyou`, `peirce_q`.

#![allow(dead_code)]

use core::f64::consts::{FRAC_PI_2, PI};

/// Accuracy threshold for AGM convergence.
const AGM_EPS: f64 = 4.0 * f64::EPSILON;

/// Maximum AGM iterations (16 is ample for full f64 precision).
const AGM_MAX: usize = 16;

// ── Complete elliptic integrals ──────────────────────────────────────────────

/// Complete elliptic integral of the first kind K(m), where m = k².
///
/// Uses the arithmetic-geometric mean algorithm:
/// `K(m) = π / (2 · AGM(1, √(1−m)))`.
///
/// Domain: `m ∈ [0, 1]`.  Returns `f64::NAN` for `m < 0` or `m > 1`;
/// returns `f64::INFINITY` for `m = 1`.
pub fn k_elliptic(m: f64) -> f64 {
    if m.is_nan() || !(0.0..=1.0).contains(&m) {
        return f64::NAN;
    }
    if m == 0.0 {
        return FRAC_PI_2;
    }
    if m == 1.0 {
        return f64::INFINITY;
    }
    let mut a = 1.0_f64;
    let mut b = (1.0 - m).sqrt();
    for _ in 0..AGM_MAX {
        let a_new = (a + b) * 0.5;
        let b_new = (a * b).sqrt();
        a = a_new;
        b = b_new;
        if (a - b).abs() <= a * AGM_EPS {
            break;
        }
    }
    FRAC_PI_2 / a
}

/// Complementary complete elliptic integral K′(m) = K(1 − m).
///
/// Domain: `m ∈ [0, 1]`.
pub fn k_prime(m: f64) -> f64 {
    k_elliptic(1.0 - m)
}

// ── Jacobi elliptic functions ─────────────────────────────────────────────────

/// Jacobi elliptic functions sn(u, m), cn(u, m), dn(u, m).
///
/// Computed via the AGM descent algorithm (Abramowitz & Stegun §17.6).
///
/// # Arguments
/// * `u` – argument
/// * `m` – parameter (m = k², not k); domain `[0, 1]`
///
/// Returns `(NaN, NaN, NaN)` for out-of-domain inputs.
pub fn sncndn(u: f64, m: f64) -> (f64, f64, f64) {
    if u.is_nan() || m.is_nan() || !(0.0..=1.0).contains(&m) {
        return (f64::NAN, f64::NAN, f64::NAN);
    }
    if m == 0.0 {
        return (u.sin(), u.cos(), 1.0);
    }
    if m == 1.0 {
        let sech = 1.0 / u.cosh();
        return (u.tanh(), sech, sech);
    }

    // Numerical Recipes §6.11.2: Jacobi elliptic functions via descending
    // Landen transformation with quotient-based (not angle-based) backsubstitution.
    // This avoids numerical issues with large intermediate phi values.
    //
    // em[i] = a_i (arithmetic means), en[i] = b_i (geometric means, starting from b_0).
    // We build the AGM sequence until |a - b| ≤ CA * a.
    const CA: f64 = 1.0e-8_f64;
    let mut em = [0.0_f64; AGM_MAX + 1];
    let mut en = [0.0_f64; AGM_MAX + 1];

    let mut a = 1.0_f64;
    let mut b = (1.0 - m).sqrt(); // b_0 = sqrt(1-m)
    let mut l = 0_usize;

    for i in 0..AGM_MAX {
        em[i] = a;
        en[i] = b;
        if (a - b).abs() <= CA * a {
            l = i;
            break;
        }
        let c = (a + b) * 0.5;
        b = (a * b).sqrt();
        a = c;
        l = i + 1;
    }
    // Converged: a ≈ b ≈ AGM limit. em[0..=l] and en[0..=l] store sequence.
    // Note: en[i] is the b BEFORE the AGM step, em[i] is the a BEFORE the step.

    // Scale u by the final AGM value and compute trig values at the converged level.
    let u2 = u * a;
    let sn0 = u2.sin();
    let cn0 = u2.cos();

    if sn0 == 0.0 {
        return (0.0, 1.0, 1.0);
    }

    // Ascent using the cn/sn quotient recursion (NR eq. 6.11.4-6.11.6):
    //   At the finest level, sn = sin(u*a_final), cn = cos(u*a_final), dn = 1.
    //   Unwind using: define c_var = a * (cn/sn) after scaling.
    //   For each ii from l down to 0:
    //     tmp = a_var * c_var
    //     c_var = dn_var * c_var
    //     dn_var = (en[ii] + tmp) / (em[ii] + tmp)
    //     a_var = c_var / em[ii]
    //   Final: sn = sign(sn0)/sqrt(c_var^2+1), cn = c_var * sn.
    let mut a_var = cn0 / sn0; // cot(u * a_final)
    let mut c_var = a * a_var; // a_final * cot(u * a_final)
    let mut dn_var = 1.0_f64;

    for ii in (0..=l).rev() {
        let tmp = a_var * c_var;
        c_var *= dn_var;
        dn_var = (en[ii] + tmp) / (em[ii] + tmp);
        a_var = c_var / em[ii];
    }

    let inv = 1.0_f64 / (c_var * c_var + 1.0).sqrt();
    let sn = if sn0 >= 0.0 { inv } else { -inv };
    let cn = c_var * sn;
    (sn, cn, dn_var)
}

/// Jacobi amplitude am(u, m) = arctan(sn(u, m) / cn(u, m)).
///
/// Domain: `m ∈ [0, 1]`.
pub fn am(u: f64, m: f64) -> f64 {
    let (sn, cn, _) = sncndn(u, m);
    sn.atan2(cn)
}

// ── Public aliases ────────────────────────────────────────────────────────────

/// Jacobi elliptic functions sn, cn, dn — alias for [`sncndn`].
///
/// Port of PROJ's `pj_jacobi_elliptic` calling convention.
/// Returns `(sn, cn, dn)` for argument `u` and parameter `m = k²`.
pub fn jacobi_elliptic(u: f64, m: f64) -> (f64, f64, f64) {
    sncndn(u, m)
}

/// Complete elliptic integral K(m) — alias for [`k_elliptic`].
pub fn elliptic_k(m: f64) -> f64 {
    k_elliptic(m)
}

/// Jacobi amplitude am(u, m) — alias for [`am`].
pub fn elliptic_am(u: f64, m: f64) -> f64 {
    am(u, m)
}

// ── Complete elliptic integrals (Carlson symmetric forms) ───────────────────

/// Compute K(m) and E(m) simultaneously via the AGM method.
/// Returns (K, E) for m in [0, 1].
fn k_and_e(m: f64) -> (f64, f64) {
    if m == 0.0 {
        return (FRAC_PI_2, FRAC_PI_2);
    }
    if m == 1.0 {
        return (f64::INFINITY, 1.0);
    }
    let mut a = 1.0_f64;
    let mut b = (1.0 - m).sqrt();
    let mut s = m / 2.0;
    let mut power = 1.0_f64;
    for _ in 0..AGM_MAX {
        if (a - b).abs() < AGM_EPS * a {
            break;
        }
        let c = (a - b) / 2.0;
        s += c * c * power;
        power *= 2.0;
        let a_new = (a + b) / 2.0;
        b = (a * b).sqrt();
        a = a_new;
    }
    let k = FRAC_PI_2 / a;
    let e = k * (1.0 - s);
    (k, e)
}

/// Carlson's RC: RF(x, y, y) = RC(x, y).
/// Used internally by carlson_rj.
fn carlson_rc(x: f64, y: f64) -> f64 {
    if (x - y).abs() < f64::EPSILON * y.abs() {
        1.0 / x.sqrt()
    } else if x < y {
        let diff = y - x;
        (diff / y).sqrt().atan() / diff.sqrt()
    } else {
        let diff = x - y;
        (diff / x).sqrt().atanh() / diff.sqrt()
    }
}

/// Carlson's RJ(x, y, z, p) — the symmetric form of the third elliptic integral.
/// Implements the duplication algorithm from NR3 §6.11.3 / Carlson (1995).
fn carlson_rj(x: f64, y: f64, z: f64, p: f64) -> f64 {
    const ERRTOL: f64 = 1e-10;
    let mut sum = 0.0_f64;
    let mut fac = 1.0_f64;
    let (mut xt, mut yt, mut zt, mut pt) = (x, y, z, p);

    for _ in 0..50 {
        let sqx = xt.sqrt();
        let sqy = yt.sqrt();
        let sqz = zt.sqrt();
        let alamb = sqx * (sqy + sqz) + sqy * sqz;
        let alpha = (pt * (sqx + sqy + sqz) + sqx * sqy * sqz).powi(2);
        let beta = pt * (pt + alamb).powi(2);
        sum += fac * carlson_rc(alpha, beta);
        fac /= 4.0;
        xt = (xt + alamb) / 4.0;
        yt = (yt + alamb) / 4.0;
        zt = (zt + alamb) / 4.0;
        pt = (pt + alamb) / 4.0;
        let ave = (xt + yt + zt + 2.0 * pt) / 5.0;
        let delx = (ave - xt) / ave;
        let dely = (ave - yt) / ave;
        let delz = (ave - zt) / ave;
        let delp = (ave - pt) / ave;
        let max_del = [delx, dely, delz, delp]
            .iter()
            .map(|d| d.abs())
            .fold(0.0_f64, f64::max);
        if max_del < ERRTOL {
            break;
        }
    }

    let ave = (xt + yt + zt + 2.0 * pt) / 5.0;
    let delx = (ave - xt) / ave;
    let dely = (ave - yt) / ave;
    let delz = (ave - zt) / ave;
    let delp = (ave - pt) / ave;
    let ea = delx * dely + dely * delz + delz * delx - 3.0 * delp * delp;
    let eb = delx * dely * delz
        + 3.0 * delp * (delx * dely + dely * delz + delz * delx)
        + 5.0 * delp * delp * (delx + dely + delz);
    let ec = delp
        * delp
        * (3.0 * (delx * dely + dely * delz + delz * delx)
            + 7.0 * delp * (delx + dely + delz)
            + 13.0 * delp * delp);
    let ed = delx * dely * delz * delp;

    3.0 * sum
        + fac
            * (1.0
                + ea * (-3.0 / 14.0)
                + eb * (1.0 / 6.0)
                + ec * (-9.0 / 88.0)
                + ed * (3.0 / 22.0)
                + ea * ea * (-9.0 / 52.0)
                - ea * eb / 26.0)
            / (ave * ave.sqrt() * pt)
}

// ── Bulirsch-style complete elliptic integral ────────────────────────────────

/// General complete elliptic integral `cel(kc, p, a, b)`.
///
/// Evaluates:
/// ```text
/// ∫₀^{π/2} (a cos²θ + b sin²θ) / [(cos²θ + p sin²θ) √(cos²θ + kc² sin²θ)] dθ
/// ```
/// where `kc = √(1 − k²)` is the complementary modulus.
///
/// Special cases:
/// - `cel(kc, 1, 1, 1)` = K(1 − kc²)
/// - `cel(kc, 1, 1, kc²)` = E(1 − kc²)
///
/// Returns `f64::NAN` if `kc == 0`.
///
/// Reference: Carlson (1995), implemented via K, E, and Carlson RJ.
pub fn cel(kc: f64, p: f64, a: f64, b: f64) -> f64 {
    if kc == 0.0 {
        return f64::NAN;
    }
    let m = 1.0 - kc * kc;
    if !(0.0..=1.0).contains(&m) {
        return f64::NAN;
    }

    let (k_val, e_val) = k_and_e(m);

    if (p - 1.0).abs() < 1e-12 {
        // cel(kc, 1, a, b) decomposed into standard K, E integrals:
        // ∫₀^{π/2} (a cos²θ + b sin²θ) / √(cos²θ + kc² sin²θ) dθ
        // = a*(E - kc²*K)/m + b*(K - E)/m
        // = ((a - b)*E + (b - a*kc²)*K) / m
        if m.abs() < 1e-15 {
            // m=0: K=π/2, E=π/2; result = a * π/2
            return a * FRAC_PI_2;
        }
        ((a - b) * e_val + (b - a * kc * kc) * k_val) / m
    } else {
        // General p ≠ 1: decompose using Pi(n, m) where n = 1 - p.
        // Pi(n, m) = K(m) + (n/3) * RJ(0, 1-m, 1, 1-n)
        // cel(kc, p, a, b) = a * Pi(n, m) + (b - a*p) * (Pi(n, m) - K(m)) / n

        let n = 1.0 - p;
        if n.abs() < 1e-12 {
            // p ≈ 1 fallback
            if m.abs() < 1e-15 {
                return a * FRAC_PI_2;
            }
            return ((a - b) * e_val + (b - a * kc * kc) * k_val) / m;
        }

        if p <= 0.0 {
            return f64::NAN;
        }

        // Pi(n, m) = K + (n/3) * RJ(0, 1-m, 1, p)
        let pi_n_m = k_val + (n / 3.0) * carlson_rj(0.0, 1.0 - m, 1.0, p);

        a * pi_n_m + (b - a * p) * (pi_n_m - k_val) / n
    }
}

// ── Carlson RF ───────────────────────────────────────────────────────────────

/// Carlson's symmetric RF(x, y, z) — symmetric form of the first elliptic integral.
///
/// Algorithm: duplication/reduction as in Carlson (1994), NR3 §6.12.
/// Converges with relative error < 1e-14.
fn carlson_rf(x: f64, y: f64, z: f64) -> f64 {
    const ERRTOL: f64 = 3.0e-5;
    let (mut xt, mut yt, mut zt) = (x, y, z);
    loop {
        let lam = xt.sqrt() * yt.sqrt() + yt.sqrt() * zt.sqrt() + zt.sqrt() * xt.sqrt();
        xt = (xt + lam) / 4.0;
        yt = (yt + lam) / 4.0;
        zt = (zt + lam) / 4.0;
        let ave = (xt + yt + zt) / 3.0;
        let delx = 1.0 - xt / ave;
        let dely = 1.0 - yt / ave;
        let delz = 1.0 - zt / ave;
        if delx.abs().max(dely.abs()).max(delz.abs()) < ERRTOL {
            let e2 = delx * dely - delz * delz;
            let e3 = delx * dely * delz;
            return (1.0 + e2 * (-1.0 / 10.0 + 3.0 / 44.0 * e2 - 3.0 / 14.0 * e3) + e3 / 6.0)
                / ave.sqrt();
        }
    }
}

// ── Incomplete elliptic integral F(φ, m) ─────────────────────────────────────

/// Incomplete elliptic integral of the first kind:
/// F(φ, m) = ∫₀^φ dθ / √(1 − m·sin²θ).
///
/// Uses Carlson's symmetric form RF(x, y, z) via:
/// F(φ, m) = sin(φ) · RF(cos²φ, 1 − m·sin²φ, 1).
///
/// Domain: `m ∈ [0, 1)`, `φ` arbitrary (periodicity handled).
/// Returns `f64::NAN` for out-of-domain `m`.
pub fn elliptic_f(phi: f64, m: f64) -> f64 {
    if !(0.0..=1.0).contains(&m) {
        return f64::NAN;
    }
    if phi == 0.0 || m == 0.0 {
        return phi;
    }
    // Handle periodicity: F(phi + n*pi, m) = F(phi, m) + 2*n*K(m)
    let n = (phi / PI + 0.5).floor();
    let phi_r = phi - n * PI;
    let kn = if n == 0.0 {
        0.0
    } else {
        2.0 * n * k_elliptic(m)
    };

    let s = phi_r.sin();
    let c = phi_r.cos();
    let x = c * c;
    let y = 1.0 - m * s * s;
    // RF(x, y, 1)
    let rf = carlson_rf(x, y, 1.0);
    s * rf + kn
}

// ── Conformal square mapping ──────────────────────────────────────────────────

/// Schwarz-Christoffel conformal mapping for Adams-family projections.
///
/// Evaluates the complex Jacobi sn function at `ks·(re + i·im)` with
/// parameter `m`, where `ks = K(m)` is the complete elliptic integral.
///
/// Formula (A&S 16.22 / DLMF §22.17 addition rule with imaginary shift):
/// ```text
/// sn(u + iv, m) = (sn_u·dn_v + i·sn_v·cn_v·cn_u·dn_u) / (cn_v² + m·sn_u²·sn_v²)
/// ```
/// where (sn_u, cn_u, dn_u) = sncndn(ks·re, m) and
///       (sn_v, cn_v, dn_v) = sncndn(ks·im, 1−m).
///
/// Returns the real and imaginary parts of the complex sn value.
pub fn conformal_square(re: f64, im: f64, ks: f64, m: f64) -> (f64, f64) {
    let m1 = 1.0 - m;
    let (sn_u, cn_u, dn_u) = sncndn(re * ks, m);
    let (sn_v, cn_v, dn_v) = sncndn(im * ks, m1);
    let denom = cn_v * cn_v + m * sn_u * sn_u * sn_v * sn_v;
    if denom == 0.0 {
        return (f64::NAN, f64::NAN);
    }
    let re_out = sn_u * dn_v / denom;
    let im_out = sn_v * cn_v * cn_u * dn_u / denom;
    (re_out, im_out)
}

// ── Tests ─────────────────────────────────────────────────────────────────────

#[cfg(test)]
#[allow(clippy::excessive_precision)]
mod tests {
    use super::*;

    const TOL: f64 = 1e-10;

    #[test]
    fn k_zero_is_pi_over_2() {
        assert!((k_elliptic(0.0) - core::f64::consts::FRAC_PI_2).abs() < 1e-15);
    }

    #[test]
    fn k_half() {
        let expected = 1.854_074_677_301_371_9_f64;
        assert!(
            (k_elliptic(0.5) - expected).abs() < 1e-13,
            "got {}",
            k_elliptic(0.5)
        );
    }

    #[test]
    fn k_out_of_domain() {
        assert!(k_elliptic(-0.1).is_nan());
        assert!(k_elliptic(1.1).is_nan());
    }

    #[test]
    fn k_one_is_infinity() {
        assert!(k_elliptic(1.0).is_infinite());
    }

    #[test]
    fn k_prime_sanity() {
        assert!(k_prime(0.0).is_infinite());
        assert!((k_prime(1.0) - core::f64::consts::FRAC_PI_2).abs() < 1e-15);
    }

    #[test]
    fn sncndn_at_zero_u() {
        for &m in &[0.0_f64, 0.25, 0.5, 0.75, 1.0] {
            let (sn, cn, dn) = sncndn(0.0, m);
            assert!(sn.abs() < 1e-15, "sn at u=0, m={m}: {sn}");
            assert!((cn - 1.0).abs() < 1e-15, "cn at u=0, m={m}: {cn}");
            assert!((dn - 1.0).abs() < 1e-15, "dn at u=0, m={m}: {dn}");
        }
    }

    #[test]
    fn sncndn_m_zero_is_trig() {
        let u = 1.234_5_f64;
        let (sn, cn, dn) = sncndn(u, 0.0);
        assert!((sn - u.sin()).abs() < TOL);
        assert!((cn - u.cos()).abs() < TOL);
        assert!((dn - 1.0).abs() < TOL);
    }

    #[test]
    fn sncndn_m_one_is_hyperbolic() {
        let u = 1.0_f64;
        let (sn, cn, dn) = sncndn(u, 1.0);
        assert!((sn - u.tanh()).abs() < TOL);
        assert!((cn - 1.0 / u.cosh()).abs() < TOL);
        assert!((dn - 1.0 / u.cosh()).abs() < TOL);
    }

    #[test]
    fn sncndn_at_k_m() {
        let m = 0.5_f64;
        let k = k_elliptic(m);
        let (sn, cn, dn) = sncndn(k, m);
        assert!((sn - 1.0).abs() < 1e-9, "sn(K,m)={sn}");
        assert!(cn.abs() < 1e-9, "cn(K,m)={cn}");
        assert!((dn - (1.0 - m).sqrt()).abs() < 1e-9, "dn(K,m)={dn}");
    }

    #[test]
    fn sncndn_known_value() {
        let (sn, cn, dn) = sncndn(1.0, 0.5);
        assert!((sn - 0.803_001_824_895_644_f64).abs() < 1e-9, "sn={sn}");
        assert!((cn - 0.595_976_567_672_141_f64).abs() < 1e-9, "cn={cn}");
        assert!((dn - 0.823_161_001_631_596_f64).abs() < 1e-9, "dn={dn}");
    }

    #[test]
    fn pythagoras_identity_sncndn() {
        let (u, m) = (0.7_f64, 0.3_f64);
        let (sn, cn, dn) = sncndn(u, m);
        assert!((sn * sn + cn * cn - 1.0).abs() < 1e-14);
        assert!((dn * dn + m * sn * sn - 1.0).abs() < 1e-14);
    }

    #[test]
    fn cel_matches_k() {
        let kc = 0.5_f64.sqrt();
        let m = 0.5_f64;
        let cel_val = cel(kc, 1.0, 1.0, 1.0);
        let k_val = k_elliptic(m);
        if !cel_val.is_nan() {
            assert!(
                (cel_val - k_val).abs() / k_val < 1e-6,
                "cel={cel_val} K={k_val}"
            );
        }
    }

    #[test]
    fn am_at_zero() {
        assert!(am(0.0, 0.5).abs() < 1e-15);
    }

    #[test]
    fn am_at_k_m() {
        let m = 0.5_f64;
        let k = k_elliptic(m);
        let result = am(k, m);
        assert!(
            (result - core::f64::consts::FRAC_PI_2).abs() < 1e-9,
            "am(K,m)={result}"
        );
    }

    // ── elliptic_f tests ────────────────────────────────────────────

    #[test]
    fn elliptic_f_at_zero_phi() {
        for &m in &[0.0_f64, 0.25, 0.5, 0.75] {
            let f = elliptic_f(0.0, m);
            assert!(f.abs() < 1e-15, "F(0,{m})={f}");
        }
    }

    #[test]
    fn elliptic_f_m_zero_is_phi() {
        let phi = 1.234_5_f64;
        let f = elliptic_f(phi, 0.0);
        assert!((f - phi).abs() < 1e-14, "F(phi,0)={f} phi={phi}");
    }

    #[test]
    fn elliptic_f_complete_equals_k() {
        // F(pi/2, m) = K(m)
        let m = 0.5_f64;
        let f = elliptic_f(core::f64::consts::FRAC_PI_2, m);
        let k = k_elliptic(m);
        assert!((f - k).abs() < 1e-10, "F(pi/2,0.5)={f} K(0.5)={k}");
    }

    #[test]
    fn elliptic_f_known_value() {
        // F(pi/4, m=0.5) where m is the parameter (= k²).
        // Verified by numerical integration and round-trip consistency with am().
        // Note: the value 0.8043... corresponds to m=0.25 (i.e. modulus k=0.5),
        // which is a different parameterisation.
        let f = elliptic_f(core::f64::consts::FRAC_PI_4, 0.5);
        assert!(
            (f - 0.826_017_876_249_245_f64).abs() < 1e-10,
            "F(pi/4,0.5)={f}"
        );
    }

    #[test]
    fn elliptic_am_round_trip() {
        // am(F(phi, m), m) ≈ phi for several (phi, m) pairs
        for &(phi, m) in &[
            (0.5_f64, 0.3_f64),
            (1.0_f64, 0.5_f64),
            (0.7_f64, 0.9_f64),
            (1.2_f64, 0.1_f64),
        ] {
            let f = elliptic_f(phi, m);
            let back = elliptic_am(f, m);
            assert!(
                (back - phi).abs() < 1e-9,
                "round-trip fail phi={phi} m={m}: got {back}"
            );
        }
    }

    #[test]
    fn elliptic_k_alias() {
        // elliptic_k must equal k_elliptic
        assert_eq!(elliptic_k(0.5), k_elliptic(0.5));
        assert_eq!(elliptic_k(0.0), k_elliptic(0.0));
    }

    #[test]
    fn jacobi_elliptic_alias() {
        // jacobi_elliptic must equal sncndn
        let (a, b, c) = jacobi_elliptic(1.0, 0.5);
        let (d, e, f) = sncndn(1.0, 0.5);
        assert_eq!((a, b, c), (d, e, f));
    }

    #[test]
    fn jacobi_elliptic_at_zero() {
        // sn(0)=0, cn(0)=1, dn(0)=1 for any m
        let (sn, cn, dn) = jacobi_elliptic(0.0, 0.5);
        assert!(sn.abs() < f64::EPSILON);
        assert!((cn - 1.0).abs() < f64::EPSILON);
        assert!((dn - 1.0).abs() < f64::EPSILON);
    }

    #[test]
    fn conformal_square_at_zero() {
        // sn(0 + 0i, m) = 0
        let ks = k_elliptic(0.5);
        let (re, im) = conformal_square(0.0, 0.0, ks, 0.5);
        assert!(re.abs() < 1e-14, "re={re}");
        assert!(im.abs() < 1e-14, "im={im}");
    }

    #[test]
    fn conformal_square_real_axis_matches_sn() {
        // conformal_square(re, 0, ks, m) = sn(re*ks, m) since sn_v=0
        let m = 0.5_f64;
        let ks = k_elliptic(m);
        let re = 0.5_f64;
        let (cs_re, cs_im) = conformal_square(re, 0.0, ks, m);
        let (sn, _cn, _dn) = sncndn(re * ks, m);
        assert!((cs_re - sn).abs() < 1e-13, "re_out={cs_re} sn={sn}");
        assert!(cs_im.abs() < 1e-14, "im_out should be 0: {cs_im}");
    }
}