/*
 * // Copyright (c) Radzivon Bartoshyk 3/2026. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

use crate::zpk::{bilinear_zpk, lp2bp_zpk, lp2hp_zpk, lp2lp_zpk, zpk2sos, zpk2tf};
use crate::{BaFilter, FilterType, ImpliirError, SosFilter, ZpkFilter};
use num_complex::Complex64;
use pxfm::{f_asin, f_cosh, f_exp, f_exp10m1, f_pow, f_sin, f_sincos, f_tanh, f_tanpi};
use std::f64::consts::PI;

// ─── builder ─────────────────────────────────────────────────────────────────

/// Builder for an elliptic (Cauer) IIR filter, matching `scipy.signal.ellip`.
///
/// Elliptic filters have equiripple in both the passband and stopband,
/// giving the sharpest possible rolloff for a given filter order and ripple
/// specification — at the cost of nonlinear phase response.
///
/// - `rp` — maximum passband ripple in dB (> 0)
/// - `rs` — minimum stopband attenuation in dB (> `rp`)
///
/// Three output formats are available:
/// - [`.sos()`][EllipticFilter::sos] — second-order sections (recommended)
/// - [`.ba()`][EllipticFilter::ba]   — transfer function `(b, a)` coefficients
/// - [`.zpk()`][EllipticFilter::zpk] — zeros, poles, gain
#[derive(Clone)]
pub struct EllipticFilter {
    order: usize,
    /// Maximum passband ripple in dB.
    rp: f64,
    /// Minimum stopband attenuation in dB.
    rs: f64,
    filter_type: FilterType,
    freqs: [f64; 2],
    /// Sample rate in Hz; defaults to 2.0 (normalised frequencies).
    fs: f64,
}

impl EllipticFilter {
    /// Create a new builder.
    ///
    /// - `order` — filter order (≥ 1)
    /// - `rp`    — maximum passband ripple in dB (> 0)
    /// - `rs`    — minimum stopband attenuation in dB (> rp)
    pub fn new(order: usize, rp: f64, rs: f64) -> Self {
        Self {
            order,
            rp,
            rs,
            filter_type: FilterType::Lowpass,
            freqs: [0.0; 2],
            fs: 2.0,
        }
    }

    pub fn lowpass(mut self, cutoff: f64) -> Self {
        self.filter_type = FilterType::Lowpass;
        self.freqs[0] = cutoff;
        self
    }

    pub fn highpass(mut self, cutoff: f64) -> Self {
        self.filter_type = FilterType::Highpass;
        self.freqs[0] = cutoff;
        self
    }

    pub fn bandpass(mut self, low: f64, high: f64) -> Self {
        self.filter_type = FilterType::Bandpass;
        self.freqs[0] = low;
        self.freqs[1] = high;
        self
    }

    /// Set the sample rate in Hz (default: `2.0`).
    pub fn sample_rate(mut self, fs: f64) -> Self {
        self.fs = fs;
        self
    }

    /// Design the filter and return second-order sections.
    pub fn sos(self) -> Result<Vec<SosFilter>, ImpliirError> {
        let zpk = elliptic_zpk(
            self.order,
            self.rp,
            self.rs,
            self.filter_type,
            self.freqs,
            self.fs,
        )?;
        Ok(zpk2sos(zpk))
    }

    /// Design the filter and return transfer function coefficients `(b, a)`.
    ///
    /// For orders above ~4 prefer [`.sos()`][EllipticFilter::sos].
    pub fn ba(self) -> Result<BaFilter, ImpliirError> {
        let zpk = elliptic_zpk(
            self.order,
            self.rp,
            self.rs,
            self.filter_type,
            self.freqs,
            self.fs,
        )?;
        Ok(zpk2tf(&zpk))
    }

    /// Design the filter and return zeros, poles, and gain.
    pub fn zpk(self) -> Result<ZpkFilter, ImpliirError> {
        elliptic_zpk(
            self.order,
            self.rp,
            self.rs,
            self.filter_type,
            self.freqs,
            self.fs,
        )
    }
}

// ─── pipeline ────────────────────────────────────────────────────────────────

fn elliptic_zpk(
    order: usize,
    rp: f64,
    rs: f64,
    ftype: FilterType,
    freqs: [f64; 2],
    fs: f64,
) -> Result<ZpkFilter, ImpliirError> {
    if order == 0 {
        return Err(ImpliirError::InvalidOrder);
    }
    if rp <= 0.0 {
        return Err(ImpliirError::InvalidRipple);
    }
    if rs <= rp {
        return Err(ImpliirError::InvalidAttenuation);
    }

    let nyquist = fs / 2.0;

    let warped: Vec<f64> = match ftype {
        FilterType::Lowpass | FilterType::Highpass => {
            let f = freqs[0];
            if f <= 0.0 || f >= nyquist {
                return Err(ImpliirError::InvalidFrequency);
            }
            let wn = f / nyquist;
            vec![4.0 * f_tanpi(wn / 2.0)]
        }
        FilterType::Bandpass => {
            let (lo, hi) = (freqs[0], freqs[1]);
            if lo <= 0.0 || hi >= nyquist || lo >= hi {
                return Err(ImpliirError::InvalidBandFrequencies);
            }
            let (wn_lo, wn_hi) = (lo / nyquist, hi / nyquist);
            vec![4.0 * f_tanpi(wn_lo / 2.0), 4.0 * f_tanpi(wn_hi / 2.0)]
        }
    };

    let (z, p, k) = ellipap(order, rp, rs)?;

    let (z, p, k) = match ftype {
        FilterType::Lowpass => lp2lp_zpk(z, p, k, warped[0]),
        FilterType::Highpass => lp2hp_zpk(z, p, k, warped[0]),
        FilterType::Bandpass => {
            let bw = warped[1] - warped[0];
            let wo = (warped[0] * warped[1]).sqrt();
            lp2bp_zpk(z, p, k, wo, bw)
        }
    };

    let (z, p, k) = bilinear_zpk(z, p, k, 2.0);

    Ok(ZpkFilter { z, p, k })
}

/// Analog elliptic lowpass prototype.
///
/// Returns `(zeros, poles, gain)` for an Nth-order elliptic filter with
/// `rp` dB passband ripple and `rs` dB stopband attenuation.
///
/// Uses the AGM (arithmetic-geometric mean) method for complete elliptic
/// integrals and the Landen ascending transformation for Jacobi elliptic
/// functions.
fn ellipap(
    order: usize,
    rp: f64,
    rs: f64,
) -> Result<(Vec<Complex64>, Vec<Complex64>, f64), ImpliirError> {
    // order=1: closed form
    if order == 1 {
        let p = -1.0 / f_exp10m1(0.1 * rp).sqrt();
        return Ok((vec![], vec![Complex64::new(p, 0.0)], -p));
    }

    let eps_sq = f_exp10m1(0.1 * rp);
    let eps = eps_sq.sqrt();
    let ck1_sq = eps_sq / f_exp10m1(0.1 * rs);
    if ck1_sq == 0.0 {
        return Err(ImpliirError::InvalidRpRsCombination);
    }

    let n = order as f64;
    let m = ellipdeg(order, ck1_sq);
    let capk = ellipk(m);

    // scipy: j = arange(1 - N%2, N, 2)
    // Same j array is used for BOTH zeros and poles — one ellipj call per entry.
    // For even N: j = [1, 3, ..., N-1]
    // For odd  N: j = [0, 2, ..., N-1]  (j=0 gives sn=0, filtered out for zeros)
    let j_vals: Vec<f64> = {
        let start = 1 - (order % 2); // 0 for odd, 1 for even
        (start..order).step_by(2).map(|v| v as f64).collect()
    };

    // Evaluate ellipj once per j value; reuse (s, c, d) for both zeros and poles.
    let scj: Vec<(f64, f64, f64)> = j_vals
        .iter()
        .map(|&jv| {
            let (s, c, d, _) = ellipj(jv * capk / n, m);
            (s, c, d)
        })
        .collect();

    // Zeros: z = i / (sqrt(m) * sn(j*K/N, m)), compress sn≈0 (happens at j=0)
    let zeros: Vec<Complex64> = scj
        .iter()
        .filter(|(s, _, _)| s.abs() > f64::EPSILON)
        .flat_map(|(s, _, _)| {
            let z = Complex64::new(0.0, 1.0 / (m.sqrt() * s));
            [z, z.conj()]
        })
        .collect();

    // v0: shift parameter that optimally places the poles
    let k1 = ellipk(ck1_sq);
    let r = arc_jac_sc1(1.0 / eps, ck1_sq);
    let v0 = capk * r / (n * k1);
    let (sv, cv, dv, _) = ellipj(v0, 1.0 - m);

    // Poles: p = -(c*d*sv*cv + i*s*dv) / (1 - (d*sv)^2)
    // Uses the same (s, c, d) from scj — NOT a separate j_pole range.
    let mut poles: Vec<Complex64> = scj
        .iter()
        .map(|(s, c, d)| {
            let num = Complex64::new(c * d * sv * cv, s * dv);
            let den = 1.0 - (d * sv).powi(2);
            -num / den
        })
        .collect();

    // Duplicate conjugates.
    // For odd order: the j=0 entry gives a real pole; only conjugate the complex ones.
    // For even order: all entries are complex; conjugate all.
    if order % 2 == 1 {
        let norm_sq: f64 = poles.iter().map(|p| p.norm_sqr()).sum::<f64>().sqrt();
        let threshold = f64::EPSILON * norm_sq;
        let extra: Vec<Complex64> = poles
            .iter()
            .filter(|p| p.im.abs() > threshold)
            .map(|p| p.conj())
            .collect();
        poles.extend(extra);
    } else {
        let extra: Vec<Complex64> = poles.iter().map(|p| p.conj()).collect();
        poles.extend(extra);
    }

    // Gain
    let prod_neg_p: Complex64 = poles
        .iter()
        .fold(Complex64::from(1.0), |acc, &p| acc * (-p));
    let prod_neg_z: Complex64 = zeros
        .iter()
        .fold(Complex64::from(1.0), |acc, &z| acc * (-z));
    let mut k = (prod_neg_p / prod_neg_z).re;
    if order.is_multiple_of(2) {
        k /= (1.0 + eps_sq).sqrt();
    }

    Ok((zeros, poles, k))
}

/// Complete elliptic integral of the first kind K(m) via AGM.
///
/// `m` is the parameter (not the modulus k); K(m) = K(k²).
/// Valid for `m` in `[0, 1)`.
fn ellipk(m: f64) -> f64 {
    if m == 0.0 {
        return PI / 2.0;
    }
    if m >= 1.0 {
        return f64::INFINITY;
    }
    let mut a = 1.0_f64;
    let mut b = (1.0 - m).sqrt();
    for _ in 0..100 {
        let a_new = (a + b) / 2.0;
        let b_new = (a * b).sqrt();
        // Stop when converged or when floats no longer change (stall detection)
        if (a_new - b_new).abs() < 1e-16 * a_new || a_new == a || b_new == b {
            return PI / (2.0 * a_new);
        }
        a = a_new;
        b = b_new;
    }
    PI / (2.0 * a)
}

/// Jacobi elliptic functions `(sn, cn, dn, phi)` for real `u` and parameter `m`.
///
/// Uses the AGM descending Landen transformation.
fn ellipj(u: f64, m: f64) -> (f64, f64, f64, f64) {
    if m == 0.0 {
        let (us, uc) = f_sincos(u);
        return (us, uc, 1.0, u);
    }
    if m == 1.0 {
        let t = f_tanh(u);
        let s = 1.0 / f_cosh(u);
        return (t, s, s, u);
    }

    // Build AGM sequence
    let mut a_seq = vec![1.0_f64];
    let mut b_seq = vec![(1.0 - m).sqrt()];
    let mut c_seq = vec![m.sqrt()];
    while c_seq.last().unwrap().abs() > 1e-16 * a_seq.last().unwrap() {
        let a_last = *a_seq.last().unwrap();
        let b_last = *b_seq.last().unwrap();
        a_seq.push((a_last + b_last) / 2.0);
        b_seq.push((a_last * b_last).sqrt());
        c_seq.push((a_last - b_last) / 2.0);
        if a_seq.len() > 100 {
            break;
        }
    }

    let n = a_seq.len() - 1;
    let mut phi = (2_f64.powi(n as i32)) * a_seq[n] * u;
    for i in (1..=n).rev() {
        let sin_phi = (c_seq[i] / a_seq[i] * f_sin(phi)).clamp(-1.0, 1.0);
        phi = (f_asin(sin_phi) + phi) / 2.0;
    }

    let (sn, cn) = f_sincos(phi);
    let dn = (1.0 - m * sn * sn).max(0.0).sqrt();
    (sn, cn, dn, phi)
}

/// Solve the degree equation `N * K(m) / K'(m) = K1(m1) / K1'(m1)` for `m`
/// using nome-based iteration (Orfanidis Eq. 49).
fn ellipdeg(order: usize, m1: f64) -> f64 {
    const MMAX: usize = 5;
    let k1 = ellipk(m1);
    let k1p = ellipk(1.0 - m1);
    let q1 = f_exp(-PI * k1p / k1);
    let q = f_pow(q1, 1.0 / order as f64);

    let num: f64 = (0..=MMAX).map(|i| q.powi((i * (i + 1)) as i32)).sum();
    let den: f64 = 1.0 + 2.0 * (1..=MMAX + 1).map(|i| q.powi((i * i) as i32)).sum::<f64>();

    16.0 * q * (num / den).powi(4)
}

/// Inverse Jacobi sn via ascending Landen transformation (Orfanidis Eq. 56).
/// Handles complex `w` — needed for `arc_jac_sc1`.
fn arc_jac_sn(w: Complex64, m: f64) -> Complex64 {
    let k = m.sqrt();
    if k > 1.0 {
        return Complex64::new(f64::NAN, 0.0);
    }
    if k == 1.0 {
        return w.atanh();
    }

    // Build Landen sequence
    let mut ks = vec![k];
    for _ in 0..20 {
        let k_last = *ks.last().unwrap();
        let kp = ((1.0 - k_last) * (1.0 + k_last)).sqrt();
        let k_next = (1.0 - kp) / (1.0 + kp);
        ks.push(k_next);
        if k_next.abs() < 1e-16 {
            break;
        }
    }

    let k_prod: f64 = ks[1..].iter().map(|&ki| 1.0 + ki).product();
    let cap_k = k_prod * PI / 2.0;

    // Ascending iteration on w
    let mut wn = w;
    for (&kn, &knext) in ks[..ks.len() - 1].iter().zip(ks[1..].iter()) {
        let denom = (Complex64::from(1.0 + knext))
            * (Complex64::from(1.0)
                + (Complex64::from(1.0) - Complex64::from(kn) * kn * wn * wn).sqrt());
        wn = Complex64::from(2.0) * wn / denom;
    }

    let u = Complex64::from(2.0 / PI) * wn.asin();
    cap_k * u
}

/// Real inverse Jacobi sc with complementary modulus.
/// Solves `w = sc(z, 1-m)` for `z`. Uses `sc(z,m) = -i * sn(i*z, 1-m)`.
fn arc_jac_sc1(w: f64, m: f64) -> f64 {
    let z_complex = arc_jac_sn(Complex64::new(0.0, w), m);
    z_complex.im
}

// ─── tests ───────────────────────────────────────────────────────────────────

#[cfg(test)]
mod elliptic_tests {
    use super::*;
    use crate::SosFilter;
    use crate::zpk::zpk2tf;

    const ATOL: f64 = 1e-7;

    fn assert_sos_close(got: &[SosFilter], expected: &[[f64; 6]], label: &str) {
        assert_eq!(
            got.len(),
            expected.len(),
            "{label}: got {} sections, expected {}",
            got.len(),
            expected.len()
        );
        for (i, (s, e)) in got.iter().zip(expected).enumerate() {
            for (j, (&g, &ex)) in s.b.iter().chain(s.a.iter()).zip(e).enumerate() {
                let err = (g - ex).abs();
                assert!(
                    err <= ATOL + ATOL * ex.abs(),
                    "{label} section {i} coeff {j}: got {g:.10e}, expected {ex:.10e}, err={err:.3e}"
                );
            }
        }
    }

    fn assert_vec_close(got: &[f64], expected: &[f64], label: &str) {
        assert_eq!(got.len(), expected.len(), "{label}: length mismatch");
        for (i, (&g, &e)) in got.iter().zip(expected).enumerate() {
            let err = (g - e).abs();
            assert!(
                err <= ATOL + ATOL * e.abs(),
                "{label}[{i}]: got {g:.10e}, expected {e:.10e}, err={err:.3e}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Reproduce in Python:
    //   from scipy.signal import ellip
    //   sos = ellip(N, rp, rs, freq, btype='...', fs=fs, output='sos')
    //   b, a = ellip(N, rp, rs, freq, btype='...', fs=fs, output='ba')
    //   z, p, k = ellip(N, rp, rs, freq, btype='...', fs=fs, output='zpk')
    // -----------------------------------------------------------------------

    /// ellip(2, 1.0, 40.0, 100.0, btype='low', fs=1000.0)
    #[test]
    fn test_lowpass_order2() {
        let f = EllipticFilter::new(2, 1.0, 40.0)
            .lowpass(100.0)
            .sample_rate(1000.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [[
            7.76415285226175339e-02,
            1.28102389780719333e-01,
            7.76415285226175200e-02,
            1.00000000000000000e+00,
            -1.20005151271612331e+00,
            5.18015213735451896e-01,
        ]];
        let expected_b = vec![
            7.76415285226175339e-02,
            1.28102389780719333e-01,
            7.76415285226175200e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.20005151271612331e+00,
            5.18015213735451896e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip lowpass n=2 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip lowpass n=2 b");
        assert_vec_close(&ba.a, &expected_a, "ellip lowpass n=2 a");
        assert_eq!(zpk.z.len(), 2);
        assert_eq!(zpk.p.len(), 2);
        assert!(
            (zpk.k - 7.76415285226175339e-02_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }

    /// ellip(3, 3.0, 60.0, 200.0, btype='high', fs=1000.0)
    #[test]
    fn test_highpass_order3() {
        let f = EllipticFilter::new(3, 3.0, 60.0)
            .highpass(200.0)
            .sample_rate(1000.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [
            [
                1.61016948254104786e-01,
                -3.06624681795000464e-01,
                1.61016948254104786e-01,
                1.00000000000000000e+00,
                4.10289446033275829e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                -1.00000000000000000e+00,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                -3.99768661152529736e-01,
                7.32323342541485234e-01,
            ],
        ];
        let expected_b = vec![
            1.61016948254104786e-01,
            -4.67641630049105306e-01,
            4.67641630049105306e-01,
            -1.61016948254104786e-01,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            1.05207848807460924e-02,
            5.68302480015749456e-01,
            3.00464538528582903e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip highpass n=3 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip highpass n=3 b");
        assert_vec_close(&ba.a, &expected_a, "ellip highpass n=3 a");
        assert_eq!(zpk.z.len(), 3);
        assert_eq!(zpk.p.len(), 3);
        assert!(
            (zpk.k - 1.61016948254104786e-01_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }

    /// ellip(4, 0.5, 80.0, [100.0, 300.0], btype='band', fs=1000.0)
    #[test]
    fn test_bandpass_order4() {
        let f = EllipticFilter::new(4, 0.5, 80.0)
            .bandpass(100.0, 300.0)
            .sample_rate(1000.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [
            [
                3.38877107791935375e-02,
                6.69392465351462823e-02,
                3.38877107791935445e-02,
                1.00000000000000000e+00,
                -8.34433384262342759e-02,
                5.20221206762832677e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.99504055026488092e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -1.02078622222406690e+00,
                6.07500461344596232e-01,
            ],
            [
                1.00000000000000000e+00,
                1.86594297631455741e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                5.84987724857857105e-01,
                8.19389347428538972e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.97244993943867031e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.53405601755750132e+00,
                8.85137193302973913e-01,
            ],
        ];
        let expected_b = vec![
            3.38877107791935375e-02,
            -4.27738778731329324e-03,
            -1.22647521321135267e-01,
            3.85043718511953294e-03,
            1.78444680721969950e-01,
            3.85043718511953901e-03,
            -1.22647521321135294e-01,
            -4.27738778731331752e-03,
            3.38877107791935583e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -2.05329785334994552e+00,
            3.06801134340021342e+00,
            -3.36329434563726304e+00,
            3.38860408701129945e+00,
            -2.46689876581910017e+00,
            1.56476989132644384e+00,
            -6.55521079070339296e-01,
            2.29211059128408584e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip bandpass n=4 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip bandpass n=4 b");
        assert_vec_close(&ba.a, &expected_a, "ellip bandpass n=4 a");
        assert_eq!(zpk.z.len(), 8);
        assert_eq!(zpk.p.len(), 8);
        assert!(
            (zpk.k - 3.38877107791935375e-02_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }

    /// ellip(2, 1.0, 40.0, 0.1, btype='low', fs=2.0) — normalised
    #[test]
    fn test_normalised_lp_order2() {
        let f = EllipticFilter::new(2, 1.0, 40.0)
            .lowpass(0.1)
            .sample_rate(2.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [[
            2.89617809426324056e-02,
            2.46089541426474274e-02,
            2.89617809426324126e-02,
            1.00000000000000000e+00,
            -1.61942888409922547e+00,
            7.12031890162515668e-01,
        ]];
        let expected_b = vec![
            2.89617809426324056e-02,
            2.46089541426474274e-02,
            2.89617809426324126e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.61942888409922547e+00,
            7.12031890162515668e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip normalised lowpass n=2 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip normalised lowpass n=2 b");
        assert_vec_close(&ba.a, &expected_a, "ellip normalised lowpass n=2 a");
        assert_eq!(zpk.z.len(), 2);
        assert_eq!(zpk.p.len(), 2);
        assert!(
            (zpk.k - 2.89617809426324056e-02_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }

    /// ellip(3, 0.5, 60.0, 0.2, btype='low', fs=2.0) — normalised, odd order
    #[test]
    fn test_normalised_lp_order3() {
        let f = EllipticFilter::new(3, 0.5, 60.0)
            .lowpass(0.2)
            .sample_rate(2.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [
            [
                1.90435342703293466e-02,
                2.42938495504018415e-02,
                1.90435342703293431e-02,
                1.00000000000000000e+00,
                -6.58662344673364819e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.33124084937645559e+00,
                6.96749413214023172e-01,
            ],
        ];
        let expected_b = vec![
            1.90435342703293466e-02,
            4.33373838207311846e-02,
            4.33373838207311846e-02,
            1.90435342703293431e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.98990319404982041e+00,
            1.57358763238928123e+00,
            -4.58922602157339599e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip normalised lowpass n=3 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip normalised lowpass n=3 b");
        assert_vec_close(&ba.a, &expected_a, "ellip normalised lowpass n=3 a");
        assert_eq!(zpk.z.len(), 3);
        assert_eq!(zpk.p.len(), 3);
        assert!(
            (zpk.k - 1.90435342703293466e-02_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }

    // ── ZPK round-trip ───────────────────────────────────────────────────────

    #[test]
    fn test_zpk_ba_roundtrip() {
        let ba = EllipticFilter::new(3, 3.0, 60.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .ba()
            .unwrap();
        let zpk = EllipticFilter::new(3, 3.0, 60.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        let ba_rt = zpk2tf(&zpk);
        assert_vec_close(&ba_rt.b, &ba.b, "elliptic roundtrip b");
        assert_vec_close(&ba_rt.a, &ba.a, "elliptic roundtrip a");
    }

    /// Elliptic filters have finite zeros (unlike Butterworth/Chebyshev I)
    #[test]
    fn test_has_finite_zeros() {
        let zpk = EllipticFilter::new(4, 1.0, 40.0)
            .lowpass(100.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        // Order-4 lowpass: 4 zeros, none should be exactly at -1
        assert_eq!(zpk.z.len(), 4);
        let all_at_minus_one = zpk
            .z
            .iter()
            .all(|z| (z.re + 1.0).abs() < 1e-6 && z.im.abs() < 1e-6);
        assert!(
            !all_at_minus_one,
            "elliptic zeros should not all be at z=-1"
        );
    }

    #[test]
    fn test_normalised_bp_order4() {
        let f = EllipticFilter::new(4, 1.0, 60.0)
            .bandpass(0.1, 0.4)
            .sample_rate(2.0);
        let sos = f.clone().sos().unwrap();
        let ba = f.clone().ba().unwrap();
        let zpk = f.zpk().unwrap();
        let expected_sos = [
            [
                1.55353890977011216e-02,
                2.65019241865987881e-02,
                1.55353890977011251e-02,
                1.00000000000000000e+00,
                -1.01185898156845511e+00,
                6.47165935354563904e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.99580104102376810e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -1.59659496996155026e+00,
                7.83953319520538261e-01,
            ],
            [
                1.00000000000000000e+00,
                8.52472860716617653e-01,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -5.83328423553997721e-01,
                8.54731655081845787e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.97880462190522444e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -1.85341401934793981e+00,
                9.49992632344384691e-01,
            ],
        ];
        let expected_b = vec![
            1.55353890977011216e-02,
            -2.20016237104663920e-02,
            -1.18847337351259859e-02,
            1.16676467497352720e-03,
            3.43830230901177175e-02,
            1.16676467497353067e-03,
            -1.18847337351259963e-02,
            -2.20016237104664059e-02,
            1.55353890977011355e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -5.04519639443194290e+00,
            1.22886520347314949e+01,
            -1.89164062199133589e+01,
            2.01400507241080255e+01,
            -1.51401312837241537e+01,
            7.84367027018817176e+00,
            -2.56798417291641590e+00,
            4.11960786240429810e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "ellip normalised bandpass n=4 sos");
        assert_vec_close(&ba.b, &expected_b, "ellip normalised bandpass n=4 b");
        assert_vec_close(&ba.a, &expected_a, "ellip normalised bandpass n=4 a");
        assert_eq!(zpk.z.len(), 8);
        assert_eq!(zpk.p.len(), 8);
        assert!(
            (zpk.k - 1.55353890977011216e-02_f64).abs() < 1e-7,
            "gain {}",
            zpk.k
        );
    }
}
