math_audio_wave/analytical/
mod.rs

1//! Analytical solutions for wave equation validation
2//!
3//! This module provides exact solutions to acoustic scattering problems
4//! used to validate numerical methods (BEM, FEM).
5//!
6//! ## Available Solutions
7//!
8//! - **1D**: Plane wave, standing wave, damped wave
9//! - **2D**: Cylinder scattering (Bessel/Hankel series)
10//! - **3D**: Sphere scattering (Mie theory)
11
12use num_complex::Complex64;
13use serde::{Deserialize, Serialize};
14
15pub mod solutions_1d;
16pub mod solutions_2d;
17pub mod solutions_3d;
18
19pub use solutions_1d::*;
20pub use solutions_2d::*;
21pub use solutions_3d::*;
22
23/// Point in space (1D, 2D, or 3D)
24#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq)]
25pub struct Point {
26    /// x-coordinate
27    pub x: f64,
28    /// y-coordinate (0 for 1D)
29    pub y: f64,
30    /// z-coordinate (0 for 1D/2D)
31    pub z: f64,
32}
33
34impl Point {
35    /// Create 1D point
36    pub fn new_1d(x: f64) -> Self {
37        Self { x, y: 0.0, z: 0.0 }
38    }
39
40    /// Create 2D point
41    pub fn new_2d(x: f64, y: f64) -> Self {
42        Self { x, y, z: 0.0 }
43    }
44
45    /// Create 3D point
46    pub fn new_3d(x: f64, y: f64, z: f64) -> Self {
47        Self { x, y, z }
48    }
49
50    /// Polar coordinates (r, θ) for 2D
51    pub fn from_polar(r: f64, theta: f64) -> Self {
52        Self::new_2d(r * theta.cos(), r * theta.sin())
53    }
54
55    /// Spherical coordinates (r, θ, φ) for 3D
56    /// θ = polar angle from z-axis
57    /// φ = azimuthal angle in xy-plane
58    pub fn from_spherical(r: f64, theta: f64, phi: f64) -> Self {
59        Self::new_3d(
60            r * theta.sin() * phi.cos(),
61            r * theta.sin() * phi.sin(),
62            r * theta.cos(),
63        )
64    }
65
66    /// Distance from origin
67    pub fn radius(&self) -> f64 {
68        (self.x * self.x + self.y * self.y + self.z * self.z).sqrt()
69    }
70
71    /// Polar angle (2D) - angle from positive x-axis
72    pub fn theta_2d(&self) -> f64 {
73        self.y.atan2(self.x)
74    }
75
76    /// Spherical polar angle (3D, from z-axis)
77    pub fn theta_3d(&self) -> f64 {
78        let r = self.radius();
79        if r < 1e-15 { 0.0 } else { (self.z / r).acos() }
80    }
81
82    /// Spherical azimuthal angle (3D)
83    pub fn phi_3d(&self) -> f64 {
84        self.y.atan2(self.x)
85    }
86
87    /// Distance to another point
88    pub fn distance_to(&self, other: &Point) -> f64 {
89        let dx = self.x - other.x;
90        let dy = self.y - other.y;
91        let dz = self.z - other.z;
92        (dx * dx + dy * dy + dz * dz).sqrt()
93    }
94}
95
96impl Default for Point {
97    fn default() -> Self {
98        Self::new_3d(0.0, 0.0, 0.0)
99    }
100}
101
102/// Analytical solution result
103#[derive(Debug, Clone, Serialize, Deserialize)]
104pub struct AnalyticalSolution {
105    /// Test name
106    pub name: String,
107    /// Dimensionality (1, 2, or 3)
108    pub dimensions: usize,
109    /// Evaluation points
110    pub positions: Vec<Point>,
111    /// Complex pressure values
112    pub pressure: Vec<Complex64>,
113    /// Wave number k = ω/c = 2πf/c
114    pub wave_number: f64,
115    /// Frequency (Hz)
116    pub frequency: f64,
117    /// Additional metadata
118    pub metadata: serde_json::Value,
119}
120
121impl AnalyticalSolution {
122    /// Create a new analytical solution
123    pub fn new(
124        name: impl Into<String>,
125        dimensions: usize,
126        positions: Vec<Point>,
127        pressure: Vec<Complex64>,
128        wave_number: f64,
129        frequency: f64,
130    ) -> Self {
131        Self {
132            name: name.into(),
133            dimensions,
134            positions,
135            pressure,
136            wave_number,
137            frequency,
138            metadata: serde_json::json!({}),
139        }
140    }
141
142    /// Compute pressure magnitude |p|
143    pub fn magnitude(&self) -> Vec<f64> {
144        self.pressure.iter().map(|p| p.norm()).collect()
145    }
146
147    /// Compute pressure phase (radians)
148    pub fn phase(&self) -> Vec<f64> {
149        self.pressure.iter().map(|p| p.arg()).collect()
150    }
151
152    /// Real part of pressure
153    pub fn real(&self) -> Vec<f64> {
154        self.pressure.iter().map(|p| p.re).collect()
155    }
156
157    /// Imaginary part of pressure
158    pub fn imag(&self) -> Vec<f64> {
159        self.pressure.iter().map(|p| p.im).collect()
160    }
161
162    /// Compute L2 error compared to another solution
163    pub fn l2_error(&self, other: &AnalyticalSolution) -> f64 {
164        assert_eq!(self.pressure.len(), other.pressure.len());
165
166        let sum_sq: f64 = self
167            .pressure
168            .iter()
169            .zip(other.pressure.iter())
170            .map(|(a, b)| (a - b).norm_sqr())
171            .sum();
172
173        sum_sq.sqrt()
174    }
175
176    /// Compute relative L2 error
177    pub fn relative_l2_error(&self, other: &AnalyticalSolution) -> f64 {
178        let l2_err = self.l2_error(other);
179        let norm: f64 = other
180            .pressure
181            .iter()
182            .map(|p| p.norm_sqr())
183            .sum::<f64>()
184            .sqrt();
185
186        if norm < 1e-15 { l2_err } else { l2_err / norm }
187    }
188
189    /// Compute max (L∞) error
190    pub fn linf_error(&self, other: &AnalyticalSolution) -> f64 {
191        assert_eq!(self.pressure.len(), other.pressure.len());
192
193        self.pressure
194            .iter()
195            .zip(other.pressure.iter())
196            .map(|(a, b)| (a - b).norm())
197            .fold(0.0, f64::max)
198    }
199}
200
201#[cfg(test)]
202mod tests {
203    use super::*;
204    use std::f64::consts::PI;
205
206    #[test]
207    fn test_point_creation() {
208        let p1d = Point::new_1d(1.0);
209        assert_eq!(p1d.x, 1.0);
210        assert_eq!(p1d.y, 0.0);
211        assert_eq!(p1d.z, 0.0);
212
213        let p2d = Point::from_polar(1.0, PI / 4.0);
214        assert!((p2d.radius() - 1.0).abs() < 1e-10);
215
216        let p3d = Point::from_spherical(1.0, PI / 2.0, 0.0);
217        assert!((p3d.radius() - 1.0).abs() < 1e-10);
218    }
219
220    #[test]
221    fn test_point_distance() {
222        let p1 = Point::new_3d(1.0, 0.0, 0.0);
223        let p2 = Point::new_3d(4.0, 4.0, 0.0);
224        assert!((p1.distance_to(&p2) - 5.0).abs() < 1e-10);
225    }
226
227    #[test]
228    fn test_solution_error() {
229        let sol1 = AnalyticalSolution::new(
230            "test1",
231            1,
232            vec![Point::new_1d(0.0)],
233            vec![Complex64::new(1.0, 0.0)],
234            1.0,
235            1.0,
236        );
237
238        let sol2 = AnalyticalSolution::new(
239            "test2",
240            1,
241            vec![Point::new_1d(0.0)],
242            vec![Complex64::new(1.1, 0.0)],
243            1.0,
244            1.0,
245        );
246
247        assert!((sol1.l2_error(&sol2) - 0.1).abs() < 1e-10);
248    }
249}
math_audio_wave/analytical/mod.rs

math_audio_wave/analytical/
mod.rs
