gravitational 0.2.6

//! Black-hole ringdown: the quasinormal modes a merger's remnant rings in, from the fits of
//! [Berti, Cardoso & Will 2006](https://arxiv.org/abs/gr-qc/0512160), Tables VIII–X.
//!
//! A Kerr black hole of mass `M` and spin `j` rings in modes `(l, m, n)` at the dimensionless
//! frequency `F = Mω = f₁ + f₂(1 − j)^f₃` (with `G = c = 1`), damped with quality factor
//! `Q = πfτ = q₁ + q₂(1 − j)^q₃`. The fits cover `l = 2, 3, 4`, every `m` (positive co-rotating)
//! and the overtones `n = 0, 1, 2`, for `0 ≤ j ≤ 0.99`.
//!
//! Accuracy, from an exact Leaver solver on a 0.05 spin grid up to 0.99 (the frontier-gaps
//! research) and tested here against the paper's Tables II–IV at all eleven tabulated spins:
//! every frequency fit stays within the maximum error the paper states (at most 4%). So do the
//! quality factors, with two exceptions. Nine `l = 4`, `m = 0–2` fits exceed their printed
//! maxima by exactly 2.00 or 3.00 points (0.26 printed for 2.26, 0.12 for 3.11), a leading digit
//! lost from Table X's error column; their true maxima are 2.26–3.53%. Four fits genuinely
//! miss: `(2, −1, 1)` 3.22% against 0.67, `(3, 3, 0)` 0.59 against 0.42, `(4, 4, 1)` 1.02 against
//! 0.77 and `(4, 4, 2)` 1.09 against 0.10.
//!
//! [`qnm_exact`] solves for the same modes exactly, by Leaver's continued fraction started from
//! the fit, to about 10⁻¹⁴.
//!
//! ```
//! use gravitational::ringdown::{damping_time, frequency, qnm};
//! use metrology::SOLAR_MASS;
//!
//! // GW150914's remnant, 68 M☉ in the detector frame with spin 0.67: LVC predicted
//! // f = 251 ± 8 Hz and τ = 4.0 ± 0.3 ms for the (2, 2, 0) mode.
//! let (f, q) = qnm(2, 2, 0, 0.67).unwrap();
//! let hz = frequency(68.0 * SOLAR_MASS, f);
//! assert!((hz - 251.0).abs() < 8.0 && (damping_time(hz, q) * 1e3 - 4.0).abs() < 0.3);
//! ```

use core::f64::consts::PI;
use core::ops::{Add, Div, Mul, Neg, Sub};

use metrology::{C, G};

/// `((l, m, n), [f₁, f₂, f₃], [q₁, q₂, q₃], F error %, Q error %)`, the coefficients and the
/// maximum errors over `j ∈ [0, 0.99]` as printed in Tables VIII–X.
type Fit = ((u32, i32, u32), [f64; 3], [f64; 3], f64, f64);

#[rustfmt::skip]
const FITS: [Fit; 63] = [
    ((2, 2, 0), [1.5251, -1.1568, 0.1292], [0.7, 1.4187, -0.499], 1.85, 0.88),
    ((2, 2, 1), [1.3673, -1.026, 0.1628], [0.1, 0.5436, -0.4731], 1.56, 1.69),
    ((2, 2, 2), [1.3223, -1.0257, 0.186], [-0.1, 0.4206, -0.4256], 1.91, 2.52),
    ((2, 1, 0), [0.6, -0.2339, 0.4175], [-0.3, 2.3561, -0.2277], 2.03, 3.65),
    ((2, 1, 1), [0.58, -0.2416, 0.4708], [-0.33, 0.9501, -0.2072], 2.4, 3.18),
    ((2, 1, 2), [0.566, -0.274, 0.496], [-0.1, 0.4173, -0.2774], 4.04, 2.46),
    ((2, 0, 0), [0.4437, -0.0739, 0.335], [4.0, -1.955, 0.142], 1.04, 2.63),
    ((2, 0, 1), [0.4185, -0.0768, 0.4355], [1.25, -0.6359, 0.1614], 1.5, 4.01),
    ((2, 0, 2), [0.3734, -0.0794, 0.6306], [0.56, -0.2589, 0.3034], 2.72, 4.33),
    ((2, -1, 0), [0.3441, 0.0293, 2.001], [2.0, 0.1078, 5.0069], 0.07, 2.82),
    ((2, -1, 1), [0.3165, 0.0301, 2.3415], [0.61, 0.0276, 13.1683], 0.05, 0.67),
    ((2, -1, 2), [0.2696, 0.0315, 2.7755], [0.29, 0.0276, 6.4715], 0.43, 2.4),
    ((2, -2, 0), [0.2938, 0.0782, 1.3546], [1.67, 0.4192, 1.47], 0.63, 0.71),
    ((2, -2, 1), [0.2528, 0.0921, 1.3344], [0.455, 0.1729, 1.3617], 0.87, 0.79),
    ((2, -2, 2), [0.1873, 0.1117, 1.3322], [0.185, 0.1266, 1.3661], 1.34, 1.16),
    ((3, 3, 0), [1.8956, -1.3043, 0.1818], [0.9, 2.343, -0.481], 1.36, 0.42),
    ((3, 3, 1), [1.8566, -1.2818, 0.1934], [0.2274, 0.8173, -0.4731], 1.35, 0.88),
    ((3, 3, 2), [1.8004, -1.2558, 0.2133], [0.04, 0.5445, -0.4539], 1.28, 1.52),
    ((3, 2, 0), [1.1481, -0.5552, 0.3002], [0.8313, 2.3773, -0.3655], 1.09, 1.28),
    ((3, 2, 1), [1.1226, -0.5471, 0.3264], [0.23, 0.8025, -0.3684], 1.23, 0.51),
    ((3, 2, 2), [1.0989, -0.555, 0.3569], [0.1, 0.4804, -0.3784], 1.41, 0.81),
    ((3, 1, 0), [0.8345, -0.2405, 0.4095], [23.845, -20.724, 0.03837], 1.12, 3.47),
    ((3, 1, 1), [0.8105, -0.2342, 0.466], [8.853, -7.8506, 0.03418], 1.55, 3.64),
    ((3, 1, 2), [0.7684, -0.2252, 0.5805], [2.18, -1.6273, 0.1136], 2.67, 4.04),
    ((3, 0, 0), [0.6873, -0.09282, 0.3479], [6.7841, -3.6112, 0.0948], 0.83, 3.99),
    ((3, 0, 1), [0.6687, -0.09155, 0.4021], [2.0075, -0.993, 0.12297], 0.95, 4.18),
    ((3, 0, 2), [0.6343, -0.08915, 0.5117], [0.9, -0.3409, 0.2679], 1.28, 2.89),
    ((3, -1, 0), [0.5751, 0.02508, 3.136], [3.0464, 0.1162, -0.2812], 0.42, 2.65),
    ((3, -1, 1), [0.5584, 0.02514, 3.4154], [1.2, -0.1928, 0.1037], 0.42, 2.75),
    ((3, -1, 2), [0.5271, 0.02561, 3.8011], [1.0, -0.4424, 0.02467], 0.29, 3.15),
    ((3, -2, 0), [0.5158, 0.08195, 1.4084], [2.9, 0.3356, 2.305], 0.35, 0.72),
    ((3, -2, 1), [0.4951, 0.08577, 1.4269], [0.9, 0.1295, 1.6142], 0.41, 0.8),
    ((3, -2, 2), [0.4567, 0.093, 1.4469], [0.49, 0.0848, 1.9737], 0.53, 0.52),
    ((3, -3, 0), [0.4673, 0.1296, 1.3255], [2.55, 0.6576, 1.3378], 0.61, 0.79),
    ((3, -3, 1), [0.4413, 0.1387, 1.3178], [0.79, 0.2381, 1.3706], 0.68, 0.73),
    ((3, -3, 2), [0.3933, 0.1555, 1.3037], [0.407, 0.1637, 1.3819], 0.82, 0.88),
    ((4, 4, 0), [2.3, -1.5056, 0.2244], [1.1929, 3.1191, -0.4825], 1.83, 0.37),
    ((4, 4, 1), [2.3, -1.5173, 0.2271], [0.3, 1.1034, -0.4703], 1.75, 0.77),
    ((4, 4, 2), [2.3, -1.5397, 0.2321], [0.11, 0.6997, -0.4607], 1.61, 0.1),
    ((4, 3, 0), [1.6869, -0.8862, 0.2822], [1.4812, 2.8096, -0.4271], 1.05, 0.14),
    ((4, 3, 1), [1.6722, -0.8843, 0.2923], [0.4451, 0.9569, -0.425], 1.1, 0.37),
    ((4, 3, 2), [1.6526, -0.8888, 0.3081], [0.22, 0.5904, -0.4236], 1.15, 0.66),
    ((4, 2, 0), [1.2702, -0.4685, 0.3835], [-3.6, 7.7749, -0.1491], 1.11, 0.97),
    ((4, 2, 1), [1.2462, -0.458, 0.4139], [-1.5, 2.8601, -0.1392], 1.39, 0.12),
    ((4, 2, 2), [1.2025, -0.4401, 0.4769], [-1.5, 2.2784, -0.1124], 2.26, 0.31),
    ((4, 1, 0), [1.0507, -0.2478, 0.4348], [14.0, -9.824, 0.09047], 0.97, 0.81),
    ((4, 1, 1), [1.0337, -0.2439, 0.4695], [4.2, -2.8399, 0.1081], 1.15, 0.91),
    ((4, 1, 2), [1.0019, -0.2374, 0.5397], [2.2, -1.4195, 0.1372], 1.53, 0.53),
    ((4, 0, 0), [0.9175, -0.1144, 0.3511], [7.0, -2.7934, 0.1708], 0.75, 0.26),
    ((4, 0, 1), [0.9028, -0.1127, 0.3843], [2.2, -0.8308, 0.2023], 0.82, 0.26),
    ((4, 0, 2), [0.8751, -0.1096, 0.4516], [1.2, -0.4159, 0.2687], 0.96, 0.6),
    ((4, -1, 0), [0.7908, 0.02024, 5.4628], [4.6, -0.4038, 0.4629], 0.96, 2.52),
    ((4, -1, 1), [0.7785, 0.02005, 5.8547], [1.6, -0.2323, 0.2306], 0.98, 2.37),
    ((4, -1, 2), [0.7549, 0.01985, 6.5272], [1.6, -0.8136, 0.03163], 0.96, 2.32),
    ((4, -2, 0), [0.7294, 0.07842, 1.5646], [4.0, 0.2777, 2.0647], 0.23, 2.11),
    ((4, -2, 1), [0.7154, 0.07979, 1.5852], [1.32, 0.08694, 4.3255], 0.25, 0.75),
    ((4, -2, 2), [0.6885, 0.08259, 1.6136], [0.75, 0.05803, 3.7971], 0.32, 0.66),
    ((4, -3, 0), [0.6728, 0.1338, 1.3413], [3.7, 0.5829, 1.6681], 0.43, 0.45),
    ((4, -3, 1), [0.6562, 0.1377, 1.3456], [1.18, 0.2111, 1.4129], 0.46, 0.7),
    ((4, -3, 2), [0.6244, 0.1454, 1.3513], [0.66, 0.1385, 1.3742], 0.52, 0.82),
    ((4, -4, 0), [0.6256, 0.18, 1.3218], [3.4, 0.8696, 1.4074], 0.62, 0.63),
    ((4, -4, 1), [0.6061, 0.1869, 1.3168], [1.08, 0.3095, 1.3279], 0.67, 0.81),
    ((4, -4, 2), [0.5686, 0.2003, 1.3068], [0.598, 0.2015, 1.3765], 0.74, 0.69),
];

/// The dimensionless frequency `Mω` and quality factor `Q` of mode `(l, m, n)` of a Kerr black
/// hole with spin `j`, or `None` for a mode outside the fits.
pub fn qnm(l: u32, m: i32, n: u32, j: f64) -> Option<(f64, f64)> {
    let (_, f, q, _, _) = FITS.iter().find(|row| row.0 == (l, m, n))?;
    Some((f[0] + f[1] * (1.0 - j).powf(f[2]), q[0] + q[1] * (1.0 - j).powf(q[2])))
}

/// The exact `(F, Q)` of mode `(l, m, n)` of a Kerr black hole with spin `0 ≤ j < 1`, as [`qnm`]
/// returns them, or `None` for a mode outside the fits or a root that does not converge to it.
///
/// The method is that of [Cook & Zalutskiy 2014](https://arxiv.org/abs/1410.7698) (`G = c = M =
/// 1`, `s = −2`): Leaver's radial continued fraction in their confluent-Heun form, inverted `n`
/// times and closed with Nollert's tail through `u₃`; the angular separation constant from their
/// spectral (spherical-harmonic) matrix, re-solved at every trial `ω`; and a complex secant
/// started from the BCW06 fit, so it covers the fits' 63 modes (`n ≤ 2`). It reproduces Berti's
/// Schwarzschild overtones and the `qnm` package's `(2, 2, 0)` at `j = 0.68` to 6×10⁻¹⁶, and all
/// 693 entries of BCW06 Tables II–IV within 5.4×10⁻⁵. Up to `j = 0.99` a fraction four times
/// deeper moves no root by more than 10⁻¹⁴. Beyond the fits' `j = 0.99` it still converges (to
/// 10⁻¹² for the slowest modes at 0.9999), but returns `None` where the fit's guess is more than
/// 10% from the root, as for `(2, 2, n)` at `j = 0.9999`.
///
/// Pitfalls, from the frontier-gaps research: the tail takes the minimal branch, `Re u₁ < 0`;
/// for `m < 0` the root wanted is the counter-rotating one with `Re ω > 0`, which the fit's
/// guess selects and a final check enforces; overtones `n ≥ 5` approach the algebraically
/// special `−2i` and are not covered; the fraction's depth grows as `1/√(1 − j²)`, which matters
/// above `j = 0.99`; Leaver's own angular continued fraction picks the wrong `l` without inversion
/// (36 of 315 cases), which the spectral matrix avoids; and Leaver and Berti use `2M = 1`
/// (`ω` doubled, `a` halved), Cook, `qnm` and BCW06's tables `M = 1`.
pub fn qnm_exact(l: u32, m: i32, n: u32, j: f64) -> Option<(f64, f64)> {
    if !(0.0..1.0).contains(&j) {
        return None;
    }
    let (f, q) = qnm(l, m, n, j)?;
    let fit = Cx(f, -f / (2.0 * q));
    // The tail converges as exp(2 Re u₁ √N) with u₁ ∝ (1 − j²)^(1/4), so the depth grows as 1/√(1 − j²).
    let depth = (600.0 / (1.0 - j * j).sqrt()) as usize;
    // ℓ(ℓ + 1) − s(s + 1) at c = 0, then carried from one trial ω to the next.
    let mut sep = Cx((l * (l + 1)) as f64 - 2.0, 0.0);
    let mut cf = |w: Cx| {
        sep = separation(w * j, l, m, sep);
        leaver(w, j, m, n as usize, sep, depth)
    };
    let (mut w0, mut w1) = (fit, fit * 1.000001);
    let (mut c0, mut c1) = (cf(w0), cf(w1));
    for _ in 0..50 {
        let w2 = w1 - c1 * (w1 - w0) / (c1 - c0);
        (w0, c0, w1) = (w1, c1, w2);
        if !w1.abs().is_finite() {
            return None;
        }
        if (w1 - w0).abs() < 1e-13 * w1.abs() {
            return (w1.0 > 0.0 && (w1 - fit).abs() < 0.1 * fit.abs()).then(|| (w1.0, -w1.0 / (2.0 * w1.1)));
        }
        c1 = cf(w1);
    }
    None
}

/// A complex number, enough for the solver.
#[derive(Clone, Copy, Debug)]
struct Cx(f64, f64);

impl Cx {
    fn abs(self) -> f64 {
        self.0.hypot(self.1)
    }

    /// The principal square root.
    fn sqrt(self) -> Cx {
        let r = self.abs();
        Cx(((r + self.0) / 2.0).sqrt(), ((r - self.0) / 2.0).sqrt().copysign(self.1))
    }
}

impl From<f64> for Cx {
    fn from(x: f64) -> Cx {
        Cx(x, 0.0)
    }
}

impl<T: Into<Cx>> Add<T> for Cx {
    type Output = Cx;
    fn add(self, z: T) -> Cx {
        let z = z.into();
        Cx(self.0 + z.0, self.1 + z.1)
    }
}

impl<T: Into<Cx>> Sub<T> for Cx {
    type Output = Cx;
    fn sub(self, z: T) -> Cx {
        let z = z.into();
        Cx(self.0 - z.0, self.1 - z.1)
    }
}

impl<T: Into<Cx>> Mul<T> for Cx {
    type Output = Cx;
    fn mul(self, z: T) -> Cx {
        let z = z.into();
        Cx(self.0 * z.0 - self.1 * z.1, self.0 * z.1 + self.1 * z.0)
    }
}

impl<T: Into<Cx>> Div<T> for Cx {
    type Output = Cx;
    fn div(self, z: T) -> Cx {
        let z = z.into();
        let d = z.0 * z.0 + z.1 * z.1;
        Cx((self.0 * z.0 + self.1 * z.1) / d, (self.1 * z.0 - self.0 * z.1) / d)
    }
}

impl Neg for Cx {
    type Output = Cx;
    fn neg(self) -> Cx {
        Cx(-self.0, -self.1)
    }
}

/// Cook & Zalutskiy's `Cf(n; N)` (their Eqs. 34–42, 48), zero at a quasinormal frequency `w` of
/// a black hole of spin `a`, for azimuthal number `m`, inversion `n`, separation constant `sep`
/// and depth `depth`, with `s = −2` and the exponents `(ζ₊, ξ₋, η₊)`.
fn leaver(w: Cx, a: f64, m: i32, n: usize, sep: Cx, depth: usize) -> Cx {
    let (s, i) = (-2.0, Cx(0.0, 1.0));
    let root = (1.0 - a * a).sqrt();
    let (rp, rm) = (1.0 + root, 1.0 - root);
    let ma = m as f64 * a;
    let (sig_p, sig_m) = ((w * 2.0 * rp - ma) / (rp - rm), (w * 2.0 * rm - ma) / (rp - rm));
    let (zeta, xi, eta) = (i * w, -(i * sig_p) - s, -(i * sig_m));
    let p = zeta * root;
    // Cook's α has + s·iω/ζ, which is s for ζ = iω.
    let al = xi + eta - zeta * 2.0 + (1.0 + 2.0 * s);
    let (ga, de) = (eta * 2.0 + (1.0 + s), xi * 2.0 + (1.0 + s));
    let gd = ga + de;
    let sigma = sep + w * w * (a * a - 8.0) + p * (al * 2.0 + ga - de) + (-(gd / 2.0) + (1.0 + s)) * (gd / 2.0 + s);
    let (d0, d1, d2) = (de, p * 4.0 - al * 2.0 + ga - de - 2.0, al * 2.0 - ga + 2.0);
    let (d3, d4) = (al * (p * 4.0 - de) - sigma, al * (al - ga + 1.0));
    let alpha = |k: f64| (d0 + 1.0) * k + d0 + k * k;
    let beta = |k: f64| (d1 + 2.0) * k + d3 - 2.0 * k * k;
    let gamma = |k: f64| (d2 - 3.0) * k + d4 - d2 + 2.0 + k * k;
    // Nollert's tail r_N = 1 + u₁/√N + u₂/N + u₃/N^(3/2) on the minimal branch, Re u₁ < 0.
    let u1 = (-(p * 4.0)).sqrt();
    let u1 = if u1.0 > 0.0 { -u1 } else { u1 };
    let u2 = (p * 8.0 - al * 4.0 + ga * 2.0 + de * 2.0 + 3.0) / -4.0;
    let u3 = (p * 32.0 * (p * 2.0 - al * 4.0 + ga + de * 3.0 + 4.0) + gd * 4.0 * (gd - 2.0) + sigma * 16.0 + 3.0) / (u1 * 32.0);
    let big = depth as f64;
    let tail = u1 / big.sqrt() + u2 / big + u3 / big.powf(1.5) + 1.0;
    let mut e = beta(big) + alpha(big) * tail;
    for k in (n + 1..depth).rev() {
        let k = k as f64;
        e = beta(k) - alpha(k) * gamma(k + 1.0) / e;
    }
    let mut result = beta(n as f64) - alpha(n as f64) * gamma(n as f64 + 1.0) / e;
    if n > 0 {
        let mut d = beta(0.0);
        for k in 1..n {
            let k = k as f64;
            d = beta(k) - alpha(k - 1.0) * gamma(k) / d;
        }
        result = result - alpha(n as f64 - 1.0) * gamma(n as f64) / d;
    }
    result
}

/// The spin-weighted (`s = −2`) spheroidal separation constant `A_lm(c)` nearest `guess`: the
/// eigenvalue of Cook & Zalutskiy's spectral matrix (their Eqs. 51–55) on `ℓ = max(|m|, 2)` to
/// `l + 12`, by Rayleigh-quotient inverse iteration. The matrix is complex symmetric, so the
/// quotient is the bilinear `vᵀMv/vᵀv`.
fn separation(c: Cx, l: u32, m: i32, guess: Cx) -> Cx {
    let (s, mf) = (-2.0, m as f64);
    let lmin = m.unsigned_abs().max(2);
    let size = (l + 13 - lmin) as usize;
    let ell = |k: usize| (lmin as usize + k) as f64;
    // F_sℓm; G_sℓm = F_s(ℓ−1)m, so F alone fills the matrix. It vanishes at ℓ_min − 1.
    let f = |l: f64| (((l + 1.0).powi(2) - mf * mf) / ((2.0 * l + 3.0) * (2.0 * l + 1.0)) * ((l + 1.0).powi(2) - s * s) / (l + 1.0).powi(2)).max(0.0).sqrt();
    let h = |l: f64| -mf * s / (l * (l + 1.0));
    let mut mat = vec![vec![Cx(0.0, 0.0); size]; size];
    for k in 0..size {
        let l = ell(k);
        mat[k][k] = c * c * -(f(l).powi(2) + f(l - 1.0).powi(2) + h(l).powi(2)) + c * (2.0 * s * h(l)) + (l * (l + 1.0) - s * (s + 1.0));
        if k + 1 < size {
            mat[k][k + 1] = (c * (2.0 * s) - c * c * (h(l) + h(l + 1.0))) * f(l);
            mat[k + 1][k] = mat[k][k + 1];
        }
        if k + 2 < size {
            mat[k][k + 2] = c * c * -(f(l) * f(l + 1.0));
            mat[k + 2][k] = mat[k][k + 2];
        }
    }
    let dot = |u: &[Cx], v: &[Cx]| u.iter().zip(v).fold(Cx(0.0, 0.0), |sum, (&a, &b)| sum + a * b);
    let mut v = vec![Cx(0.0, 0.0); size];
    v[(l - lmin) as usize] = Cx(1.0, 0.0);
    let mut lambda = guess;
    for _ in 0..20 {
        let Some(w) = solve(&mat, lambda, &v) else { break };
        let scale = w.iter().fold(0.0f64, |big, z| big.max(z.abs()));
        v = w.iter().map(|&z| z / scale).collect();
        let mv: Vec<Cx> = mat.iter().map(|row| dot(row, &v)).collect();
        let next = dot(&v, &mv) / dot(&v, &v);
        let step = (next - lambda).abs();
        lambda = next;
        if step < 1e-15 * lambda.abs() {
            break;
        }
    }
    lambda
}

/// `(M − λI)⁻¹v` by Gaussian elimination with partial pivoting, or `None` when it is singular.
fn solve(mat: &[Vec<Cx>], lambda: Cx, v: &[Cx]) -> Option<Vec<Cx>> {
    let size = v.len();
    let mut a: Vec<Vec<Cx>> = mat.iter().enumerate().map(|(k, row)| {
        let mut row = row.clone();
        row[k] = row[k] - lambda;
        row.push(v[k]);
        row
    }).collect();
    for k in 0..size {
        let pivot = (k..size).max_by(|&x, &y| a[x][k].abs().total_cmp(&a[y][k].abs()))?;
        a.swap(k, pivot);
        if a[k][k].abs() == 0.0 {
            return None;
        }
        let (top, rest) = a.split_at_mut(k + 1);
        let pivot_row = &top[k];
        for row in rest {
            let factor = row[k] / pivot_row[k];
            for (x, &p) in row[k..].iter_mut().zip(&pivot_row[k..]) {
                *x = *x - factor * p;
            }
        }
    }
    let mut x = vec![Cx(0.0, 0.0); size];
    for k in (0..size).rev() {
        let tail = (k + 1..size).fold(a[k][size], |sum, col| sum - a[k][col] * x[col]);
        x[k] = tail / a[k][k];
    }
    Some(x)
}

/// The ringing frequency (Hz) of a black hole of mass `mass` (kg, in the detector frame for an
/// observed signal) for the dimensionless frequency `f = Mω`: `fc³/(2πGM)`.
pub fn frequency(mass: f64, f: f64) -> f64 {
    f * C.powi(3) / (2.0 * PI * G * mass)
}

/// The damping time (s) of a mode ringing at `frequency` (Hz) with quality factor `q`: `Q/(πf)`.
pub fn damping_time(frequency: f64, q: f64) -> f64 {
    q / (PI * frequency)
}

#[cfg(test)]
mod tests {
    use super::*;
    use metrology::SOLAR_MASS;

    /// Tables II–IV: one line per `(l, n, j)`, then `Mω_R,Mω_I` for `m = l` down to `−l`. Two
    /// misprints are corrected against the exact solver: Table II's `(2, −2, 0)` at `j = 0.6`
    /// (printed .0890, exact 0.0888917) and Table IV's `(4, 2, 0)` at `j = 0` (printed .0941,
    /// exact 0.0941640, as the other eight `m` print).
    const TABLES: &str = "\
2 0 0.0 0.3737,0.089 0.3737,0.089 0.3737,0.089 0.3737,0.089 0.3737,0.089
2 0 0.1 0.387,0.0887 0.3804,0.0888 0.374,0.0889 0.3678,0.089 0.3618,0.0891
2 0 0.2 0.4021,0.0883 0.3882,0.0885 0.3751,0.0887 0.3627,0.0889 0.3511,0.0892
2 0 0.3 0.4195,0.0877 0.3973,0.088 0.377,0.0884 0.3584,0.0888 0.3413,0.0892
2 0 0.4 0.4398,0.0869 0.408,0.0873 0.3797,0.0878 0.3546,0.0885 0.3325,0.0891
2 0 0.5 0.4641,0.0856 0.4206,0.0862 0.3833,0.0871 0.3515,0.0881 0.3243,0.089
2 0 0.6 0.494,0.0838 0.436,0.0846 0.3881,0.086 0.3489,0.0876 0.3168,0.0889
2 0 0.7 0.5326,0.0808 0.4551,0.0821 0.3941,0.0845 0.3469,0.0869 0.3098,0.0887
2 0 0.8 0.586,0.0756 0.4802,0.078 0.4019,0.0822 0.3454,0.086 0.3033,0.0885
2 0 0.9 0.6716,0.0649 0.5163,0.0698 0.412,0.0785 0.3444,0.0849 0.2972,0.0883
2 0 0.98 0.8254,0.0386 0.5642,0.0516 0.4223,0.0735 0.3439,0.0837 0.2927,0.0881
2 1 0.0 0.3467,0.2739 0.3467,0.2739 0.3467,0.2739 0.3467,0.2739 0.3467,0.2739
2 1 0.1 0.3619,0.2725 0.3545,0.2731 0.3472,0.2737 0.34,0.2744 0.333,0.275
2 1 0.2 0.379,0.2705 0.3635,0.2717 0.3486,0.273 0.3344,0.2744 0.3206,0.2759
2 1 0.3 0.3984,0.268 0.374,0.2698 0.3511,0.2718 0.3296,0.2741 0.3093,0.2765
2 1 0.4 0.4208,0.2647 0.3863,0.267 0.3547,0.27 0.3256,0.2734 0.2989,0.2769
2 1 0.5 0.4474,0.2602 0.4009,0.2631 0.3594,0.2674 0.3225,0.2723 0.2893,0.2772
2 1 0.6 0.4798,0.2538 0.4183,0.2575 0.3655,0.2638 0.3201,0.2708 0.2803,0.2773
2 1 0.7 0.5212,0.2442 0.4399,0.2492 0.3732,0.2585 0.3184,0.2686 0.272,0.2773
2 1 0.8 0.5779,0.2281 0.4676,0.2358 0.3826,0.2507 0.3173,0.2658 0.2643,0.2772
2 1 0.9 0.6677,0.1953 0.5059,0.2097 0.3935,0.2385 0.3167,0.262 0.257,0.277
2 1 0.98 0.8249,0.1159 0.5477,0.1509 0.4014,0.2231 0.3164,0.2581 0.2515,0.2768
2 2 0.0 0.3011,0.4783 0.3011,0.4783 0.3011,0.4783 0.3011,0.4783 0.3011,0.4783
2 2 0.1 0.3192,0.4735 0.3104,0.4756 0.3017,0.4778 0.2932,0.4801 0.2846,0.4825
2 2 0.2 0.3393,0.4679 0.3214,0.4719 0.3038,0.4764 0.2866,0.4811 0.2697,0.4862
2 2 0.3 0.3619,0.4613 0.3342,0.4671 0.3074,0.4739 0.2813,0.4814 0.2559,0.4895
2 2 0.4 0.3878,0.4533 0.3492,0.4607 0.3124,0.4701 0.2772,0.4808 0.2433,0.4925
2 2 0.5 0.4179,0.4433 0.3669,0.4522 0.319,0.4647 0.2741,0.4794 0.2316,0.4952
2 2 0.6 0.4542,0.4303 0.3878,0.4407 0.3273,0.4571 0.2721,0.4768 0.2207,0.4977
2 2 0.7 0.4999,0.4123 0.4133,0.4241 0.3374,0.4464 0.2709,0.4729 0.2107,0.4999
2 2 0.8 0.5622,0.3839 0.4451,0.3984 0.3488,0.4307 0.2703,0.4674 0.2013,0.5019
2 2 0.9 0.6598,0.3275 0.4867,0.3502 0.3591,0.4067 0.2697,0.46 0.1925,0.5038
2 2 0.98 0.8238,0.1933 0.5201,0.2331 0.3599,0.3808 0.2686,0.4527 0.1858,0.5051
3 0 0.0 0.5994,0.0927 0.5994,0.0927 0.5994,0.0927 0.5994,0.0927 0.5994,0.0927 0.5994,0.0927 0.5994,0.0927
3 0 0.1 0.6208,0.0924 0.6137,0.0925 0.6067,0.0926 0.5999,0.0926 0.5932,0.0927 0.5867,0.0928 0.5802,0.0928
3 0 0.2 0.6448,0.092 0.6297,0.0921 0.6153,0.0923 0.6014,0.0924 0.588,0.0926 0.5752,0.0927 0.5628,0.0929
3 0 0.3 0.6721,0.0913 0.648,0.0915 0.6252,0.0918 0.6038,0.0921 0.5837,0.0923 0.5647,0.0926 0.5469,0.0928
3 0 0.4 0.7037,0.0902 0.6689,0.0906 0.6369,0.0911 0.6074,0.0915 0.5802,0.092 0.5553,0.0924 0.5323,0.0927
3 0 0.5 0.7409,0.0887 0.6934,0.0893 0.6506,0.09 0.6121,0.0908 0.5776,0.0915 0.5467,0.092 0.5188,0.0925
3 0 0.6 0.7862,0.0864 0.7228,0.0873 0.667,0.0885 0.6183,0.0897 0.5758,0.0908 0.5388,0.0916 0.5063,0.0922
3 0 0.7 0.8437,0.0829 0.7592,0.0842 0.687,0.0861 0.6261,0.0882 0.5749,0.0899 0.5316,0.0912 0.4946,0.0919
3 0 0.8 0.9219,0.077 0.8068,0.0791 0.7121,0.0824 0.636,0.086 0.5749,0.0888 0.525,0.0906 0.4837,0.0916
3 0 0.9 1.0446,0.0655 0.8762,0.0689 0.7455,0.0758 0.6487,0.0826 0.5758,0.0872 0.5191,0.09 0.4735,0.0913
3 0 0.98 1.2602,0.0387 0.9769,0.0453 0.7833,0.0643 0.6615,0.0782 0.5773,0.0856 0.5146,0.0894 0.4657,0.091
3 1 0.0 0.5826,0.2813 0.5826,0.2813 0.5826,0.2813 0.5826,0.2813 0.5826,0.2813 0.5826,0.2813 0.5826,0.2813
3 1 0.1 0.6053,0.2802 0.5978,0.2805 0.5904,0.2808 0.5832,0.2811 0.5761,0.2814 0.5691,0.2817 0.5623,0.282
3 1 0.2 0.6306,0.2785 0.6148,0.2791 0.5995,0.2798 0.5848,0.2804 0.5706,0.2811 0.5569,0.2817 0.5437,0.2823
3 1 0.3 0.6593,0.2761 0.6341,0.2771 0.6102,0.2782 0.5876,0.2793 0.5662,0.2804 0.5459,0.2814 0.5268,0.2823
3 1 0.4 0.6924,0.2726 0.6563,0.2741 0.6227,0.2758 0.5916,0.2776 0.5627,0.2793 0.536,0.2808 0.5112,0.2822
3 1 0.5 0.7312,0.2676 0.6821,0.2698 0.6375,0.2724 0.5969,0.2751 0.5602,0.2777 0.527,0.28 0.4967,0.2818
3 1 0.6 0.7782,0.2604 0.713,0.2634 0.655,0.2674 0.6038,0.2717 0.5587,0.2756 0.5188,0.2789 0.4833,0.2813
3 1 0.7 0.8375,0.2494 0.751,0.2538 0.6762,0.26 0.6124,0.2668 0.558,0.2728 0.5113,0.2775 0.4708,0.2807
3 1 0.8 0.9176,0.2315 0.8004,0.2379 0.7025,0.2484 0.623,0.2598 0.5582,0.2692 0.5045,0.2758 0.459,0.28
3 1 0.9 1.0425,0.1966 0.8714,0.2068 0.7362,0.2277 0.6357,0.2491 0.5593,0.2644 0.4983,0.2739 0.448,0.2792
3 1 0.98 1.2599,0.1161 0.9706,0.1341 0.7679,0.193 0.6468,0.2359 0.5606,0.2593 0.4938,0.272 0.4396,0.2784
3 2 0.0 0.5517,0.4791 0.5517,0.4791 0.5517,0.4791 0.5517,0.4791 0.5517,0.4791 0.5517,0.4791 0.5517,0.4791
3 2 0.1 0.5766,0.4763 0.5684,0.4771 0.5603,0.4779 0.5523,0.4787 0.5445,0.4795 0.5368,0.4803 0.5292,0.4812
3 2 0.2 0.6043,0.4725 0.5872,0.4741 0.5705,0.4758 0.5543,0.4775 0.5386,0.4792 0.5234,0.4809 0.5086,0.4826
3 2 0.3 0.6356,0.4674 0.6084,0.4699 0.5825,0.4725 0.5577,0.4753 0.534,0.4781 0.5114,0.4809 0.4897,0.4836
3 2 0.4 0.6714,0.4605 0.6328,0.464 0.5966,0.4679 0.5625,0.4721 0.5305,0.4763 0.5005,0.4804 0.4723,0.4842
3 2 0.5 0.7131,0.4511 0.6612,0.4558 0.6131,0.4614 0.5689,0.4675 0.5282,0.4736 0.4907,0.4793 0.4561,0.4845
3 2 0.6 0.7632,0.4379 0.6948,0.4441 0.6327,0.4522 0.5769,0.4611 0.527,0.4698 0.4819,0.4778 0.441,0.4846
3 2 0.7 0.8258,0.4185 0.7358,0.4268 0.6561,0.4387 0.5868,0.4521 0.5267,0.4648 0.4739,0.4757 0.4269,0.4843
3 2 0.8 0.9096,0.3874 0.7884,0.3988 0.6844,0.4177 0.5984,0.4391 0.5272,0.4581 0.4667,0.473 0.4136,0.4839
3 2 0.9 1.0383,0.3282 0.8622,0.3449 0.7177,0.3806 0.6106,0.4197 0.5282,0.4492 0.46,0.4696 0.4011,0.4834
3 2 0.98 1.2592,0.1935 0.9605,0.2181 0.7349,0.3224 0.6176,0.3976 0.5286,0.4403 0.455,0.4665 0.3916,0.4828
4 0 0.0 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942 0.8092,0.0942
4 0 0.1 0.8387,0.094 0.8313,0.094 0.824,0.094 0.8168,0.0941 0.8098,0.0941 0.8028,0.0941 0.796,0.0942 0.7893,0.0942 0.7826,0.0942
4 0 0.2 0.8717,0.0935 0.856,0.0936 0.8407,0.0937 0.8259,0.0938 0.8116,0.0939 0.7977,0.094 0.7842,0.0941 0.7712,0.0941 0.7585,0.0942
4 0 0.3 0.9092,0.0929 0.8839,0.093 0.8597,0.0932 0.8366,0.0934 0.8146,0.0935 0.7937,0.0937 0.7737,0.0938 0.7546,0.094 0.7365,0.0941
4 0 0.4 0.9525,0.0918 0.9159,0.0921 0.8815,0.0924 0.8493,0.0927 0.8191,0.093 0.7908,0.0933 0.7643,0.0935 0.7395,0.0937 0.7162,0.0939
4 0 0.5 1.0034,0.0903 0.9532,0.0907 0.9069,0.0912 0.8642,0.0917 0.825,0.0922 0.789,0.0927 0.7559,0.0931 0.7255,0.0934 0.6974,0.0936
4 0 0.6 1.065,0.0879 0.9978,0.0885 0.9369,0.0893 0.882,0.0902 0.8326,0.0911 0.7883,0.0919 0.7485,0.0925 0.7125,0.093 0.68,0.0933
4 0 0.7 1.1427,0.0842 1.053,0.0852 0.9734,0.0865 0.9035,0.0881 0.8423,0.0896 0.7889,0.0909 0.742,0.0919 0.7005,0.0926 0.6637,0.093
4 0 0.8 1.2475,0.0781 1.1254,0.0796 1.0198,0.082 0.9301,0.0848 0.8546,0.0874 0.7907,0.0895 0.7363,0.0911 0.6894,0.0921 0.6485,0.0926
4 0 0.9 1.4104,0.0662 1.2322,0.0685 1.0836,0.0733 0.9647,0.0791 0.8702,0.0841 0.7941,0.0877 0.7314,0.0901 0.679,0.0915 0.6343,0.0922
4 0 0.98 1.6919,0.0388 1.3944,0.0425 1.1631,0.0558 1.0019,0.0703 0.886,0.08 0.7979,0.0858 0.7282,0.0892 0.6712,0.091 0.6235,0.0919
4 1 0.0 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843 0.7966,0.2843
4 1 0.1 0.8271,0.2835 0.8194,0.2837 0.8119,0.2838 0.8046,0.284 0.7973,0.2841 0.7901,0.2843 0.783,0.2844 0.7761,0.2846 0.7692,0.2847
4 1 0.2 0.8611,0.2821 0.8449,0.2824 0.8293,0.2828 0.814,0.2831 0.7992,0.2835 0.7848,0.2838 0.7709,0.2841 0.7574,0.2844 0.7442,0.2847
4 1 0.3 0.8997,0.2799 0.8738,0.2805 0.8489,0.2811 0.8252,0.2817 0.8025,0.2823 0.7808,0.2829 0.7601,0.2835 0.7403,0.284 0.7214,0.2845
4 1 0.4 0.9441,0.2766 0.9067,0.2775 0.8715,0.2785 0.8384,0.2796 0.8072,0.2806 0.778,0.2816 0.7505,0.2826 0.7247,0.2834 0.7003,0.284
4 1 0.5 0.9962,0.2717 0.9451,0.2731 0.8977,0.2747 0.8539,0.2765 0.8136,0.2782 0.7763,0.2799 0.742,0.2813 0.7103,0.2825 0.6809,0.2834
4 1 0.6 1.0591,0.2645 0.9909,0.2665 0.9287,0.2691 0.8725,0.272 0.8217,0.2748 0.7759,0.2775 0.7345,0.2797 0.6969,0.2814 0.6628,0.2826
4 1 0.7 1.1382,0.2532 1.0474,0.2562 0.9663,0.2604 0.8947,0.2653 0.8319,0.2701 0.7767,0.2743 0.7279,0.2777 0.6846,0.2801 0.6459,0.2817
4 1 0.8 1.2445,0.2347 1.1212,0.2391 1.0137,0.2464 0.9221,0.2551 0.8446,0.2634 0.7788,0.2702 0.7222,0.2752 0.6732,0.2787 0.63,0.2807
4 1 0.9 1.409,0.1986 1.2293,0.2055 1.0779,0.2199 0.9565,0.2376 0.8602,0.2532 0.7822,0.2646 0.7174,0.2723 0.6625,0.277 0.6152,0.2796
4 1 0.98 1.6917,0.1165 1.3917,0.1268 1.1519,0.1663 0.9899,0.2113 0.8746,0.2407 0.7857,0.2586 0.714,0.2694 0.6545,0.2756 0.6039,0.2787
4 2 0.0 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799 0.7727,0.4799
4 2 0.1 0.8049,0.478 0.7969,0.4784 0.789,0.4788 0.7812,0.4791 0.7734,0.4795 0.7658,0.4799 0.7583,0.4803 0.7509,0.4807 0.7436,0.4811
4 2 0.2 0.8409,0.4751 0.8239,0.4759 0.8074,0.4767 0.7913,0.4775 0.7756,0.4784 0.7604,0.4792 0.7455,0.48 0.7311,0.4808 0.717,0.4816
4 2 0.3 0.8815,0.4708 0.8544,0.4721 0.8283,0.4735 0.8033,0.4749 0.7794,0.4763 0.7563,0.4778 0.7342,0.4792 0.713,0.4805 0.6927,0.4817
4 2 0.4 0.928,0.4647 0.8892,0.4666 0.8524,0.4687 0.8176,0.471 0.7847,0.4733 0.7536,0.4756 0.7243,0.4777 0.6965,0.4797 0.6702,0.4814
4 2 0.5 0.9824,0.4559 0.9296,0.4586 0.8802,0.4618 0.8344,0.4654 0.7918,0.469 0.7522,0.4725 0.7155,0.4757 0.6813,0.4785 0.6493,0.4808
4 2 0.6 1.0478,0.4431 0.9775,0.4469 0.9131,0.4518 0.8543,0.4573 0.8008,0.463 0.7522,0.4683 0.7078,0.473 0.6672,0.4769 0.6299,0.48
4 2 0.7 1.1295,0.4237 1.0365,0.4291 0.9526,0.4367 0.878,0.4456 0.8119,0.4546 0.7534,0.4627 0.7012,0.4696 0.6542,0.475 0.6117,0.4789
4 2 0.8 1.2387,0.3921 1.1129,0.3998 1.002,0.4123 0.9065,0.4276 0.8253,0.4426 0.7558,0.4554 0.6954,0.4654 0.6422,0.4727 0.5946,0.4777
4 2 0.9 1.4061,0.3313 1.2236,0.3428 1.0667,0.3665 0.9402,0.3972 0.8404,0.4247 0.7591,0.4455 0.6904,0.4602 0.631,0.4701 0.5786,0.4764
4 2 0.98 1.6913,0.1942 1.3871,0.2096 1.1283,0.273 0.9655,0.3535 0.852,0.4037 0.7619,0.4352 0.6868,0.4552 0.6225,0.4678 0.5664,0.4752";

    /// The true maximum Q errors (%) where they exceed the printed ones, from the exact solver:
    /// nine `l = 4` rows that lost a leading digit, and four genuine misses. `(3, 3, 0)` and
    /// `(4, 4, 1)` peak between the tabulated spins, so the tables alone would not show them.
    const Q_TRUE: [((u32, i32, u32), f64); 13] = [
        ((4, 2, 0), 2.97), ((4, 2, 1), 3.11), ((4, 2, 2), 3.31), ((4, 1, 0), 2.81), ((4, 1, 1), 2.91),
        ((4, 1, 2), 3.53), ((4, 0, 0), 2.26), ((4, 0, 1), 2.26), ((4, 0, 2), 2.60),
        ((2, -1, 1), 3.22), ((3, 3, 0), 0.59), ((4, 4, 1), 1.02), ((4, 4, 2), 1.09),
    ];

    #[test]
    fn every_fit_against_the_tables() {
        let (mut checked, mut worst) = (0, 0.0f64);
        for line in TABLES.lines() {
            let mut words = line.split_whitespace();
            let (l, n, j): (u32, u32, f64) = (words.next().unwrap().parse().unwrap(), words.next().unwrap().parse().unwrap(), words.next().unwrap().parse().unwrap());
            for (k, pair) in words.enumerate() {
                let m = l as i32 - k as i32;
                let (wr, wi) = pair.split_once(',').unwrap();
                let (wr, wi): (f64, f64) = (wr.parse().unwrap(), wi.parse().unwrap());
                let (f, q) = qnm(l, m, n, j).unwrap();
                let (_, _, _, f_error, q_error) = FITS.iter().find(|row| row.0 == (l, m, n)).unwrap();
                // The tables round to four decimals, which allows 5×10⁻⁵ in each part.
                let (f_round, q_round) = (5e-5 / wr, 5e-5 / wr + 5e-5 / wi);
                assert!((f / wr - 1.0).abs() <= f_error / 100.0 + f_round, "F({l},{m},{n}) at j = {j}");
                let q_tol = Q_TRUE.iter().find(|row| row.0 == (l, m, n)).map_or(*q_error, |row| row.1) / 100.0;
                assert!((q / (wr / (2.0 * wi)) - 1.0).abs() <= q_tol + q_round, "Q({l},{m},{n}) at j = {j}");
                // The exact solver against the table, with its two misprints corrected above: the worst
                // is 5.4×10⁻⁵. Five l = 3, 4, n = 0 imaginary parts are rounded up from just past the
                // half (0.0790459 printed .0791), so BCW06's own values carry about 5×10⁻⁶.
                let (f, q) = qnm_exact(l, m, n, j).unwrap();
                worst = worst.max((f - wr).abs()).max((f / (2.0 * q) - wi).abs());
                checked += 1;
            }
        }
        assert_eq!(checked, 693);
        assert!(worst <= 6e-5, "{worst}");
        // The lost digit: (4, 2, 0) at j = 0 is 2.8% off, against a printed maximum of 0.97%
        // and a true one of 2.97%.
        let (_, q) = qnm(4, 2, 0, 0.0).unwrap();
        let off = (q / (0.8092 / (2.0 * 0.0942)) - 1.0).abs();
        assert!(off > 0.0097 + 0.01 && off < 0.0297, "{off}");
        assert!(qnm(5, 5, 0, 0.5).is_none());
    }

    #[test]
    fn exact_modes() {
        // Berti's s2l2.dat, Schwarzschild l = 2, printed in units 2M = 1 (halved here), for the
        // overtones the fits seed, n = 0–2.
        let berti = [
            (0.7473433688360838, 0.1779246313778714),
            (0.6934219937583268, 0.5478297505824697),
            (0.6021069092247328, 0.9565539664461437),
        ];
        let close = |(f, q): (f64, f64), (re, im): (f64, f64), tol: f64| (f - re).abs() < tol && (f / (2.0 * q) - im).abs() < tol;
        for (n, (re, im)) in berti.into_iter().enumerate() {
            for m in [2, -2] {
                let mode = qnm_exact(2, m, n as u32, 0.0).unwrap();
                assert!(close(mode, (re / 2.0, im / 2.0), 1e-15), "n = {n}, m = {m}: {mode:?}");
            }
        }
        // The qnm package's README: (2, 2, 0) at a = 0.68.
        let mode = qnm_exact(2, 2, 0, 0.68).unwrap();
        assert!(close(mode, (0.5239751042900845, 0.08151262363119974), 1e-15), "{mode:?}");
        // Near extremality the solver still converges; outside the fits or spins it declines.
        assert!(qnm_exact(2, 2, 0, 0.999).is_some());
        assert!(qnm_exact(5, 5, 0, 0.5).is_none() && qnm_exact(2, 2, 0, 1.0).is_none() && qnm_exact(2, 2, 0, -0.1).is_none());
    }

    #[test]
    fn gw150914() {
        // 68 M☉ and j = 0.67 (LVC, Tests of GR with GW150914, 2016): 248.4 Hz and 4.06 ms from
        // the central values, inside the predicted 251 ± 8 Hz and 4.0 ± 0.3 ms.
        let (f, q) = qnm(2, 2, 0, 0.67).unwrap();
        let hz = frequency(68.0 * SOLAR_MASS, f);
        assert!((hz / 248.4 - 1.0).abs() < 1e-3 && (damping_time(hz, q) / 4.06e-3 - 1.0).abs() < 2e-3, "{hz} {}", damping_time(hz, q));
        // At j = 0 the fit gives 0.3683 against Table II's 0.3737, inside its stated 1.85%.
        assert!((qnm(2, 2, 0, 0.0).unwrap().0 - 0.3683).abs() < 1e-4);
    }
}