use std::{f64::consts::PI};
use num_complex::Complex;
use crate::filters::iir::lib::hyperbolic::{arcsinh, cosh, sinh};

pub fn cheb2ap(n: usize, rs: f64) -> Result<(Vec<Complex<f64>>, Vec<Complex<f64>>, f64), String> {
    if n == 0 {
        return Ok((vec![], vec![], 1.0));
    }

    let rs_linear = 10f64.powf(0.1 * rs);
    let de = 1.0 / (rs_linear - 1.0).sqrt();
    let mu = arcsinh(1.0 / de) / n as f64;

    // === Construct m ===
    let m: Vec<f64> = if n % 2 == 1 {
        let mut left: Vec<_> = (-((n as i32) - 1)..0).step_by(2).map(|x| x as f64).collect();
        let right: Vec<_> = (2..(n as i32)).step_by(2).map(|x| x as f64).collect();
        left.extend(right);
        left
    } else {
        (-((n as i32) - 1)..=(n as i32 - 1)).step_by(2).map(|x| x as f64).collect()
    };

    // === z = -conj(1j / sin(m*pi/(2N))) ===
    let z: Vec<Complex<f64>> = m.iter()
        .map(|&mi| {
            let theta = PI * mi / (2.0 * n as f64);
            let denom = theta.sin();
            let j_over_sin = Complex::new(0.0, 1.0) / denom;
            -j_over_sin.conj()
        })
        .collect();

    // === Poles like Butterworth ===
    let mut p: Vec<Complex<f64>> = (-((n as i32) - 1)..=(n as i32 - 1))
        .step_by(2)
        .map(|k| {
            let angle = PI * k as f64 / (2.0 * n as f64);
            let exp_jtheta = Complex::from_polar(1.0, angle);
            -exp_jtheta
        })
        .collect();

    // === Warp into Chebyshev II ===
    p = p.iter()
        .map(|&pole| {
            let warped = Complex::new(sinh(mu) * pole.re, cosh(mu) * pole.im);
            1.0 / warped
        })
        .collect();

    // === Gain k = prod(-p) / prod(-z) ===
    let p_prod: Complex<f64> = p.iter().map(|pi| -pi).product();
    let z_prod: Complex<f64> = z.iter().map(|zi| -zi).product();
    let k = (p_prod / z_prod).re;

    // Debug print (optional)
    // println!("z: {:?}", z);
    // println!("p: {:?}", p);
    // println!("k: {:?}", k);

    Ok((z, p, k))
}