/*
 * // 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_exp10, f_sincospi, f_sinh};

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

/// Builder for a Chebyshev Type I IIR filter.
///
/// Chebyshev Type I filters have equiripple in the passband and a monotone
/// stopband.  The ripple `rp` controls the allowed passband variation in dB.
///
#[derive(Copy, Clone)]
pub struct Chebyshev1Filter {
    order: usize,
    /// Maximum passband ripple in dB.
    rp: f64,
    filter_type: FilterType,
    freqs: [f64; 2],
    /// Sample rate in Hz; defaults to 2.0 (normalised frequencies).
    fs: f64,
}

impl Chebyshev1Filter {
    /// Create a new builder.
    ///
    /// - `order` — filter order (≥ 1).
    /// - `rp`    — maximum passband ripple in dB (> 0).
    pub fn new(order: usize, rp: f64) -> Self {
        Self {
            order,
            rp,
            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 = filter_zpk(
            self.order,
            self.rp,
            self.filter_type,
            self.freqs,
            self.fs,
            cheb1ap,
        )?;
        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.rp,
            self.filter_type,
            self.freqs,
            self.fs,
            cheb1ap,
        )?;
        Ok(zpk2tf(&zpk))
    }

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

// ─── step 1: Chebyshev Type I analog prototype ───────────────────────────────

/// Analog Chebyshev Type I lowpass prototype (`scipy.signal.cheb1ap`).
///
/// Returns `(zeros=[], poles, gain)`.  Poles lie on an ellipse in the left
/// half-plane; the passband edge is at 1 rad/s with `rp` dB of ripple.
fn cheb1ap(order: usize, rp: f64) -> ZpkFilter {
    let n = order as f64;

    // Ripple factor ε and the hyperbolic parameter μ
    let eps = (f_exp10(0.1 * rp) - 1.0).sqrt();
    let mu = f_asinh(1.0 / eps) / n;

    // Poles: p = -sinh(mu + i * pi * m / (2N))  for m in -N+1, -N+3, ..., N-1
    let poles: Vec<Complex64> = (0..order)
        .map(|k| {
            let m = -(n 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_sincos = f_sincospi(theta);
            let re = -f_sinh(mu) * v_sincos.1;
            let im = -f_cosh(mu) * v_sincos.0;
            Complex64::new(re, im)
        })
        .collect();

    // Gain: real(prod(-p))
    // For even order, divide by sqrt(1 + eps^2) to normalise DC gain to -rp dB
    let prod_neg_p: Complex64 = poles
        .iter()
        .fold(Complex64::from(1.0), |acc, &pi| acc * (-pi));
    let mut k = prod_neg_p.re;
    if order.is_multiple_of(2) {
        k /= (1.0 + eps * eps).sqrt();
    }

    ZpkFilter {
        z: vec![],
        p: poles,
        k,
    }
}

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

    const ATOL: f64 = 1e-8;

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

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

    // -----------------------------------------------------------------------
    // Reproduce in Python:
    //   from scipy.signal import cheby1
    //   sos = cheby1(N, rp, freq, btype='...', fs=fs, output='sos')
    //   for row in sos: print(list(row))
    // -----------------------------------------------------------------------

    /// cheby1(2, 1.0, 100.0, btype='lowpass', fs=1000.0)
    #[test]
    fn test_lowpass_order2_1db() {
        let sos = Chebyshev1Filter::new(2, 1.0)
            .lowpass(100.0)
            .sample_rate(1000.0)
            .sos()
            .unwrap();
        let expected = [[
            7.04224575725001384e-02,
            1.40844915145000277e-01,
            7.04224575725001384e-02,
            1.00000000000000000e+00,
            -1.19967756820218363e+00,
            5.15738756176752577e-01,
        ]];
        assert_sos_close(&sos, &expected, "cheby1 lowpass order=2 rp=1");
    }

    /// cheby1(3, 3.0, 200.0, btype='highpass', fs=1000.0)
    #[test]
    fn test_highpass_order3_3db() {
        let sos = Chebyshev1Filter::new(3, 3.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .sos()
            .unwrap();
        let expected = [
            [
                1.54320524554328958e-01,
                -3.08641049108657917e-01,
                1.54320524554328958e-01,
                1.00000000000000000e+00,
                4.17418917594625072e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                -1.00000000000000000e+00,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                -3.93068569598192863e-01,
                7.26059966703282877e-01,
            ],
        ];
        assert_sos_close(&sos, &expected, "cheby1 highpass order=3 rp=3");
    }

    /// cheby1(4, 0.5, [100.0, 300.0], btype='bandpass', fs=1000.0)
    #[test]
    fn test_bandpass_order4_05db() {
        let sos = Chebyshev1Filter::new(4, 0.5)
            .bandpass(100.0, 300.0)
            .sample_rate(1000.0)
            .sos()
            .unwrap();
        let expected = [
            [
                3.04448700370295554e-02,
                6.08897400740591108e-02,
                3.04448700370295554e-02,
                1.00000000000000000e+00,
                -9.19508632663331715e-02,
                5.22430246448362556e-01,
            ],
            [
                1.00000000000000000e+00,
                -2.00000000000000000e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.01591342954771990e+00,
                6.08209120698525174e-01,
            ],
            [
                1.00000000000000000e+00,
                2.00000000000000000e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                5.83260079113768426e-01,
                8.15277357162817129e-01,
            ],
            [
                1.00000000000000000e+00,
                -2.00000000000000000e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.53167973823163917e+00,
                8.82411787238062129e-01,
            ],
        ];
        assert_sos_close(&sos, &expected, "cheby1 bandpass order=4 rp=0.5");
    }

    /// cheby1(2, 1.0, 0.1, btype='lowpass', fs=2.0)
    #[test]
    fn test_normalised_lp_order2() {
        let sos = Chebyshev1Filter::new(2, 1.0).lowpass(0.1).sos().unwrap();
        let expected = [[
            2.05152236317149468e-02,
            4.10304472634298936e-02,
            2.05152236317149468e-02,
            1.00000000000000000e+00,
            -1.61851963861553316e+00,
            7.10593476651196876e-01,
        ]];
        assert_sos_close(&sos, &expected, "cheby1 normalised lowpass order=2");
    }

    /// cheby1(5, 2.0, 0.2, btype='lowpass', fs=2.0)
    #[test]
    fn test_normalised_lp_order5() {
        let sos = Chebyshev1Filter::new(5, 2.0).lowpass(0.2).sos().unwrap();
        let expected = [
            [
                2.08849950408453047e-04,
                4.17699900816906095e-04,
                2.08849950408453047e-04,
                1.00000000000000000e+00,
                -8.67531031170745814e-01,
                0.00000000000000000e+00,
            ],
            [
                1.00000000000000000e+00,
                2.00000000000000000e+00,
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.65789632003404552e+00,
                8.01481074195483867e-01,
            ],
            [
                1.00000000000000000e+00,
                1.00000000000000000e+00,
                0.00000000000000000e+00,
                1.00000000000000000e+00,
                -1.57201523547386857e+00,
                9.23382963274086643e-01,
            ],
        ];
        assert_sos_close(&sos, &expected, "cheby1 normalised lowpass order=5");
    }

    /// cheby1(2, 1.0, 100.0, btype='lowpass', fs=1000.0, output='ba')
    #[test]
    fn test_ba_lowpass_order2_1db() {
        let ba = Chebyshev1Filter::new(2, 1.0)
            .lowpass(100.0)
            .sample_rate(1000.0)
            .ba()
            .unwrap();
        let expected_b = vec![
            7.04224575725001384e-02,
            1.40844915145000277e-01,
            7.04224575725001384e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.19967756820218363e+00,
            5.15738756176752577e-01,
        ];
        assert_vec_close(&ba.b, &expected_b, "cheby1 ba lowpass order=2 rp=1 b");
        assert_vec_close(&ba.a, &expected_a, "cheby1 ba lowpass order=2 rp=1 a");
    }

    /// cheby1(3, 3.0, 200.0, btype='highpass', fs=1000.0, output='ba')
    #[test]
    fn test_ba_highpass_order3_3db() {
        let ba = Chebyshev1Filter::new(3, 3.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .ba()
            .unwrap();
        let expected_b = vec![
            1.54320524554328958e-01,
            -4.62961573662986847e-01,
            4.62961573662986847e-01,
            -1.54320524554328958e-01,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            2.43503479964322089e-02,
            5.61985709841137715e-01,
            3.03071165410073839e-01,
        ];
        assert_vec_close(&ba.b, &expected_b, "cheby1 ba highpass order=3 rp=3 b");
        assert_vec_close(&ba.a, &expected_a, "cheby1 ba highpass order=3 rp=3 a");
    }

    /// cheby1(2, 1.0, 0.1, output='ba') — normalised
    #[test]
    fn test_ba_normalised_lp_order2() {
        let ba = Chebyshev1Filter::new(2, 1.0).lowpass(0.1).ba().unwrap();
        let expected_b = vec![
            2.05152236317149468e-02,
            4.10304472634298936e-02,
            2.05152236317149468e-02,
        ];
        let expected_a = vec![
            1.00000000000000000e+00,
            -1.61851963861553316e+00,
            7.10593476651196876e-01,
        ];
        assert_vec_close(&ba.b, &expected_b, "cheby1 ba normalised lowpass order=2 b");
        assert_vec_close(&ba.a, &expected_a, "cheby1 ba normalised lowpass order=2 a");
    }

    // ── ZPK output tests ─────────────────────────────────────────────────────

    /// cheby1(2, 1.0, 100.0, output='zpk', fs=1000.0): zeros at -1 for lowpass
    #[test]
    fn test_zpk_lowpass_order2_1db() {
        let zpk = Chebyshev1Filter::new(2, 1.0)
            .lowpass(100.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        assert_eq!(zpk.z.len(), 2);
        assert_eq!(zpk.p.len(), 2);
        for z in &zpk.z {
            assert!(
                (z.re + 1.0).abs() < 1e-10,
                "zero re should be -1, got {}",
                z.re
            );
            assert!(z.im.abs() < 1e-10, "zero im should be 0, got {}", z.im);
        }
        assert!(
            (zpk.k - 7.04224575725001384e-02).abs() < 1e-8,
            "gain mismatch: {}",
            zpk.k
        );
    }

    /// ZPK → BA round-trip: zpk2tf(cheby1_zpk) == cheby1_ba
    #[test]
    fn test_zpk_ba_roundtrip() {
        let ba = Chebyshev1Filter::new(3, 3.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .ba()
            .unwrap();
        let zpk = Chebyshev1Filter::new(3, 3.0)
            .highpass(200.0)
            .sample_rate(1000.0)
            .zpk()
            .unwrap();
        let ba_from_zpk = zpk2tf(&zpk);
        assert_vec_close(&ba_from_zpk.b, &ba.b, "cheby1 roundtrip b");
        assert_vec_close(&ba_from_zpk.a, &ba.a, "cheby1 roundtrip a");
    }
}
