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)]
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)]
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        self.state - a
115    }
116
117    fn update_state(&mut self) {
118        self.state = -self.current * self.port_resistance;
119
120        let t = T::ONE / self.sample_rate;
121        self.voltage += self.current * t / self.capacitance;
122    }
123
124    fn voltage(&self) -> T {
125        self.voltage
126    }
127
128    fn current(&self) -> T {
129        self.current
130    }
131
132    fn reset(&mut self) {
133        self.voltage = T::ZERO;
134        self.current = T::ZERO;
135        self.state = T::ZERO;
136    }
137}
138
139/// Inductor WDF element (trapezoidal integration)
140#[derive(Debug, Clone)]
141pub struct Inductor<T: Transcendental> {
142    inductance: T,
143    sample_rate: T,
144    port_resistance: T,
145    voltage: T,
146    current: T,
147    state: T,
148}
149
150impl<T: Transcendental> Inductor<T> {
151    /// Create a new inductor with given inductance in henries and sample rate
152    pub fn new(inductance: T, sample_rate: T) -> Self {
153        let two = T::from_f32(2.0);
154        let t = T::ONE / sample_rate;
155        let port_resistance = two * inductance / t;
156
157        Self {
158            inductance,
159            sample_rate,
160            port_resistance,
161            voltage: T::ZERO,
162            current: T::ZERO,
163            state: T::ZERO,
164        }
165    }
166}
167
168impl<T: Transcendental> WdfElement<T> for Inductor<T> {
169    fn port_resistance(&self) -> T {
170        self.port_resistance
171    }
172
173    fn process_incident(&mut self, _a: T) -> T {
174        -self.state
175    }
176
177    fn update_state(&mut self) {
178        self.state = self.current * self.port_resistance;
179
180        let t = T::ONE / self.sample_rate;
181        self.current += self.voltage * t / self.inductance;
182    }
183
184    fn voltage(&self) -> T {
185        self.voltage
186    }
187
188    fn current(&self) -> T {
189        self.current
190    }
191
192    fn reset(&mut self) {
193        self.voltage = T::ZERO;
194        self.current = T::ZERO;
195        self.state = T::ZERO;
196    }
197}
198
199/// Diode WDF element (nonlinear, Newton-Raphson solution)
200#[derive(Debug, Clone)]
201pub struct Diode<T: Transcendental> {
202    saturation_current: T,
203    thermal_voltage: T,
204    ideality_factor: T,
205    port_resistance: T,
206    voltage: T,
207    current: T,
208    last_b: T,
209}
210
211impl<T: Transcendental> Diode<T> {
212    /// Create a new diode with Shockley parameters
213    ///
214    /// * `saturation_current` - Reverse saturation current Is (amperes)
215    /// * `ideality_factor` - Ideality factor n (1-2)
216    /// * `temperature_k` - Temperature in Kelvin
217    pub fn new(saturation_current: T, ideality_factor: T, temperature_k: T) -> Self {
218        let k = T::from_f64(BOLTZMANN);
219        let q = T::from_f64(ELECTRON_CHARGE);
220        let thermal_voltage = (k * temperature_k) / q;
221        let port_resistance = thermal_voltage / saturation_current;
222
223        Self {
224            saturation_current,
225            thermal_voltage,
226            ideality_factor,
227            port_resistance,
228            voltage: T::ZERO,
229            current: T::ZERO,
230            last_b: T::ZERO,
231        }
232    }
233
234    /// Get saturation current
235    pub fn saturation_current(&self) -> T {
236        self.saturation_current
237    }
238
239    /// Get thermal voltage
240    pub fn thermal_voltage(&self) -> T {
241        self.thermal_voltage
242    }
243
244    fn diode_equation(&self, v: T) -> T {
245        let vt = self.thermal_voltage * self.ideality_factor;
246        self.saturation_current * ((v / vt).exp() - T::ONE)
247    }
248
249    fn diode_derivative(&self, v: T) -> T {
250        let vt = self.thermal_voltage * self.ideality_factor;
251        self.saturation_current * (v / vt).exp() / vt
252    }
253
254    fn solve_newton(&self, a: T, r: T) -> T {
255        let mut v = T::ZERO;
256        let tolerance = T::from_f64(NEWTON_TOLERANCE);
257
258        for _ in 0..10 {
259            let i = self.diode_equation(v);
260            let g = self.diode_derivative(v);
261
262            let f = v + r * i - a;
263
264            if f.abs() < tolerance {
265                break;
266            }
267
268            let df = T::ONE + r * g;
269            v -= f / df;
270        }
271
272        v
273    }
274}
275
276impl<T: Transcendental> WdfElement<T> for Diode<T> {
277    fn port_resistance(&self) -> T {
278        self.port_resistance
279    }
280
281    fn process_incident(&mut self, a: T) -> T {
282        let v = self.solve_newton(a, self.port_resistance);
283        let i = self.diode_equation(v);
284
285        self.voltage = v;
286        self.current = i;
287
288        T::from_f32(2.0) * v - a
289    }
290
291    fn update_state(&mut self) {
292        let g = self.diode_derivative(self.voltage);
293        if g > T::ZERO {
294            self.port_resistance = T::ONE / g;
295        }
296    }
297
298    fn voltage(&self) -> T {
299        self.voltage
300    }
301
302    fn current(&self) -> T {
303        self.current
304    }
305
306    fn reset(&mut self) {
307        self.voltage = T::ZERO;
308        self.current = T::ZERO;
309        self.last_b = T::ZERO;
310    }
311}
312
313#[cfg(test)]
314mod tests {
315    use super::*;
316
317    #[test]
318    fn test_resistor_wdf() {
319        let mut resistor: Resistor<f64> = Resistor::new(1000.0);
320        assert_eq!(resistor.port_resistance(), 1000.0);
321
322        let b = resistor.process_incident(1.0);
323        assert!((b - 0.0).abs() < 1e-10);
324    }
325
326    #[test]
327    fn test_capacitor_wdf() {
328        let sample_rate = 44100.0;
329        let capacitance = 1e-6;
330        let capacitor: Capacitor<f64> = Capacitor::new(capacitance, sample_rate);
331
332        let expected_r = 1.0 / (sample_rate * 2.0 * capacitance);
333        assert!((capacitor.port_resistance() - expected_r).abs() < 1e-10);
334    }
335
336    #[test]
337    fn test_inductor_wdf() {
338        let sample_rate = 44100.0;
339        let inductance = 100e-6;
340        let inductor: Inductor<f64> = Inductor::new(inductance, sample_rate);
341
342        let t = 1.0 / sample_rate;
343        let expected_r = 2.0 * inductance / t;
344        assert!((inductor.port_resistance() - expected_r).abs() < 1e-10);
345    }
346
347    #[test]
348    fn test_diode_thermal_voltage() {
349        let diode: Diode<f64> = Diode::new(1e-9, 1.0, 300.0);
350        let expected_vt = 1.380649e-23 * 300.0 / 1.60217662e-19;
351        assert!((diode.thermal_voltage() - expected_vt).abs() < 1e-15);
352    }
353}
rill_core_wdf/elements.rs

rill_core_wdf/
elements.rs
