mkaudiolibrary 2.1.0

//! Nonlinear WDF ports for elements at the tree root.
//!
//! In a Wave Digital Filter tree, nonlinear elements (diodes, transistors,
//! tube stages) are placed at the root. The tree presents a single port
//! (with known port resistance R) to the nonlinear element, which must
//! solve for the reflected wave given the incident wave.
//!
//! # Wave variable convention
//!
//! This module uses the same convention as the rest of the WDF framework:
//! ```text
//! V = a + b
//! I = (a - b) / (2R)
//! ```
//!
//! From these: `b = V - a` and `I = (2a - V) / (2R)`.
//!
//! Given the element's I-V characteristic `I = f(V)`:
//! ```text
//! (2a - V) / (2R) = f(V)
//! V + 2R * f(V) = 2a
//!
//! Solve: g(V) = V + 2R * f(V) - 2a = 0
//!        g'(V) = 1 + 2R * f'(V)
//!
//! Once V is found: b = V - a.
//! ```

use crate::sim::core::solver::NewtonSolver;

/// Trait for nonlinear one-port WDF elements at the tree root.
///
/// Given the incident wave `a` from the WDF tree (with port resistance `R`),
/// compute and return the reflected wave `b`.
pub trait NonlinearPort {
    /// Solve for the reflected wave given the incident wave from the tree.
    ///
    /// `a` is the incident wave, `port_resistance` is the port resistance R
    /// of the tree looking into the root.
    fn solve(&mut self, a: f32, port_resistance: f32) -> f32;
}

// ---------------------------------------------------------------------------
// Diode
// ---------------------------------------------------------------------------

/// WDF diode port using the Shockley diode equation.
///
/// `I = Is * (exp(V / (n * Vt)) - 1)`
///
/// where:
/// - `Is`: saturation current (typically 1e-12 to 1e-6 A)
/// - `n`: ideality factor (1.0 for ideal, ~1.8 for silicon)
/// - `Vt`: thermal voltage (kT/q ~ 25.85 mV at room temperature)
///
/// The port solves the implicit equation `V + R * Is * (exp(V/(n*Vt)) - 1) = a`
/// using Newton-Raphson iteration.
pub struct DiodePort {
    /// Saturation current (A).
    pub is: f32,
    /// Ideality factor (dimensionless).
    pub n: f32,
    /// Thermal voltage (V). Default: 25.85 mV (room temperature).
    pub vt: f32,
    /// Newton solver.
    solver: NewtonSolver,
    /// Previous voltage — warm-start for Newton iteration.
    prev_v: f32,
}

impl DiodePort {
    /// Create a diode port with given Shockley parameters.
    pub fn new(is: f32, n: f32) -> Self {
        Self {
            is,
            n,
            vt: 0.02585, // 25.85 mV at room temperature
            solver: NewtonSolver {
                max_iterations: 32,
                tolerance: 1e-6,
                voltage_clamp: 20.0,
            },
            prev_v: 0.0,
        }
    }

    /// Create a silicon diode with typical parameters.
    pub fn silicon() -> Self {
        Self::new(1e-12, 1.8)
    }

    /// Create a germanium diode with typical parameters.
    pub fn germanium() -> Self {
        Self::new(1e-6, 1.3)
    }

    /// Create an LED with typical parameters.
    pub fn led() -> Self {
        Self::new(1e-18, 2.0)
    }

    /// Compute diode current for a given voltage.
    #[inline]
    fn diode_current(&self, v: f32) -> f32 {
        let exponent = v / (self.n * self.vt);
        // Clamp exponent to avoid overflow.
        let exponent = exponent.min(80.0);
        self.is * (exponent.exp() - 1.0)
    }

    /// Compute derivative of diode current w.r.t. voltage.
    #[inline]
    fn diode_current_derivative(&self, v: f32) -> f32 {
        let nvt = self.n * self.vt;
        let exponent = (v / nvt).min(80.0);
        self.is * exponent.exp() / nvt
    }

    /// Reset the warm-start voltage.
    pub fn reset(&mut self) {
        self.prev_v = 0.0;
    }
}

impl NonlinearPort for DiodePort {
    #[inline]
    fn solve(&mut self, a: f32, port_resistance: f32) -> f32 {
        // Using V = (a + b) convention:
        //   V + 2R * f(V) = 2a
        //
        // Solve: g(V) = V + 2R * I_d(V) - 2a = 0
        //        g'(V) = 1 + 2R * dI_d/dV
        let r2 = 2.0 * port_resistance;
        let nvt = self.n * self.vt;

        // For the diode, the forward voltage is typically 0.2-0.8V.
        // In reverse bias, V ≈ 2*a (all voltage drops across R).
        // Use a smart initial guess to help Newton converge.
        let init = if self.prev_v.abs() < 1e-6 {
            // First call: estimate based on whether forward or reverse biased
            if a > 0.0 {
                // Forward: V somewhere around n*Vt * ln(a/(R*Is))
                let rough = nvt * (a / (port_resistance * self.is).max(1.0)).ln();
                rough.clamp(0.0, 1.5)
            } else {
                // Reverse: diode passes ~zero current, V ≈ 2*a (all voltage on R)
                2.0 * a
            }
        } else {
            self.prev_v
        };

        // Newton iteration with adaptive step limiting for the exponential diode
        let mut v = init;
        for _ in 0..self.solver.max_iterations {
            let id = self.diode_current(v);
            let did_dv = self.diode_current_derivative(v);

            let g = v + r2 * id - 2.0 * a;
            let g_prime = (1.0 + r2 * did_dv).max(1e-6);

            if g.abs() < self.solver.tolerance {
                break;
            }

            let dv = g / g_prime;
            // Limit step to a few n*Vt to prevent exponential explosion
            let max_step = 10.0 * nvt;
            let dv = dv.clamp(-max_step, max_step);
            v -= dv;
            v = v.clamp(-20.0, 2.0); // forward diode < 2V, reverse can go further
        }

        self.prev_v = v;

        // Reflected wave: b = V - a
        v - a
    }
}

// ---------------------------------------------------------------------------
// Diode Pair (anti-parallel)
// ---------------------------------------------------------------------------

/// WDF anti-parallel diode pair (clipper).
///
/// Two diodes in anti-parallel: `I = Is * (exp(V/(n*Vt)) - exp(-V/(n*Vt)))`
/// which simplifies to `I = 2 * Is * sinh(V / (n * Vt))`.
///
/// This is the classic symmetrical soft-clipper used in overdrive pedals.
pub struct DiodePairPort {
    /// Saturation current (A).
    pub is: f32,
    /// Ideality factor (dimensionless).
    pub n: f32,
    /// Thermal voltage (V).
    pub vt: f32,
    /// Newton solver.
    solver: NewtonSolver,
    /// Previous voltage — warm-start.
    prev_v: f32,
}

impl DiodePairPort {
    /// Create an anti-parallel diode pair with given parameters.
    pub fn new(is: f32, n: f32) -> Self {
        Self {
            is,
            n,
            vt: 0.02585,
            solver: NewtonSolver {
                max_iterations: 8,
                tolerance: 1e-6,
                voltage_clamp: 100.0,
            },
            prev_v: 0.0,
        }
    }

    /// Create a silicon diode pair (1N4148-like).
    pub fn silicon() -> Self {
        Self::new(2.52e-9, 1.752)
    }

    /// Compute anti-parallel diode pair current.
    #[inline]
    fn pair_current(&self, v: f32) -> f32 {
        let x = v / (self.n * self.vt);
        let x_clamped = x.clamp(-80.0, 80.0);
        2.0 * self.is * x_clamped.sinh()
    }

    /// Derivative of anti-parallel diode pair current.
    #[inline]
    fn pair_current_derivative(&self, v: f32) -> f32 {
        let nvt = self.n * self.vt;
        let x = (v / nvt).clamp(-80.0, 80.0);
        2.0 * self.is * x.cosh() / nvt
    }

    /// Reset the warm-start voltage.
    pub fn reset(&mut self) {
        self.prev_v = 0.0;
    }
}

impl NonlinearPort for DiodePairPort {
    #[inline]
    fn solve(&mut self, a: f32, port_resistance: f32) -> f32 {
        // Using V = a + b convention:
        //   g(V) = V + 2R * I(V) - 2a = 0
        //   g'(V) = 1 + 2R * dI/dV
        let r2 = 2.0 * port_resistance;

        let v = self.solver.solve(self.prev_v, |v| {
            let i = self.pair_current(v);
            let di_dv = self.pair_current_derivative(v);

            let g = v + r2 * i - 2.0 * a;
            let g_prime = 1.0 + r2 * di_dv;

            (g, g_prime)
        });

        self.prev_v = v;
        // b = V - a
        v - a
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::sim::wdf::components::WdfComponent;

    #[test]
    fn test_diode_forward_bias() {
        let mut diode = DiodePort::silicon();

        // With a positive incident wave and moderate port resistance,
        // the diode should conduct and the reflected wave should be
        // less than the incident (energy absorbed).
        let a = 1.0;
        let r = 1000.0;
        let b = diode.solve(a, r);

        // V = a + b should be positive but small (diode forward voltage ~0.6V)
        let v = a + b;
        assert!(v > 0.0, "Diode voltage should be positive: {v}");
        assert!(v < 2.0, "Diode voltage should be bounded: {v}");

        // Reflected wave should have smaller magnitude than incident (diode absorbs power)
        assert!(b.abs() < a.abs(), "Diode should absorb power: a={a}, b={b}");
    }

    #[test]
    fn test_diode_reverse_bias() {
        let mut diode = DiodePort::silicon();

        // With a negative incident wave, the diode should block.
        let a = -1.0;
        let r = 1000.0;
        let b = diode.solve(a, r);

        // In reverse bias, very little current flows.
        // I = (a - b) / (2R) should be near zero.
        let i = (a - b) / (2.0 * r);
        assert!(i.abs() < 1e-6, "Reverse bias current should be tiny: {i}");
    }

    #[test]
    fn test_diode_pair_symmetry() {
        let mut pair = DiodePairPort::silicon();
        let r = 10000.0;

        // Anti-parallel pair should clip symmetrically.
        let b_pos = pair.solve(1.0, r);
        pair.reset();
        let b_neg = pair.solve(-1.0, r);

        let v_pos = 1.0 + b_pos; // V = a + b
        let v_neg = -1.0 + b_neg;

        // Should be symmetric: V(+a) = -V(-a)
        assert!(
            (v_pos + v_neg).abs() < 1e-4,
            "Diode pair not symmetric: v_pos={v_pos}, v_neg={v_neg}"
        );
    }

    #[test]
    fn test_diode_newton_convergence() {
        let mut diode = DiodePort::silicon();
        let r = 1000.0;

        // Run several samples to verify Newton converges consistently.
        for i in 0..100 {
            let a = 2.0 * (2.0 * std::f32::consts::PI * i as f32 / 100.0).sin();
            let b = diode.solve(a, r);

            // Basic sanity: b should be finite.
            assert!(
                b.is_finite(),
                "Diode solve returned non-finite: a={a}, b={b}"
            );

            // Verify the solution: V = a + b, I = (a - b)/(2R)
            // The diode equation: I_d(V) should equal I_port.
            let v = a + b;
            let i_diode = diode.diode_current(v);
            let i_port = (a - b) / (2.0 * r);
            let error = (i_diode - i_port).abs();
            assert!(
                error < 1e-3,
                "Diode equation not satisfied: I_diode={i_diode}, I_port={i_port}, error={error}"
            );
        }
    }

    #[test]
    fn test_diode_with_wdf_tree() {
        // Build a simple diode clipper circuit:
        // Vs --[R]--+--[diode]-- GND
        //           |
        //          Vout
        //
        // Using WDF: series(R, diode_port_as_leaf).
        // But since the diode is nonlinear, it goes at the root.
        // So the tree is just the resistor, and the diode is solved at the root.
        //
        // For this test, the "tree" is a single resistor.
        // The source drives through the resistor, and the diode clips.

        use crate::sim::wdf::components::WdfResistor;

        let r_val = 4700.0; // 4.7k — typical for a diode clipper
        let resistor = WdfResistor::new(r_val);
        // Use germanium diode (Is=1e-6) which has enough current to clip
        // Silicon (Is=1e-12) can't pass enough current to drop voltage across R
        let mut diode = DiodePort::germanium();

        let fs = 44100.0;
        let freq = 440.0;
        let amplitude = 5.0; // 5V peak — enough to forward-bias the diode

        let num_samples = (fs / freq * 2.0) as usize; // 2 cycles
        let mut max_output = 0.0f32;

        for n in 0..num_samples {
            let t = n as f32 / fs;
            let vs = amplitude * (2.0 * std::f32::consts::PI * freq * t).sin();

            let _b_tree = resistor.reflected();

            // The tree port resistance seen by the diode = R (the resistor).
            // The incident wave to the diode is simply Vs (Thevenin source).
            let a_to_diode = vs;

            // Diode solves and returns reflected wave.
            let b_from_diode = diode.solve(a_to_diode, r_val);

            // Output voltage (across diode) = a + b
            let v_out = a_to_diode + b_from_diode;

            // Only track forward-biased (positive) peaks for clipping check.
            // In reverse bias the diode blocks current, so V ≈ 2*a which is
            // correct physics but not relevant to the clipping assertion.
            if n > num_samples / 2 && v_out > 0.0 {
                max_output = max_output.max(v_out);
            }
        }

        // Forward-biased germanium diode should clip the positive peaks
        // well below the 5V source amplitude (~0.2-0.4V).
        assert!(
            max_output < 1.0,
            "Diode should clip forward output: max_output={max_output}"
        );
        assert!(
            max_output > 0.1,
            "Diode output should have some signal: max_output={max_output}"
        );
    }
}