dune_core/
srcfn.rs

1//! Time-dependent source functions (`dev/vsrc/vsrcload.c`,
2//! `vsrcacct.c`; the I-source mirrors them). Arguments are stored as
3//! given on the card; the defaults that depend on TSTEP and TSTOP
4//! are resolved at evaluation time, exactly like the C macros.
5
6use crate::matrix::AcceptCtx;
7
8/// A transient source function: `PULSE SIN EXP PWL SFFM`, and
9/// ngspice-47's transient-noise functions `TRNOISE TRRANDOM`, raw
10/// args. The noise functions are random: a source's realization is
11/// fixed by `.options seed` and the source name (see
12/// `docs/dune-core.md` §Ruído transiente), and is added to the
13/// source's `DC` value, as in ngspice.
14#[derive(Clone, Debug)]
15pub enum Wave {
16    Pulse(Vec<f64>),
17    Sin(Vec<f64>),
18    Exp(Vec<f64>),
19    Pwl(Vec<f64>),
20    Sffm(Vec<f64>),
21    /// `TRNOISE(NA NT <NALPHA NAMP <RTSAM RTSCAPT RTSEMT>>)`: white
22    /// gaussian noise of rms `NA` sampled every `NT` seconds (linear
23    /// between samples), plus `1/f^NALPHA` noise driven by white noise
24    /// of rms `NAMP`, plus random telegraph noise of amplitude
25    /// `RTSAM` with mean capture/emission times `RTSCAPT`/`RTSEMT`.
26    TrNoise(Vec<f64>),
27    /// `TRRANDOM(TYPE TS <TD <PARAM1 <PARAM2>>>)`: a random value held
28    /// for `TS` seconds from `TD` on — `TYPE` 1 uniform in
29    /// `[-PARAM1, PARAM1) + PARAM2`, 2 gaussian (σ `PARAM1`, mean
30    /// `PARAM2`), 3 exponential (mean `PARAM1`, offset `PARAM2`), 4
31    /// poisson (λ `PARAM1`, offset `PARAM2`); `PARAM1` defaults to 1,
32    /// `PARAM2` to 0, and the value is `PARAM2` up to `TD`.
33    TrRandom(Vec<f64>),
34}
35
36fn arg(c: &[f64], i: usize) -> Option<f64> {
37    c.get(i).copied()
38}
39
40/// An argument that falls back to a default when missing **or zero**
41/// (most of the C's macros test `&& coeffs[i] != 0`).
42fn arg_nz(c: &[f64], i: usize, default: f64) -> f64 {
43    match c.get(i) {
44        Some(&v) if v != 0.0 => v,
45        _ => default,
46    }
47}
48
49impl Wave {
50    /// Source value at `time` (`VSRCload`'s transient switch). The
51    /// noise functions give their deterministic `t = 0` value here
52    /// (0, or `PARAM2` for `TRRANDOM`); a source evaluates their
53    /// realization through [`crate::trnoise::NoiseWave`].
54    pub(crate) fn value(&self, time: f64, step: f64, final_time: f64) -> f64 {
55        match self {
56            Wave::TrNoise(_) => 0.0,
57            Wave::TrRandom(c) => arg(c, 4).unwrap_or(0.0),
58            Wave::Pulse(c) => {
59                let v1 = arg(c, 0).unwrap_or(0.0);
60                let v2 = arg(c, 1).unwrap_or(0.0);
61                let td = arg(c, 2).unwrap_or(0.0);
62                let tr = arg_nz(c, 3, step);
63                let tf = arg_nz(c, 4, step);
64                let pw = arg_nz(c, 5, final_time);
65                let per = arg_nz(c, 6, final_time);
66                let mut time = time - td;
67                if time > per {
68                    // Repeating signal - figure out where we are in
69                    // the period.
70                    let basetime = per * (time / per).floor();
71                    time -= basetime;
72                }
73                if time <= 0.0 || time >= tr + pw + tf {
74                    v1
75                } else if time >= tr && time <= tr + pw {
76                    v2
77                } else if time > 0.0 && time < tr {
78                    v1 + (v2 - v1) * time / tr
79                } else {
80                    v2 + (v1 - v2) * (time - (tr + pw)) / tf
81                }
82            }
83            Wave::Sin(c) => {
84                let vo = arg(c, 0).unwrap_or(0.0);
85                let va = arg(c, 1).unwrap_or(0.0);
86                let freq = arg_nz(c, 2, 1.0 / final_time);
87                let td = arg(c, 3).unwrap_or(0.0);
88                let theta = arg(c, 4).unwrap_or(0.0);
89                let time = time - td;
90                if time <= 0.0 {
91                    vo
92                } else {
93                    vo + va
94                        * (freq * time * 2.0 * std::f64::consts::PI).sin()
95                        * (-(time * theta)).exp()
96                }
97            }
98            Wave::Exp(c) => {
99                let v1 = arg(c, 0).unwrap_or(0.0);
100                let v2 = arg(c, 1).unwrap_or(0.0);
101                let td1 = arg(c, 2).unwrap_or(step);
102                let tau1 = arg_nz(c, 3, step);
103                let td2 = arg_nz(c, 4, td1 + step);
104                let tau2 = arg_nz(c, 5, step);
105                if time <= td1 {
106                    v1
107                } else if time <= td2 {
108                    v1 + (v2 - v1) * (1.0 - (-(time - td1) / tau1).exp())
109                } else {
110                    v1 + (v2 - v1) * (1.0 - (-(time - td1) / tau1).exp())
111                        + (v1 - v2) * (1.0 - (-(time - td2) / tau2).exp())
112                }
113            }
114            Wave::Sffm(c) => {
115                let vo = arg(c, 0).unwrap_or(0.0);
116                let va = arg(c, 1).unwrap_or(0.0);
117                let fc = arg_nz(c, 2, 1.0 / final_time);
118                let mdi = arg(c, 3).unwrap_or(0.0);
119                let fs = arg_nz(c, 4, 1.0 / final_time);
120                // The C spells pi 3.141592654 here (and only here).
121                #[allow(clippy::approx_constant)]
122                const PI: f64 = 3.141_592_654;
123                vo + va * ((2.0 * PI * fc * time) + mdi * (2.0 * PI * fs * time).sin()).sin()
124            }
125            Wave::Pwl(c) => {
126                let pairs = c.len() / 2;
127                if pairs == 0 {
128                    return 0.0;
129                }
130                if time < c[0] {
131                    return c[1];
132                }
133                for i in 0..pairs.saturating_sub(1) {
134                    let (t0, v0) = (c[2 * i], c[2 * i + 1]);
135                    let (t1, v1) = (c[2 * i + 2], c[2 * i + 3]);
136                    if t0 == time {
137                        return v0;
138                    }
139                    if t0 < time && t1 > time {
140                        return v0 + (time - t0) / (t1 - t0) * (v1 - v0);
141                    }
142                }
143                c[2 * pairs - 1]
144            }
145        }
146    }
147
148    /// The waveform's corners inside `(t0, t1)`: every PULSE corner
149    /// (`TD + k·PER + {0, TR, TR+PW, TR+PW+TF}`, defaults as in
150    /// [`Wave::value`]) and every PWL vertex time — the breakpoints
151    /// [`Wave::accept`] would have planted by `t0` had the transient
152    /// started at 0. At most 4 corners per period are listed, and a
153    /// pulse with more than 10 000 periods in the window lists none
154    /// (the step control then meets its corners unaided).
155    pub(crate) fn corners(&self, t0: f64, t1: f64, step: f64, final_time: f64, out: &mut Vec<f64>) {
156        let inside = |t: f64| t > t0 && t < t1;
157        match self {
158            Wave::Pulse(c) => {
159                let td = arg(c, 2).unwrap_or(0.0);
160                let tr = arg_nz(c, 3, step);
161                let tf = arg_nz(c, 4, step);
162                let pw = arg_nz(c, 5, final_time);
163                let per = arg_nz(c, 6, final_time);
164                let offsets = [0.0, tr, tr + pw, tr + pw + tf];
165                if !(per > 0.0 && per.is_finite()) {
166                    return;
167                }
168                let k0 = ((t0 - td) / per).floor() - 1.0;
169                let k1 = ((t1 - td) / per).ceil() + 1.0;
170                if !(k0.is_finite() && k1.is_finite()) || k1 - k0 > 10_000.0 {
171                    return;
172                }
173                let mut k = k0.max(0.0);
174                while k <= k1 {
175                    for o in offsets {
176                        let t = td + k * per + o;
177                        if inside(t) {
178                            out.push(t);
179                        }
180                    }
181                    k += 1.0;
182                }
183            }
184            Wave::Pwl(c) => {
185                // Every even index is a PWL corner time.
186                for &t in c.iter().step_by(2) {
187                    if inside(t) {
188                        out.push(t);
189                    }
190                }
191            }
192            _ => {}
193        }
194    }
195
196    /// Breakpoint generation at an accepted timepoint
197    /// (`VSRCaccept`); only PULSE and PWL make breakpoints in the
198    /// 3f5 ("no breakpoints (yet)" for the rest).
199    pub(crate) fn accept(&self, ctx: &mut AcceptCtx) {
200        match self {
201            Wave::Pulse(c) => {
202                let td = arg(c, 2).unwrap_or(0.0);
203                let tr = arg_nz(c, 3, ctx.step);
204                let tf = arg_nz(c, 4, ctx.step);
205                let pw = arg_nz(c, 5, ctx.final_time);
206                let per = arg_nz(c, 6, ctx.final_time);
207                let timetol = 1e-7 * pw;
208                let same = |a: f64, b: f64| (a - b).abs() <= timetol;
209
210                let mut time = ctx.time - td;
211                let mut basetime = 0.0;
212                if time >= per {
213                    basetime = per * (time / per).floor();
214                    time -= basetime;
215                }
216                if !ctx.brk {
217                    return;
218                }
219                if time <= 0.0 || time >= tr + pw + tf {
220                    if same(time, 0.0) {
221                        ctx.set_break(basetime + tr + td);
222                    } else if same(tr + pw + tf, time) {
223                        ctx.set_break(basetime + per + td);
224                    } else if time == -td {
225                        ctx.set_break(basetime + td);
226                    } else if same(per, time) {
227                        ctx.set_break(basetime + td + tr + per);
228                    }
229                } else if time >= tr && time <= tr + pw {
230                    if same(time, tr) {
231                        ctx.set_break(basetime + td + tr + pw);
232                    } else if same(tr + pw, time) {
233                        ctx.set_break(basetime + td + tr + pw + tf);
234                    }
235                } else if time > 0.0 && time < tr {
236                    if same(time, 0.0) {
237                        ctx.set_break(basetime + td + tr);
238                    } else if same(time, tr) {
239                        ctx.set_break(basetime + td + tr + pw);
240                    }
241                } else {
242                    // tr + pw < time < tr + pw + tf
243                    if same(time, tr + pw) {
244                        ctx.set_break(basetime + td + tr + pw + tf);
245                    } else if same(time, tr + pw + tf) {
246                        ctx.set_break(basetime + td + per);
247                    }
248                }
249            }
250            Wave::Pwl(c) => {
251                let pairs = c.len() / 2;
252                if pairs == 0 || !ctx.brk {
253                    return;
254                }
255                if ctx.time < c[0] {
256                    ctx.set_break(c[0]);
257                    return;
258                }
259                for i in 0..pairs.saturating_sub(1) {
260                    if c[2 * i] == ctx.time {
261                        ctx.set_break(c[2 * i + 2]);
262                        return;
263                    }
264                }
265            }
266            _ => {}
267        }
268    }
269}
270
271#[cfg(test)]
272mod tests {
273    use super::*;
274
275    #[test]
276    fn pulse_shape_and_defaults() {
277        // pulse(0 1) with tstep 0.1, tstop 7 (rc.cir): TR=TF=0.1,
278        // PW=PER=7.
279        let w = Wave::Pulse(vec![0.0, 1.0]);
280        assert_eq!(w.value(0.0, 0.1, 7.0), 0.0);
281        assert!((w.value(0.05, 0.1, 7.0) - 0.5).abs() < 1e-12);
282        assert_eq!(w.value(0.1, 0.1, 7.0), 1.0);
283        assert_eq!(w.value(5.0, 0.1, 7.0), 1.0);
284    }
285
286    #[test]
287    fn pulse_full_period() {
288        let w = Wave::Pulse(vec![-1.0, 1.0, 1.0, 0.5, 0.5, 2.0, 5.0]);
289        assert_eq!(w.value(0.0, 0.1, 100.0), -1.0); // before TD
290        assert_eq!(w.value(1.25, 0.1, 100.0), 0.0); // mid rise
291        assert_eq!(w.value(2.0, 0.1, 100.0), 1.0); // on
292        assert_eq!(w.value(3.75, 0.1, 100.0), 0.0); // mid fall
293        assert_eq!(w.value(5.0, 0.1, 100.0), -1.0); // off
294        // Second period.
295        assert_eq!(w.value(7.0, 0.1, 100.0), 1.0);
296    }
297
298    #[test]
299    fn pwl_segments() {
300        let w = Wave::Pwl(vec![0.0, -7.0, 10e-9, -7.0, 11e-9, -3.0]);
301        assert_eq!(w.value(-1.0, 0.0, 1.0), -7.0);
302        assert_eq!(w.value(5e-9, 0.0, 1.0), -7.0);
303        assert!((w.value(10.5e-9, 0.0, 1.0) - (-5.0)).abs() < 1e-9);
304        assert_eq!(w.value(2e-8, 0.0, 1.0), -3.0); // past the end
305    }
306
307    #[test]
308    fn sin_and_exp() {
309        let w = Wave::Sin(vec![1.0, 2.0, 1.0]); // 1 Hz around 1 V
310        assert_eq!(w.value(0.0, 0.0, 1.0), 1.0);
311        assert!((w.value(0.25, 0.0, 1.0) - 3.0).abs() < 1e-9);
312        let w = Wave::Exp(vec![0.0, 1.0, 0.1, 0.2, 0.5, 0.3]);
313        assert_eq!(w.value(0.05, 0.0, 1.0), 0.0);
314        let v = w.value(0.3, 0.0, 1.0);
315        assert!((v - (1.0 - (-0.2f64 / 0.2).exp())).abs() < 1e-12);
316    }
317}
dune_core/srcfn.rs

dune_core/
srcfn.rs
