1mod fft_backend;
4mod oscillator;
5mod response;
6
7use std::error::Error;
8use std::fmt;
9
10#[derive(Debug, Clone, Copy, PartialEq, Eq)]
11pub struct ShapeLayout {
12 dof_per_node: usize,
13 n_nodes: usize,
14}
15
16impl ShapeLayout {
17 pub fn new(dof_per_node: usize, n_nodes: usize) -> Result<Self, ModalSetError> {
18 if dof_per_node == 0 {
19 return Err(ModalSetError::InvalidDofPerNode { dof_per_node });
20 }
21 if n_nodes == 0 {
22 return Err(ModalSetError::InvalidNodeCount { n_nodes });
23 }
24 Ok(Self {
25 dof_per_node,
26 n_nodes,
27 })
28 }
29
30 pub fn dof_per_node(&self) -> usize {
31 self.dof_per_node
32 }
33
34 pub fn n_nodes(&self) -> usize {
35 self.n_nodes
36 }
37
38 pub fn shape_len(&self) -> usize {
39 self.dof_per_node * self.n_nodes
40 }
41}
42
43#[derive(Debug, Clone)]
44pub struct ModalRecord {
45 freq: f64,
46 weight: f64,
47 shape: Vec<f64>,
48 observed_damping: Option<f64>,
49 quality: Option<f64>,
50}
51
52impl ModalRecord {
53 pub fn new(
54 freq: f64,
55 weight: f64,
56 shape: Vec<f64>,
57 observed_damping: Option<f64>,
58 quality: Option<f64>,
59 ) -> Self {
60 Self {
61 freq,
62 weight,
63 shape,
64 observed_damping,
65 quality,
66 }
67 }
68
69 pub fn freq(&self) -> f64 {
70 self.freq
71 }
72
73 pub fn weight(&self) -> f64 {
74 self.weight
75 }
76
77 pub fn shape(&self) -> &[f64] {
78 &self.shape
79 }
80
81 pub fn observed_damping(&self) -> Option<f64> {
82 self.observed_damping
83 }
84
85 pub fn quality(&self) -> Option<f64> {
86 self.quality
87 }
88}
89
90#[derive(Debug, Clone, PartialEq)]
91pub enum ModalSetError {
92 InvalidDofPerNode {
93 dof_per_node: usize,
94 },
95 InvalidNodeCount {
96 n_nodes: usize,
97 },
98 InvalidShapeLength {
99 index: usize,
100 expected: usize,
101 actual: usize,
102 },
103 InvalidFrequency {
104 index: usize,
105 freq: f64,
106 },
107 InvalidWeight {
108 index: usize,
109 weight: f64,
110 },
111 InvalidObservedDamping {
112 index: usize,
113 value: f64,
114 },
115 InvalidQuality {
116 index: usize,
117 value: f64,
118 },
119 LayoutMismatch {
120 lhs: ShapeLayout,
121 rhs: ShapeLayout,
122 },
123 IndexOutOfBounds {
124 index: usize,
125 len: usize,
126 },
127 InvalidShapeNorm {
128 index: usize,
129 norm: f64,
130 },
131 InvalidComponentIndex {
132 index: usize,
133 dof_per_node: usize,
134 },
135 GammaLengthMismatch {
136 expected: usize,
137 actual: usize,
138 },
139 InvalidGamma {
140 index: usize,
141 value: f64,
142 },
143 InvalidNodeIndex {
144 index: usize,
145 n_nodes: usize,
146 },
147 InvalidDuration {
148 duration: f64,
149 },
150 InvalidSampleRate {
151 sample_rate: f64,
152 },
153 InvalidExcitationDelay {
154 delay: f64,
155 },
156 InvalidExcitationGain {
157 gain: f64,
158 },
159 InvalidExcitationPhase {
160 phase: f64,
161 },
162 SourceAmplitudeLengthMismatch {
163 expected: usize,
164 actual: usize,
165 },
166}
167
168impl fmt::Display for ModalSetError {
169 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
170 match self {
171 Self::InvalidDofPerNode { dof_per_node } => write!(f, "dof_per_node must be > 0, got {dof_per_node}"),
172 Self::InvalidNodeCount { n_nodes } => write!(f, "n_nodes must be > 0, got {n_nodes}"),
173 Self::InvalidShapeLength { index, expected, actual } => {
174 write!(f, "mode {index} shape length mismatch: expected {expected}, got {actual}")
175 }
176 Self::InvalidFrequency { index, freq } => write!(f, "mode {index} has invalid frequency {freq}"),
177 Self::InvalidWeight { index, weight } => write!(f, "mode {index} has invalid weight {weight}"),
178 Self::InvalidObservedDamping { index, value } => {
179 write!(f, "mode {index} has invalid observed damping {value}")
180 }
181 Self::InvalidQuality { index, value } => write!(f, "mode {index} has invalid quality {value}"),
182 Self::LayoutMismatch { lhs, rhs } => write!(
183 f,
184 "shape layout mismatch: lhs(dof_per_node={}, n_nodes={}), rhs(dof_per_node={}, n_nodes={})",
185 lhs.dof_per_node(),
186 lhs.n_nodes(),
187 rhs.dof_per_node(),
188 rhs.n_nodes()
189 ),
190 Self::IndexOutOfBounds { index, len } => {
191 write!(f, "index {index} out of bounds for modal set of length {len}")
192 }
193 Self::InvalidShapeNorm { index, norm } => {
194 write!(f, "mode {index} has invalid L2 shape norm {norm}")
195 }
196 Self::InvalidComponentIndex {
197 index,
198 dof_per_node,
199 } => write!(
200 f,
201 "component index {index} out of bounds for layout with dof_per_node={dof_per_node}"
202 ),
203 Self::GammaLengthMismatch { expected, actual } => {
204 write!(f, "gamma length mismatch: expected {expected}, got {actual}")
205 }
206 Self::InvalidGamma { index, value } => {
207 write!(f, "mode {index} has invalid gamma {value}")
208 }
209 Self::InvalidNodeIndex { index, n_nodes } => {
210 write!(f, "node index {index} out of bounds for modal layout with n_nodes={n_nodes}")
211 }
212 Self::InvalidDuration { duration } => {
213 write!(f, "duration must be finite and > 0, got {duration}")
214 }
215 Self::InvalidSampleRate { sample_rate } => {
216 write!(f, "sample_rate must be finite and > 0, got {sample_rate}")
217 }
218 Self::InvalidExcitationDelay { delay } => {
219 write!(f, "excitation delay must be finite and >= 0, got {delay}")
220 }
221 Self::InvalidExcitationGain { gain } => {
222 write!(f, "excitation gain must be finite, got {gain}")
223 }
224 Self::InvalidExcitationPhase { phase } => {
225 write!(f, "excitation phase must be finite, got {phase}")
226 }
227 Self::SourceAmplitudeLengthMismatch { expected, actual } => {
228 write!(f, "source amplitude length mismatch: expected {expected}, got {actual}")
229 }
230 }
231 }
232}
233
234impl Error for ModalSetError {}
235
236#[derive(Debug, Clone)]
237pub struct ModalSet {
238 modes: Vec<ModalRecord>,
239 layout: ShapeLayout,
240}
241
242impl ModalSet {
243 pub fn new(modes: Vec<ModalRecord>, layout: ShapeLayout) -> Result<Self, ModalSetError> {
244 for (index, mode) in modes.iter().enumerate() {
245 if mode.shape.len() != layout.shape_len() {
246 return Err(ModalSetError::InvalidShapeLength {
247 index,
248 expected: layout.shape_len(),
249 actual: mode.shape.len(),
250 });
251 }
252 if !mode.freq.is_finite() || mode.freq < 0.0 {
253 return Err(ModalSetError::InvalidFrequency {
254 index,
255 freq: mode.freq,
256 });
257 }
258 if !mode.weight.is_finite() {
259 return Err(ModalSetError::InvalidWeight {
260 index,
261 weight: mode.weight,
262 });
263 }
264 if let Some(value) = mode.observed_damping {
265 if !value.is_finite() {
266 return Err(ModalSetError::InvalidObservedDamping { index, value });
267 }
268 }
269 if let Some(value) = mode.quality {
270 if !value.is_finite() {
271 return Err(ModalSetError::InvalidQuality { index, value });
272 }
273 }
274 }
275 Ok(Self { modes, layout })
276 }
277
278 pub fn len(&self) -> usize {
279 self.modes.len()
280 }
281
282 pub fn is_empty(&self) -> bool {
283 self.modes.is_empty()
284 }
285
286 pub fn layout(&self) -> ShapeLayout {
287 self.layout
288 }
289
290 pub fn modes(&self) -> &[ModalRecord] {
291 &self.modes
292 }
293
294 pub fn iter(&self) -> impl Iterator<Item = &ModalRecord> {
295 self.modes.iter()
296 }
297
298 pub fn freqs(&self) -> Vec<f64> {
299 self.modes.iter().map(ModalRecord::freq).collect()
300 }
301
302 pub fn weights(&self) -> Vec<f64> {
303 self.modes.iter().map(ModalRecord::weight).collect()
304 }
305
306 pub fn shape(&self, index: usize) -> Option<&[f64]> {
307 self.modes.get(index).map(ModalRecord::shape)
308 }
309
310 pub fn shape_value(
311 &self,
312 mode_index: usize,
313 node_index: usize,
314 component_index: usize,
315 ) -> Result<f64, ModalSetError> {
316 if mode_index >= self.modes.len() {
317 return Err(ModalSetError::IndexOutOfBounds {
318 index: mode_index,
319 len: self.modes.len(),
320 });
321 }
322 if node_index >= self.layout.n_nodes() {
323 return Err(ModalSetError::InvalidNodeIndex {
324 index: node_index,
325 n_nodes: self.layout.n_nodes(),
326 });
327 }
328 if component_index >= self.layout.dof_per_node() {
329 return Err(ModalSetError::InvalidComponentIndex {
330 index: component_index,
331 dof_per_node: self.layout.dof_per_node(),
332 });
333 }
334 let offset = node_index * self.layout.dof_per_node() + component_index;
335 Ok(self.modes[mode_index].shape[offset])
336 }
337
338 pub fn filter_freq_min(&self, f_min: f64) -> Result<Self, ModalSetError> {
339 let modes = self
340 .modes
341 .iter()
342 .filter(|mode| mode.freq() >= f_min)
343 .cloned()
344 .collect();
345 Self::new(modes, self.layout)
346 }
347
348 pub fn subset(&self, indices: &[usize]) -> Result<Self, ModalSetError> {
349 let mut modes = Vec::with_capacity(indices.len());
350 for &index in indices {
351 let mode = self
352 .modes
353 .get(index)
354 .cloned()
355 .ok_or(ModalSetError::IndexOutOfBounds {
356 index,
357 len: self.modes.len(),
358 })?;
359 modes.push(mode);
360 }
361 Self::new(modes, self.layout)
362 }
363
364 pub fn sorted_by_freq(&self) -> Result<Self, ModalSetError> {
365 let mut modes = self.modes.clone();
366 modes.sort_by(|a, b| a.freq().total_cmp(&b.freq()));
367 Self::new(modes, self.layout)
368 }
369
370 pub fn merge(&self, other: &Self) -> Result<Self, ModalSetError> {
371 if self.layout != other.layout {
372 return Err(ModalSetError::LayoutMismatch {
373 lhs: self.layout,
374 rhs: other.layout,
375 });
376 }
377 let mut modes = self.modes.clone();
378 modes.extend(other.modes.iter().cloned());
379 modes.sort_by(|a, b| a.freq().total_cmp(&b.freq()));
380 Self::new(modes, self.layout)
381 }
382
383 pub fn normalized_shapes_l2(&self) -> Result<Self, ModalSetError> {
384 let mut modes = Vec::with_capacity(self.modes.len());
385 for (index, mode) in self.modes.iter().enumerate() {
386 let norm = mode
387 .shape
388 .iter()
389 .map(|value| value * value)
390 .sum::<f64>()
391 .sqrt();
392 if !norm.is_finite() || norm <= 0.0 {
393 return Err(ModalSetError::InvalidShapeNorm { index, norm });
394 }
395 let shape = mode.shape.iter().map(|value| value / norm).collect();
396 modes.push(ModalRecord::new(
397 mode.freq,
398 mode.weight,
399 shape,
400 mode.observed_damping,
401 mode.quality,
402 ));
403 }
404 Self::new(modes, self.layout)
405 }
406
407 pub fn component(&self, index: usize) -> Result<Self, ModalSetError> {
408 self.components(&[index])
409 }
410
411 pub fn components(&self, indices: &[usize]) -> Result<Self, ModalSetError> {
412 for &index in indices {
413 if index >= self.layout.dof_per_node() {
414 return Err(ModalSetError::InvalidComponentIndex {
415 index,
416 dof_per_node: self.layout.dof_per_node(),
417 });
418 }
419 }
420
421 let layout = ShapeLayout::new(indices.len(), self.layout.n_nodes())?;
422 let stride = self.layout.dof_per_node();
423 let mut modes = Vec::with_capacity(self.modes.len());
424 for mode in &self.modes {
425 let mut shape = Vec::with_capacity(layout.shape_len());
426 for node in 0..self.layout.n_nodes() {
427 let base = node * stride;
428 for &index in indices {
429 shape.push(mode.shape[base + index]);
430 }
431 }
432 modes.push(ModalRecord::new(
433 mode.freq,
434 mode.weight,
435 shape,
436 mode.observed_damping,
437 mode.quality,
438 ));
439 }
440 Self::new(modes, layout)
441 }
442
443 pub fn with_gammas(&self, gammas: &[f64]) -> Result<DampedModalSet, ModalSetError> {
444 if gammas.len() != self.modes.len() {
445 return Err(ModalSetError::GammaLengthMismatch {
446 expected: self.modes.len(),
447 actual: gammas.len(),
448 });
449 }
450 let modes = self
451 .modes
452 .iter()
453 .cloned()
454 .zip(gammas.iter().copied())
455 .map(|(base, gamma)| DampedModalRecord { base, gamma })
456 .collect();
457 DampedModalSet::new(modes, self.layout)
458 }
459}
460
461#[derive(Debug, Clone)]
462pub struct DampedModalRecord {
463 base: ModalRecord,
464 gamma: f64,
465}
466
467impl DampedModalRecord {
468 pub fn new(base: ModalRecord, gamma: f64) -> Self {
469 Self { base, gamma }
470 }
471
472 pub fn base(&self) -> &ModalRecord {
473 &self.base
474 }
475
476 pub fn gamma(&self) -> f64 {
477 self.gamma
478 }
479}
480
481#[derive(Debug, Clone)]
482pub struct DampedModalSet {
483 modes: Vec<DampedModalRecord>,
484 layout: ShapeLayout,
485}
486
487impl DampedModalSet {
488 pub fn new(modes: Vec<DampedModalRecord>, layout: ShapeLayout) -> Result<Self, ModalSetError> {
489 let base_modes: Vec<ModalRecord> = modes.iter().map(|mode| mode.base.clone()).collect();
490 let _ = ModalSet::new(base_modes, layout)?;
491 for (index, mode) in modes.iter().enumerate() {
492 if !mode.gamma.is_finite() {
493 return Err(ModalSetError::InvalidGamma {
494 index,
495 value: mode.gamma,
496 });
497 }
498 }
499 Ok(Self { modes, layout })
500 }
501
502 pub fn len(&self) -> usize {
503 self.modes.len()
504 }
505
506 pub fn is_empty(&self) -> bool {
507 self.modes.is_empty()
508 }
509
510 pub fn layout(&self) -> ShapeLayout {
511 self.layout
512 }
513
514 pub fn modes(&self) -> &[DampedModalRecord] {
515 &self.modes
516 }
517
518 pub fn iter(&self) -> impl Iterator<Item = &DampedModalRecord> {
519 self.modes.iter()
520 }
521
522 pub fn freqs(&self) -> Vec<f64> {
523 self.modes.iter().map(|mode| mode.base().freq()).collect()
524 }
525
526 pub fn weights(&self) -> Vec<f64> {
527 self.modes.iter().map(|mode| mode.base().weight()).collect()
528 }
529
530 pub fn gammas(&self) -> Vec<f64> {
531 self.modes.iter().map(DampedModalRecord::gamma).collect()
532 }
533
534 pub fn shape(&self, index: usize) -> Option<&[f64]> {
535 self.modes.get(index).map(|mode| mode.base().shape())
536 }
537
538 pub fn shape_value(
539 &self,
540 mode_index: usize,
541 node_index: usize,
542 component_index: usize,
543 ) -> Result<f64, ModalSetError> {
544 if mode_index >= self.modes.len() {
545 return Err(ModalSetError::IndexOutOfBounds {
546 index: mode_index,
547 len: self.modes.len(),
548 });
549 }
550 if node_index >= self.layout.n_nodes() {
551 return Err(ModalSetError::InvalidNodeIndex {
552 index: node_index,
553 n_nodes: self.layout.n_nodes(),
554 });
555 }
556 if component_index >= self.layout.dof_per_node() {
557 return Err(ModalSetError::InvalidComponentIndex {
558 index: component_index,
559 dof_per_node: self.layout.dof_per_node(),
560 });
561 }
562 let offset = node_index * self.layout.dof_per_node() + component_index;
563 Ok(self.modes[mode_index].base().shape()[offset])
564 }
565}
566
567pub use oscillator::OscillatorBank;
568pub use response::{
569 compute_source_amps, compute_source_amps_multi, generate_ir, generate_ir_multi,
570 SourceExcitation,
571};
572
573#[derive(Debug, Clone)]
574pub struct ModeInfo {
575 pub freq: f64,
576 pub damping_rate: f64,
577 pub amplitude: f64,
578 pub phase: f64,
579 pub quality: f64,
580}
581
582pub fn extract_modes(
583 readout: &[f64],
584 dt: f64,
585 n_modes_max: usize,
586 threshold_db: f64,
587) -> Vec<ModeInfo> {
588 let n = readout.len();
589 if n < 4 || !dt.is_finite() || dt <= 0.0 {
590 return Vec::new();
591 }
592
593 let spectrum = fft_backend::hann_window_spectrum(readout);
594 let power = spectrum.power();
595 let nfft = n.next_power_of_two();
596 let n_pos = nfft / 2;
597 let df = 1.0 / (nfft as f64 * dt);
598 let max_power = power.iter().cloned().fold(0.0_f64, f64::max);
599 if max_power < 1e-30 {
600 return Vec::new();
601 }
602
603 let threshold = max_power * 10.0_f64.powf(threshold_db / 10.0);
604 let min_bin = 3;
605 let mut peaks: Vec<(usize, f64)> = Vec::new();
606 for i in min_bin..n_pos - 1 {
607 if power[i] > power[i - 1] && power[i] > power[i + 1] && power[i] > threshold {
608 peaks.push((i, power[i]));
609 }
610 }
611 peaks.sort_by(|a, b| b.1.total_cmp(&a.1));
612 peaks.truncate(n_modes_max);
613
614 let mut modes = Vec::new();
615 for (bin, peak_power) in &peaks {
616 let bin = *bin;
617 let alpha_val = power[bin - 1];
618 let beta = power[bin];
619 let gamma = power[bin + 1];
620 let denom = alpha_val - 2.0 * beta + gamma;
621 let delta_bin = if denom.abs() > 1e-30 {
622 0.5 * (alpha_val - gamma) / denom
623 } else {
624 0.0
625 };
626 let freq_refined = (bin as f64 + delta_bin) * df;
627 let phase = spectrum.phase_at(bin);
628 let half_power = peak_power * 0.5;
629 let mut left_bin = bin;
630 while left_bin > min_bin && power[left_bin] > half_power {
631 left_bin -= 1;
632 }
633 let mut right_bin = bin;
634 while right_bin < n_pos - 2 && power[right_bin] > half_power {
635 right_bin += 1;
636 }
637 let bandwidth = (right_bin as f64 - left_bin as f64) * df;
638 let window_bw = 1.44 * df;
639 let signal_bw = (bandwidth - window_bw).max(df * 0.1);
640 let damping_rate = std::f64::consts::PI * signal_bw;
641 let amplitude = peak_power.sqrt();
642 let q = if damping_rate > 1e-10 {
643 std::f64::consts::PI * freq_refined / damping_rate
644 } else {
645 f64::INFINITY
646 };
647 modes.push(ModeInfo {
648 freq: freq_refined,
649 damping_rate,
650 amplitude,
651 phase,
652 quality: q,
653 });
654 }
655 modes.sort_by(|a, b| a.freq.total_cmp(&b.freq));
656 modes
657}
658
659pub fn extract_spatial_modes(
660 snapshots: &[&[f64]],
661 snapshot_times: &[f64],
662 mode_freqs: &[f64],
663 nx: usize,
664 ny: usize,
665) -> Vec<SpatialMode> {
666 let n_snaps = snapshots.len();
667 if n_snaps < 2 || snapshot_times.len() != n_snaps || mode_freqs.is_empty() || nx == 0 || ny == 0
668 {
669 return Vec::new();
670 }
671 let duration = snapshot_times[n_snaps - 1] - snapshot_times[0];
672 if !duration.is_finite() || duration <= 0.0 {
673 return Vec::new();
674 }
675 let n_cells = nx * ny;
676 let mut modes = Vec::new();
677 for &freq in mode_freqs {
678 if !freq.is_finite() {
679 continue;
680 }
681 let omega = 2.0 * std::f64::consts::PI * freq;
682 let mut shape_real = vec![0.0_f64; n_cells];
683 let mut shape_imag = vec![0.0_f64; n_cells];
684 for (k, snap) in snapshots.iter().enumerate() {
685 let t = snapshot_times[k] - snapshot_times[0];
686 let hann = 0.5 * (1.0 - (2.0 * std::f64::consts::PI * t / duration).cos());
687 let cos_val = (omega * t).cos() * hann;
688 let sin_val = (omega * t).sin() * hann;
689 for c in 0..n_cells.min(snap.len()) {
690 shape_real[c] += snap[c] * cos_val;
691 shape_imag[c] += snap[c] * sin_val;
692 }
693 }
694 let shape: Vec<f64> = shape_real
695 .iter()
696 .zip(shape_imag.iter())
697 .map(|(r, i)| (r * r + i * i).sqrt())
698 .collect();
699 let m_estimated = estimate_circumferential_order(&shape, nx, ny);
700 modes.push(SpatialMode {
701 freq,
702 shape,
703 m_estimated,
704 });
705 }
706 modes
707}
708
709#[derive(Debug, Clone)]
710pub struct SpatialMode {
711 pub freq: f64,
712 pub shape: Vec<f64>,
713 pub m_estimated: i32,
714}
715
716fn estimate_circumferential_order(shape: &[f64], nx: usize, ny: usize) -> i32 {
717 if nx == 0 || ny == 0 || shape.len() < nx.saturating_mul(ny) {
718 return 0;
719 }
720 let cx = nx as f64 / 2.0;
721 let cy = ny as f64 / 2.0;
722 let r_sample = (cx.min(cy)) * 0.6;
723 let n_samples = 64;
724 let mut values = Vec::with_capacity(n_samples);
725 for k in 0..n_samples {
726 let theta = 2.0 * std::f64::consts::PI * k as f64 / n_samples as f64;
727 let fi = cx + r_sample * theta.cos();
728 let fj = cy + r_sample * theta.sin();
729 let i = (fi.round() as usize).min(nx - 1);
730 let j = (fj.round() as usize).min(ny - 1);
731 values.push(shape[i * ny + j]);
732 }
733 let mean: f64 = values.iter().sum::<f64>() / values.len() as f64;
734 let mut crossings = 0;
735 for k in 0..values.len() {
736 let a = values[k] - mean;
737 let b = values[(k + 1) % values.len()] - mean;
738 if a * b < 0.0 {
739 crossings += 1;
740 }
741 }
742 crossings / 2
743}
744
745#[cfg(test)]
746mod tests {
747 use super::*;
748
749 fn sample_layout() -> ShapeLayout {
750 ShapeLayout::new(2, 3).unwrap()
751 }
752
753 fn sample_modes() -> Vec<ModalRecord> {
754 vec![
755 ModalRecord::new(
756 220.0,
757 0.8,
758 vec![1.0, 0.0, 0.5, 0.1, 0.2, 0.3],
759 Some(0.5),
760 Some(8.0),
761 ),
762 ModalRecord::new(110.0, 1.0, vec![0.2, 0.1, 0.3, 0.4, 0.0, 0.5], None, None),
763 ModalRecord::new(
764 330.0,
765 0.4,
766 vec![0.0, 0.3, 0.4, 0.2, 0.1, 0.6],
767 Some(0.2),
768 Some(12.0),
769 ),
770 ]
771 }
772
773 #[test]
774 fn modal_set_validates_shape_length() {
775 let layout = ShapeLayout::new(2, 2).unwrap();
776 let err = ModalSet::new(
777 vec![ModalRecord::new(
778 100.0,
779 1.0,
780 vec![1.0, 2.0, 3.0],
781 None,
782 None,
783 )],
784 layout,
785 )
786 .unwrap_err();
787 assert_eq!(
788 err,
789 ModalSetError::InvalidShapeLength {
790 index: 0,
791 expected: 4,
792 actual: 3
793 }
794 );
795 }
796
797 #[test]
798 fn shape_layout_rejects_zero_dimensions() {
799 assert_eq!(
800 ShapeLayout::new(0, 3).unwrap_err(),
801 ModalSetError::InvalidDofPerNode { dof_per_node: 0 }
802 );
803 assert_eq!(
804 ShapeLayout::new(2, 0).unwrap_err(),
805 ModalSetError::InvalidNodeCount { n_nodes: 0 }
806 );
807 }
808
809 #[test]
810 fn modal_set_rejects_invalid_frequency() {
811 let err = ModalSet::new(
812 vec![ModalRecord::new(-1.0, 1.0, vec![1.0; 6], None, None)],
813 sample_layout(),
814 )
815 .unwrap_err();
816 assert_eq!(
817 err,
818 ModalSetError::InvalidFrequency {
819 index: 0,
820 freq: -1.0
821 }
822 );
823 }
824
825 #[test]
826 fn modal_set_rejects_invalid_weight() {
827 let err = ModalSet::new(
828 vec![ModalRecord::new(10.0, f64::NAN, vec![1.0; 6], None, None)],
829 sample_layout(),
830 )
831 .unwrap_err();
832 assert!(matches!(
833 err,
834 ModalSetError::InvalidWeight {
835 index: 0,
836 weight
837 } if weight.is_nan()
838 ));
839 }
840
841 #[test]
842 fn modal_set_rejects_invalid_optional_quantities() {
843 let damping_err = ModalSet::new(
844 vec![ModalRecord::new(
845 10.0,
846 1.0,
847 vec![1.0; 6],
848 Some(f64::INFINITY),
849 None,
850 )],
851 sample_layout(),
852 )
853 .unwrap_err();
854 assert!(matches!(
855 damping_err,
856 ModalSetError::InvalidObservedDamping {
857 index: 0,
858 value
859 } if value.is_infinite()
860 ));
861
862 let quality_err = ModalSet::new(
863 vec![ModalRecord::new(
864 10.0,
865 1.0,
866 vec![1.0; 6],
867 None,
868 Some(f64::NAN),
869 )],
870 sample_layout(),
871 )
872 .unwrap_err();
873 assert!(matches!(
874 quality_err,
875 ModalSetError::InvalidQuality {
876 index: 0,
877 value
878 } if value.is_nan()
879 ));
880 }
881
882 #[test]
883 fn modal_set_sorts_by_frequency() {
884 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
885 let sorted = set.sorted_by_freq().unwrap();
886 assert_eq!(sorted.freqs(), vec![110.0, 220.0, 330.0]);
887 }
888
889 #[test]
890 fn modal_set_filters_and_subsets() {
891 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
892 let filtered = set.filter_freq_min(200.0).unwrap();
893 assert_eq!(filtered.freqs(), vec![220.0, 330.0]);
894 let subset = set.subset(&[2, 0]).unwrap();
895 assert_eq!(subset.freqs(), vec![330.0, 220.0]);
896 assert_eq!(subset.shape(1).unwrap(), &[1.0, 0.0, 0.5, 0.1, 0.2, 0.3]);
897 }
898
899 #[test]
900 fn modal_set_subset_rejects_out_of_bounds_index() {
901 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
902 let err = set.subset(&[0, 3]).unwrap_err();
903 assert_eq!(err, ModalSetError::IndexOutOfBounds { index: 3, len: 3 });
904 }
905
906 #[test]
907 fn modal_set_merge_requires_matching_layout() {
908 let lhs = ModalSet::new(sample_modes(), sample_layout()).unwrap();
909 let rhs = ModalSet::new(
910 vec![ModalRecord::new(
911 440.0,
912 0.2,
913 vec![1.0, 2.0, 3.0],
914 None,
915 None,
916 )],
917 ShapeLayout::new(1, 3).unwrap(),
918 )
919 .unwrap();
920 let err = lhs.merge(&rhs).unwrap_err();
921 assert_eq!(
922 err,
923 ModalSetError::LayoutMismatch {
924 lhs: sample_layout(),
925 rhs: ShapeLayout::new(1, 3).unwrap()
926 }
927 );
928 }
929
930 #[test]
931 fn modal_set_merge_keeps_frequency_order() {
932 let lhs = ModalSet::new(sample_modes()[..2].to_vec(), sample_layout()).unwrap();
933 let rhs = ModalSet::new(vec![sample_modes()[2].clone()], sample_layout()).unwrap();
934 let merged = lhs.merge(&rhs).unwrap();
935 assert_eq!(merged.freqs(), vec![110.0, 220.0, 330.0]);
936 assert_eq!(merged.weights(), vec![1.0, 0.8, 0.4]);
937 }
938
939 #[test]
940 fn modal_set_normalizes_shapes_without_changing_weight() {
941 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
942 let normalized = set.normalized_shapes_l2().unwrap();
943
944 assert_eq!(normalized.weights(), set.weights());
945 for mode in normalized.modes() {
946 let norm = mode
947 .shape()
948 .iter()
949 .map(|value| value * value)
950 .sum::<f64>()
951 .sqrt();
952 assert!((norm - 1.0).abs() < 1e-12, "norm={norm}");
953 }
954 }
955
956 #[test]
957 fn modal_set_normalization_rejects_zero_shape() {
958 let set = ModalSet::new(
959 vec![ModalRecord::new(10.0, 1.0, vec![0.0; 6], None, None)],
960 sample_layout(),
961 )
962 .unwrap();
963
964 assert_eq!(
965 set.normalized_shapes_l2().unwrap_err(),
966 ModalSetError::InvalidShapeNorm {
967 index: 0,
968 norm: 0.0,
969 }
970 );
971 }
972
973 #[test]
974 fn modal_set_component_extracts_one_component_per_node() {
975 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
976 let component = set.component(1).unwrap();
977
978 assert_eq!(component.layout(), ShapeLayout::new(1, 3).unwrap());
979 assert_eq!(component.shape(0).unwrap(), &[0.0, 0.1, 0.3]);
980 assert_eq!(component.shape(1).unwrap(), &[0.1, 0.4, 0.5]);
981 }
982
983 #[test]
984 fn modal_set_components_rebuilds_multi_component_shapes() {
985 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
986 let components = set.components(&[1, 0]).unwrap();
987
988 assert_eq!(components.layout(), sample_layout());
989 assert_eq!(
990 components.shape(0).unwrap(),
991 &[0.0, 1.0, 0.1, 0.5, 0.3, 0.2]
992 );
993 }
994
995 #[test]
996 fn modal_set_component_rejects_out_of_bounds_index() {
997 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
998 assert_eq!(
999 set.component(2).unwrap_err(),
1000 ModalSetError::InvalidComponentIndex {
1001 index: 2,
1002 dof_per_node: 2,
1003 }
1004 );
1005 }
1006
1007 #[test]
1008 fn modal_set_component_handles_scalar_layout() {
1009 let layout = ShapeLayout::new(1, 3).unwrap();
1010 let set = ModalSet::new(
1011 vec![ModalRecord::new(10.0, 1.0, vec![1.0, 2.0, 3.0], None, None)],
1012 layout,
1013 )
1014 .unwrap();
1015
1016 let component = set.component(0).unwrap();
1017 assert_eq!(component.layout(), layout);
1018 assert_eq!(component.shape(0).unwrap(), &[1.0, 2.0, 3.0]);
1019 }
1020
1021 #[test]
1022 fn modal_set_with_gammas_wraps_core_modes() {
1023 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
1024 let damped = set.with_gammas(&[0.1, 0.2, 0.3]).unwrap();
1025
1026 assert_eq!(damped.layout(), sample_layout());
1027 assert_eq!(damped.len(), 3);
1028 assert_eq!(damped.modes()[0].base().freq(), 220.0);
1029 assert_eq!(damped.modes()[2].gamma(), 0.3);
1030 }
1031
1032 #[test]
1033 fn modal_set_with_gammas_rejects_length_mismatch_and_invalid_values() {
1034 let set = ModalSet::new(sample_modes(), sample_layout()).unwrap();
1035 assert_eq!(
1036 set.with_gammas(&[0.1, 0.2]).unwrap_err(),
1037 ModalSetError::GammaLengthMismatch {
1038 expected: 3,
1039 actual: 2,
1040 }
1041 );
1042
1043 assert!(matches!(
1044 set.with_gammas(&[0.1, f64::NAN, 0.3]).unwrap_err(),
1045 ModalSetError::InvalidGamma {
1046 index: 1,
1047 value
1048 } if value.is_nan()
1049 ));
1050 }
1051
1052 #[test]
1053 fn extract_modes_detects_a_single_tone() {
1054 let dt = 1.0 / 8_000.0;
1055 let freq = 440.0;
1056 let samples = 4096;
1057 let readout: Vec<f64> = (0..samples)
1058 .map(|n| (2.0 * std::f64::consts::PI * freq * n as f64 * dt).sin())
1059 .collect();
1060 let modes = extract_modes(&readout, dt, 4, -40.0);
1061 assert!(!modes.is_empty());
1062 assert!((modes[0].freq - freq).abs() < 2.0, "freq={}", modes[0].freq);
1063 }
1064
1065 #[test]
1066 fn extract_modes_returns_empty_for_zero_signal() {
1067 let readout = vec![0.0; 1024];
1068 let modes = extract_modes(&readout, 1.0 / 44_100.0, 4, -40.0);
1069 assert!(modes.is_empty());
1070 }
1071
1072 #[test]
1073 fn extract_modes_returns_empty_for_short_signal_or_zero_mode_limit() {
1074 let short = vec![0.0, 1.0, 0.0];
1075 assert!(extract_modes(&short, 0.01, 4, -40.0).is_empty());
1076
1077 let dt = 1.0 / 8_000.0;
1078 let freq = 220.0;
1079 let readout: Vec<f64> = (0..1024)
1080 .map(|n| (2.0 * std::f64::consts::PI * freq * n as f64 * dt).sin())
1081 .collect();
1082 assert!(extract_modes(&readout, dt, 0, -40.0).is_empty());
1083 }
1084
1085 #[test]
1086 fn extract_spatial_modes_returns_one_shape_per_frequency() {
1087 let nx = 4;
1088 let ny = 3;
1089 let snapshots = [
1090 vec![0.0, 1.0, 0.5, 0.0, 0.2, 0.4, 0.0, 0.1, 0.3, 0.0, 0.2, 0.1],
1091 vec![0.2, 0.8, 0.7, 0.1, 0.3, 0.3, 0.1, 0.2, 0.2, 0.0, 0.1, 0.2],
1092 vec![0.4, 0.5, 0.8, 0.2, 0.4, 0.1, 0.2, 0.3, 0.1, 0.1, 0.0, 0.3],
1093 ];
1094 let refs: Vec<&[f64]> = snapshots.iter().map(Vec::as_slice).collect();
1095 let times = [0.0, 0.01, 0.02];
1096 let modes = extract_spatial_modes(&refs, ×, &[120.0, 240.0], nx, ny);
1097 assert_eq!(modes.len(), 2);
1098 assert!(modes.iter().all(|mode| mode.shape.len() == nx * ny));
1099 }
1100
1101 #[test]
1102 fn extract_spatial_modes_returns_empty_for_too_few_snapshots_or_no_frequencies() {
1103 let single = [vec![1.0, 0.0, 0.0, 1.0]];
1104 let single_refs: Vec<&[f64]> = single.iter().map(Vec::as_slice).collect();
1105 assert!(extract_spatial_modes(&single_refs, &[0.0], &[100.0], 2, 2).is_empty());
1106
1107 let snapshots = [vec![1.0, 0.0, 0.0, 1.0], vec![0.5, 0.0, 0.0, 0.5]];
1108 let refs: Vec<&[f64]> = snapshots.iter().map(Vec::as_slice).collect();
1109 assert!(extract_spatial_modes(&refs, &[0.0, 0.1], &[], 2, 2).is_empty());
1110 }
1111}