rill_core_wdf/
elements.rs

1use crate::constants::{BOLTZMANN, ELECTRON_CHARGE, NEWTON_TOLERANCE};
2use crate::WdfElement;
3use rill_core::Transcendental;
4
5/// Resistor WDF element
6#[derive(Debug, Clone, Copy)]
7pub struct Resistor<T: Transcendental> {
8    resistance: T,
9    port_resistance: T,
10    voltage: T,
11    current: T,
12}
13
14impl<T: Transcendental> Resistor<T> {
15    /// Create a new resistor with given resistance in ohms
16    pub fn new(resistance: T) -> Self {
17        Self {
18            port_resistance: resistance,
19            resistance,
20            voltage: T::ZERO,
21            current: T::ZERO,
22        }
23    }
24
25    /// Get resistance value
26    pub fn resistance(&self) -> T {
27        self.resistance
28    }
29}
30
31impl<T: Transcendental> WdfElement<T> for Resistor<T> {
32    fn port_resistance(&self) -> T {
33        self.port_resistance
34    }
35
36    fn process_incident(&mut self, _a: T) -> T {
37        T::ZERO
38    }
39
40    fn update_state(&mut self) {
41        self.voltage = self.current * self.resistance;
42    }
43
44    fn voltage(&self) -> T {
45        self.voltage
46    }
47
48    fn current(&self) -> T {
49        self.current
50    }
51
52    fn reset(&mut self) {
53        self.voltage = T::ZERO;
54        self.current = T::ZERO;
55    }
56}
57
58/// Capacitor WDF element (trapezoidal integration)
59#[derive(Debug, Clone, Copy)]
60pub struct Capacitor<T: Transcendental> {
61    capacitance: T,
62    sample_rate: T,
63    port_resistance: T,
64    voltage: T,
65    current: T,
66    state: T,
67}
68
69impl<T: Transcendental> Capacitor<T> {
70    /// Create a new capacitor with given capacitance in farads and sample rate
71    pub fn new(capacitance: T, sample_rate: T) -> Self {
72        let two = T::from_f32(2.0);
73        let t = T::ONE / sample_rate;
74        let port_resistance = t / (two * capacitance);
75
76        Self {
77            capacitance,
78            sample_rate,
79            port_resistance,
80            voltage: T::ZERO,
81            current: T::ZERO,
82            state: T::ZERO,
83        }
84    }
85
86    /// Get capacitance value
87    pub fn capacitance(&self) -> T {
88        self.capacitance
89    }
90
91    /// Set capacitance and recompute port resistance
92    pub fn set_capacitance(&mut self, capacitance: T) {
93        self.capacitance = capacitance;
94        let two = T::from_f32(2.0);
95        let t = T::ONE / self.sample_rate;
96        self.port_resistance = t / (two * capacitance);
97    }
98
99    /// Set sample rate and recompute port resistance
100    pub fn set_sample_rate(&mut self, sample_rate: T) {
101        self.sample_rate = sample_rate;
102        let two = T::from_f32(2.0);
103        let t = T::ONE / sample_rate;
104        self.port_resistance = t / (two * self.capacitance);
105    }
106}
107
108impl<T: Transcendental> WdfElement<T> for Capacitor<T> {
109    fn port_resistance(&self) -> T {
110        self.port_resistance
111    }
112
113    fn process_incident(&mut self, a: T) -> T {
114        let two = T::from_f32(2.0);
115        let b = self.state - a;
116        self.voltage = (a + b) / two;
117        self.current = (a - b) / (two * self.port_resistance);
118        let next_state = self.voltage + self.port_resistance * self.current;
119        self.state = next_state;
120        b
121    }
122
123    fn update_state(&mut self) {
124        // state already updated in process_incident
125    }
126
127    fn voltage(&self) -> T {
128        self.voltage
129    }
130
131    fn current(&self) -> T {
132        self.current
133    }
134
135    fn reset(&mut self) {
136        self.voltage = T::ZERO;
137        self.current = T::ZERO;
138        self.state = T::ZERO;
139    }
140}
141
142/// Inductor WDF element (trapezoidal integration)
143#[derive(Debug, Clone, Copy)]
144pub struct Inductor<T: Transcendental> {
145    inductance: T,
146    sample_rate: T,
147    port_resistance: T,
148    voltage: T,
149    current: T,
150    state: T,
151}
152
153impl<T: Transcendental> Inductor<T> {
154    /// Create a new inductor with given inductance in henries and sample rate
155    pub fn new(inductance: T, sample_rate: T) -> Self {
156        let two = T::from_f32(2.0);
157        let t = T::ONE / sample_rate;
158        let port_resistance = two * inductance / t;
159
160        Self {
161            inductance,
162            sample_rate,
163            port_resistance,
164            voltage: T::ZERO,
165            current: T::ZERO,
166            state: T::ZERO,
167        }
168    }
169}
170
171impl<T: Transcendental> WdfElement<T> for Inductor<T> {
172    fn port_resistance(&self) -> T {
173        self.port_resistance
174    }
175
176    fn process_incident(&mut self, _a: T) -> T {
177        -self.state
178    }
179
180    fn update_state(&mut self) {
181        self.state = self.current * self.port_resistance;
182
183        let t = T::ONE / self.sample_rate;
184        self.current += self.voltage * t / self.inductance;
185    }
186
187    fn voltage(&self) -> T {
188        self.voltage
189    }
190
191    fn current(&self) -> T {
192        self.current
193    }
194
195    fn reset(&mut self) {
196        self.voltage = T::ZERO;
197        self.current = T::ZERO;
198        self.state = T::ZERO;
199    }
200}
201
202/// Diode WDF element (nonlinear, Newton-Raphson solution)
203#[derive(Debug, Clone, Copy)]
204pub struct Diode<T: Transcendental> {
205    pub(crate) saturation_current: T,
206    pub(crate) thermal_voltage: T,
207    pub(crate) ideality_factor: T,
208    pub(crate) port_resistance: T,
209    pub(crate) voltage: T,
210    pub(crate) current: T,
211    last_b: T,
212}
213
214impl<T: Transcendental> Diode<T> {
215    /// Create a new diode with Shockley parameters
216    ///
217    /// * `saturation_current` - Reverse saturation current Is (amperes)
218    /// * `ideality_factor` - Ideality factor n (1-2)
219    /// * `temperature_k` - Temperature in Kelvin
220    pub fn new(saturation_current: T, ideality_factor: T, temperature_k: T) -> Self {
221        let k = T::from_f64(BOLTZMANN);
222        let q = T::from_f64(ELECTRON_CHARGE);
223        let thermal_voltage = (k * temperature_k) / q;
224        let port_resistance = thermal_voltage / saturation_current;
225
226        Self {
227            saturation_current,
228            thermal_voltage,
229            ideality_factor,
230            port_resistance,
231            voltage: T::ZERO,
232            current: T::ZERO,
233            last_b: T::ZERO,
234        }
235    }
236
237    /// Get saturation current
238    pub fn saturation_current(&self) -> T {
239        self.saturation_current
240    }
241
242    /// Get thermal voltage
243    pub fn thermal_voltage(&self) -> T {
244        self.thermal_voltage
245    }
246
247    pub(crate) fn diode_equation(&self, v: T) -> T {
248        let vt = self.thermal_voltage * self.ideality_factor;
249        self.saturation_current * ((v / vt).exp() - T::ONE)
250    }
251
252    pub(crate) fn diode_derivative(&self, v: T) -> T {
253        let vt = self.thermal_voltage * self.ideality_factor;
254        self.saturation_current * (v / vt).exp() / vt
255    }
256
257    pub(crate) fn solve_newton(&self, a: T, r: T) -> T {
258        let vt = self.thermal_voltage * self.ideality_factor;
259        // Improved initial guess using simplified diode equation.
260        // For small a: v ≈ a / (1 + r*Is/vt)
261        // For large a: v ≈ vt * ln(a / (r*Is))
262        // Using a smoother approximation: v ≈ vt * ln(1 + a/(r*Is))
263        let guess = vt * (T::ONE + a / (r * self.saturation_current)).ln();
264        let mut v = guess.max(T::ZERO);
265        let tolerance = T::from_f64(NEWTON_TOLERANCE);
266
267        for _ in 0..10 {
268            let i = self.diode_equation(v);
269            let g = self.diode_derivative(v);
270
271            let f = v + r * i - a;
272
273            if f.abs() < tolerance {
274                break;
275            }
276
277            let df = T::ONE + r * g;
278            v -= f / df;
279        }
280
281        v
282    }
283}
284
285impl<T: Transcendental> WdfElement<T> for Diode<T> {
286    fn port_resistance(&self) -> T {
287        self.port_resistance
288    }
289
290    fn process_incident(&mut self, a: T) -> T {
291        let v = self.solve_newton(a, self.port_resistance);
292        let i = self.diode_equation(v);
293
294        self.voltage = v;
295        self.current = i;
296
297        T::from_f32(2.0) * v - a
298    }
299
300    fn update_state(&mut self) {
301        let g = self.diode_derivative(self.voltage);
302        if g > T::ZERO {
303            self.port_resistance = T::ONE / g;
304        }
305    }
306
307    fn voltage(&self) -> T {
308        self.voltage
309    }
310
311    fn current(&self) -> T {
312        self.current
313    }
314
315    fn reset(&mut self) {
316        self.voltage = T::ZERO;
317        self.current = T::ZERO;
318        self.last_b = T::ZERO;
319    }
320}
321
322#[cfg(test)]
323mod tests {
324    use super::*;
325
326    #[test]
327    fn test_resistor_wdf() {
328        let mut resistor: Resistor<f64> = Resistor::new(1000.0);
329        assert_eq!(resistor.port_resistance(), 1000.0);
330
331        let b = resistor.process_incident(1.0);
332        assert!((b - 0.0).abs() < 1e-10);
333    }
334
335    #[test]
336    fn test_capacitor_wdf() {
337        let sample_rate = 44100.0;
338        let capacitance = 1e-6;
339        let capacitor: Capacitor<f64> = Capacitor::new(capacitance, sample_rate);
340
341        let expected_r = 1.0 / (sample_rate * 2.0 * capacitance);
342        assert!((capacitor.port_resistance() - expected_r).abs() < 1e-10);
343    }
344
345    #[test]
346    fn test_inductor_wdf() {
347        let sample_rate = 44100.0;
348        let inductance = 100e-6;
349        let inductor: Inductor<f64> = Inductor::new(inductance, sample_rate);
350
351        let t = 1.0 / sample_rate;
352        let expected_r = 2.0 * inductance / t;
353        assert!((inductor.port_resistance() - expected_r).abs() < 1e-10);
354    }
355
356    #[test]
357    fn test_diode_thermal_voltage() {
358        let diode: Diode<f64> = Diode::new(1e-9, 1.0, 300.0);
359        let expected_vt = 1.380649e-23 * 300.0 / 1.60217662e-19;
360        assert!((diode.thermal_voltage() - expected_vt).abs() < 1e-15);
361    }
362}
rill_core_wdf/elements.rs

rill_core_wdf/
elements.rs
