/*
 * // 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::filter_zpk::filter_zpk;
use crate::zpk::{zpk2sos, zpk2tf};
use crate::{BaFilter, FilterType, ImpliirError, SosFilter, ZpkFilter};
use num_complex::Complex64;
use pxfm::{f_asinh, f_cosh, f_csc, f_exp10, f_sincospi, f_sinh};
use std::f64::consts::PI;

/// Builder for a Chebyshev Type II IIR filter.
///
/// Chebyshev Type II filters have a monotone passband and equiripple in the
/// stopband.  The attenuation `rs` specifies the minimum stopband rejection
/// in dB.  Unlike Type I, the passband is maximally flat.
///
/// Three output formats are available:
/// - [`.sos()`][Chebyshev2Filter::sos] — second-order sections (recommended)
/// - [`.ba()`][Chebyshev2Filter::ba]   — transfer function `(b, a)` polynomial coefficients
/// - [`.zpk()`][Chebyshev2Filter::zpk] — zeros, poles, gain
///
#[derive(Copy, Clone)]
pub struct Chebyshev2Filter {
    order: usize,
    /// 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 Chebyshev2Filter {
    /// Create a new builder.
    ///
    /// - `order` — filter order (≥ 1).
    /// - `rs`    — minimum stopband attenuation in dB (> 0).
    pub fn new(order: usize, rs: f64) -> Self {
        Self {
            order,
            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.
    ///
    /// Preferred for high-order filters due to superior numerical stability.
    pub fn sos(self) -> Result<Vec<SosFilter>, ImpliirError> {
        let zpk = filter_zpk(
            self.order,
            self.rs,
            self.filter_type,
            self.freqs,
            self.fs,
            cheb2ap,
        )?;
        Ok(zpk2sos(zpk))
    }

    /// Design the filter and return transfer function coefficients `(b, a)`.
    pub fn ba(self) -> Result<BaFilter, ImpliirError> {
        let zpk = filter_zpk(
            self.order,
            self.rs,
            self.filter_type,
            self.freqs,
            self.fs,
            cheb2ap,
        )?;
        Ok(zpk2tf(&zpk))
    }

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

// ─── Chebyshev Type II analog prototype ──────────────────────────────────────

/// Analog Chebyshev Type II lowpass prototype (`scipy.signal.cheb2ap`).
///
/// Returns `(zeros, poles, gain)`.  Unlike Type I, this prototype has finite
/// zeros on the imaginary axis (at `1 / sin(m*pi / 2N)` for the m values
/// that exclude 0) giving the equiripple stopband.  Poles lie on a
/// reciprocal-ellipse in the left half-plane.
fn cheb2ap(order: usize, rs: f64) -> ZpkFilter {
    let n = order as f64;

    // Stopband ripple factor and hyperbolic parameter
    let de = 1.0 / (f_exp10(0.1 * rs) - 1.0).sqrt();
    let mu = f_asinh(1.0 / de) / n;

    // Zeros: purely imaginary, at 1/sin(m*pi/(2N))
    // For odd order, m skips 0 (which would give infinity); for even it doesn't arise.
    // scipy: odd  → m = arange(-N+1, 0, 2) ++ arange(2, N, 2)
    //        even → m = arange(-N+1, N, 2)
    let m_zeros: Vec<i32> = if order % 2 == 1 {
        // odd: skip m=0
        let neg: Vec<i32> = (0..)
            .map(|i| -(order as i32) + 1 + 2 * i)
            .take_while(|&v| v < 0)
            .collect();
        let pos: Vec<i32> = (0..)
            .map(|i| 2 + 2 * i)
            .take_while(|&v| v < order as i32)
            .collect();
        neg.into_iter().chain(pos).collect()
    } else {
        // even: m = -N+1, -N+3, ..., N-1  (same as Butterworth / cheby1)
        (0..order)
            .map(|k| -(order as i32) + 1 + 2 * k as i32)
            .collect()
    };

    let zeros: Vec<Complex64> = m_zeros
        .iter()
        .map(|&m| {
            let angle = m as f64 * PI / (2.0 * n);
            // z = i / sin(m*pi/(2N))
            Complex64::new(0.0, f_csc(angle))
        })
        .collect();

    // Poles: p = -1 / sinh(mu + i*theta)  for m in arange(-N+1, N, 2)
    let poles: Vec<Complex64> = (0..order)
        .map(|k| {
            let m = -(order as i32) + 1 + 2 * k as i32;
            let theta = m as f64 / (2.0 * n);
            // sinh(mu + i*theta) = sinh(mu)cos(theta) + i*cosh(mu)sin(theta)
            let v_sincospi = f_sincospi(theta);
            let s_re = f_sinh(mu) * v_sincospi.1;
            let s_im = f_cosh(mu) * v_sincospi.0;
            let s = Complex64::new(s_re, s_im);
            -Complex64::from(1.0) / s
        })
        .collect();

    // Gain: real(prod(-p) / prod(-z))
    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 k = (prod_neg_p / prod_neg_z).re;

    ZpkFilter {
        z: zeros,
        p: poles,
        k,
    }
}

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

#[cfg(test)]
mod cheby2_tests {
    use super::*;

    const ATOL: f64 = 1e-8;

    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}"
            );
        }
    }

    /// cheby2(2, 20.0, 100.0, btype='lowpass', fs=1000.0)
    #[test]
    fn test_lowpass_order2_20db() {
        let f = Chebyshev2Filter::new(2, 20.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 = [[
            9.95953541992127961e-02,
            -1.29738704280005807e-01,
            9.95953541992127822e-02,
            1.00000000000000000e+00,
            -1.60992106133294777e+00,
            6.79373065451367442e-01,
        ]];
        let expected_b = vec![
            9.95953541992127961e-02,
            -1.29738704280005807e-01,
            9.95953541992127822e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.60992106133294777e+00,
            6.79373065451367442e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "cheby2 lowpass n=2 rs=20 sos");
        assert_vec_close(&ba.b, &expected_b, "cheby2 lowpass n=2 rs=20 b");
        assert_vec_close(&ba.a, &expected_a, "cheby2 lowpass n=2 rs=20 a");
        assert_eq!(zpk.z.len(), 2);
        assert_eq!(zpk.p.len(), 2);
        assert!(
            (zpk.k - 9.95953541992127961e-02_f64).abs() < 1e-8,
            "gain {} != 9.959535e-02",
            zpk.k
        );
    }

    /// cheby2(3, 40.0, 200.0, btype='highpass', fs=1000.0)
    #[test]
    fn test_highpass_order3_40db() {
        let f = Chebyshev2Filter::new(3, 40.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 = [
            [
                5.91133177592007217e-02,
                -5.11648711101531650e-02,
                5.91133177592007356e-02,
                1.00000000000000000e+00,
                3.46893245770543535e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                -1.00000000000000000e+00,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                9.46394556427146871e-01,
                4.65119881288944992e-01,
            ],
        ];
        let expected_b = vec![
            5.91133177592007217e-02,
            -1.10278188869353880e-01,
            1.10278188869353894e-01,
            -5.91133177592007356e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            1.29328780219769035e+00,
            7.93417760747531786e-01,
            1.61346945292732025e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "cheby2 highpass n=3 rs=40 sos");
        assert_vec_close(&ba.b, &expected_b, "cheby2 highpass n=3 rs=40 b");
        assert_vec_close(&ba.a, &expected_a, "cheby2 highpass n=3 rs=40 a");
        assert_eq!(zpk.z.len(), 3);
        assert_eq!(zpk.p.len(), 3);
        assert!(
            (zpk.k - 5.91133177592007217e-02_f64).abs() < 1e-8,
            "gain {} != 5.911332e-02",
            zpk.k
        );
    }

    /// cheby2(4, 30.0, [100.0, 300.0], btype='bandpass', fs=1000.0)
    #[test]
    fn test_bandpass_order4_30db() {
        let f = Chebyshev2Filter::new(4, 30.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 = [
            [
                6.74936575437130315e-02,
                1.06327580506793917e-01,
                6.74936575437130315e-02,
                1.00000000000000000e+00,
                -9.92677493575896963e-02,
                3.25871957848596239e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.90719254926909088e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -8.90551832192124859e-01,
                4.29716134350616286e-01,
            ],
            [
                1.00000000000000000e+00,
                7.07251932077786583e-01,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                1.23490866673509170e-01,
                7.54464387373831036e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.65129209139081312e+00,
                1.00000000000000022e+00,
                1.00000000000000000e+00,
                -1.28045531712827332e+00,
                8.18782887525892678e-01,
            ],
        ];
        let expected_b = vec![
            6.74936575437130315e-02,
            -8.61125435036808468e-02,
            9.50555408456695179e-03,
            -4.07427789865790271e-02,
            1.20854775472336431e-01,
            -4.07427789865790965e-02,
            9.50555408456692230e-03,
            -8.61125435036808329e-02,
            6.74936575437130593e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -2.14678403200447843e+00,
            3.40429997453508282e+00,
            -3.57499252524238553e+00,
            3.19337673180381953e+00,
            -1.97451450931632588e+00,
            1.00744060283057024e+00,
            -3.26743981124902438e-01,
            8.65039924926256643e-02,
        ];
        assert_sos_close(&sos, &expected_sos, "cheby2 bandpass n=4 rs=30 sos");
        assert_vec_close(&ba.b, &expected_b, "cheby2 bandpass n=4 rs=30 b");
        assert_vec_close(&ba.a, &expected_a, "cheby2 bandpass n=4 rs=30 a");
        assert_eq!(zpk.z.len(), 8);
        assert_eq!(zpk.p.len(), 8);
        assert!(
            (zpk.k - 6.74936575437130315e-02_f64).abs() < 1e-8,
            "gain {} != 6.749366e-02",
            zpk.k
        );
    }

    /// cheby2(2, 20.0, 0.1, btype='lowpass', fs=2.0) — normalised
    #[test]
    fn test_normalised_lp_order2() {
        let f = Chebyshev2Filter::new(2, 20.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 = [[
            9.54659671203059584e-02,
            -1.72688631608675869e-01,
            9.54659671203059307e-02,
            1.00000000000000000e+00,
            -1.80898117793047186e+00,
            8.27224480562407827e-01,
        ]];
        let expected_b = vec![
            9.54659671203059584e-02,
            -1.72688631608675869e-01,
            9.54659671203059307e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.80898117793047186e+00,
            8.27224480562407827e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "cheby2 normalised lowpass n=2 sos");
        assert_vec_close(&ba.b, &expected_b, "cheby2 normalised lowpass n=2 b");
        assert_vec_close(&ba.a, &expected_a, "cheby2 normalised lowpass n=2 a");
        assert_eq!(zpk.z.len(), 2);
        assert_eq!(zpk.p.len(), 2);
        assert!(
            (zpk.k - 9.54659671203059584e-02_f64).abs() < 1e-8,
            "gain {} != 9.546597e-02",
            zpk.k
        );
    }

    /// cheby2(5, 60.0, 0.2, btype='lowpass', fs=2.0) — normalised, odd order
    #[test]
    fn test_normalised_lp_order5() {
        let f = Chebyshev2Filter::new(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.50921188606421451e-03,
                1.50921188606421451e-03,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                -7.40283301445591180e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                -1.06378929801971900e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.57203360371094969e+00,
                6.38684872471909193e-01,
            ],
            [
                1.00000000000000000e+00,
                -1.58192352099433586e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.79043019785819801e+00,
                8.58679977595713351e-01,
            ],
        ];
        let expected_b = vec![
            1.50921188606421451e-03,
            -2.48372934750425720e-03,
            1.56523457517359363e-03,
            1.56523457517359428e-03,
            -2.48372934750425633e-03,
            1.50921188606421451e-03,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -4.10274710301473888e+00,
            6.80115709021645642e+00,
            -5.68548220452255837e+00,
            2.39424419627079876e+00,
            -4.05990544722491076e-01,
        ];
        assert_sos_close(&sos, &expected_sos, "cheby2 normalised lowpass n=5 sos");
        assert_vec_close(&ba.b, &expected_b, "cheby2 normalised lowpass n=5 b");
        assert_vec_close(&ba.a, &expected_a, "cheby2 normalised lowpass n=5 a");
        assert_eq!(zpk.z.len(), 5);
        assert_eq!(zpk.p.len(), 5);
        assert!(
            (zpk.k - 1.50921188606421451e-03_f64).abs() < 1e-8,
            "gain {} != 1.509212e-03",
            zpk.k
        );
    }

    /// zpk2tf(cheby2_zpk) must equal cheby2_ba for any filter
    #[test]
    fn test_zpk_ba_roundtrip() {
        let ba = Chebyshev2Filter::new(3, 40.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .ba()
            .unwrap();
        let zpk = Chebyshev2Filter::new(3, 40.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        let ba_rt = zpk2tf(&zpk);
        assert_vec_close(&ba_rt.b, &ba.b, "cheby2 roundtrip b");
        assert_vec_close(&ba_rt.a, &ba.a, "cheby2 roundtrip a");
    }

    /// Type II has finite zeros (not all at ±1), verify they are non-trivial
    #[test]
    fn test_zeros_are_finite_and_imaginary() {
        // cheby2 lowpass zeros are purely imaginary: z = ±i/sin(m*pi/2N)
        let zpk = Chebyshev2Filter::new(4, 40.0)
            .lowpass(100.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        assert_eq!(zpk.z.len(), 4);
        for z in &zpk.z {
            // After bilinear transform, zeros are no longer purely imaginary,
            // but they must NOT all be at -1 (which would be the Butterworth case)
            assert!(
                (z.re + 1.0).abs() > 1e-3 || z.im.abs() > 1e-3,
                "zero {z:?} looks like a Butterworth zero at -1"
            );
        }
    }
}
