Skip to main content

oxiphysics_materials/acoustics/
types.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5#[allow(unused_imports)]
6use super::functions::*;
7#[allow(unused_imports)]
8use super::functions_2::*;
9use std::f64::consts::PI;
10
11/// Frequency weighting filters for noise measurement.
12#[derive(Debug, Clone, Copy, PartialEq, Eq)]
13#[allow(dead_code)]
14pub enum WeightingFilter {
15    /// A-weighting (most common, approximates human hearing sensitivity)
16    A,
17    /// B-weighting (intermediate, less used)
18    B,
19    /// C-weighting (flat response above ~30 Hz, used for peak measurements)
20    C,
21}
22/// Porous absorber model: rigid-frame porous material (Delany-Bazley model).
23///
24/// Empirical complex impedance and propagation constant for fibrous materials
25/// as a function of resistivity R_f (Pa·s/m²) and frequency f (Hz).
26///
27/// Reference: Delany & Bazley, Appl. Acoust. 3, 105 (1970).
28#[allow(dead_code)]
29#[derive(Debug, Clone)]
30pub struct PorousAbsorber {
31    /// Flow resistivity (Pa·s/m²), typically 5000–50000 for common absorbers.
32    pub flow_resistivity: f64,
33    /// Thickness of the absorber layer (m).
34    pub thickness: f64,
35}
36impl PorousAbsorber {
37    /// Create a new porous absorber.
38    pub fn new(flow_resistivity: f64, thickness: f64) -> Self {
39        Self {
40            flow_resistivity,
41            thickness,
42        }
43    }
44    /// Dimensionless frequency parameter X = ρ_air * f / R_f.
45    fn x_param(&self, freq: f64) -> f64 {
46        1.204 * freq / self.flow_resistivity
47    }
48    /// Real part of normalised characteristic impedance (Delany-Bazley).
49    pub fn impedance_real(&self, freq: f64) -> f64 {
50        let x = self.x_param(freq);
51        1.0 + 9.08 * x.powf(-0.75)
52    }
53    /// Imaginary part of normalised characteristic impedance (Delany-Bazley).
54    pub fn impedance_imag(&self, freq: f64) -> f64 {
55        let x = self.x_param(freq);
56        -11.9 * x.powf(-0.73)
57    }
58    /// Real part of propagation constant k_r = ω/c * α_r.
59    pub fn propagation_real(&self, freq: f64) -> f64 {
60        let x = self.x_param(freq);
61        let omega = 2.0 * PI * freq;
62        let c_air = 343.0;
63        (omega / c_air) * (1.0 + 10.8 * x.powf(-0.70))
64    }
65    /// Imaginary (attenuation) part of propagation constant.
66    pub fn propagation_imag(&self, freq: f64) -> f64 {
67        let x = self.x_param(freq);
68        let omega = 2.0 * PI * freq;
69        let c_air = 343.0;
70        (omega / c_air) * (10.3 * x.powf(-0.59))
71    }
72    /// Normal incidence absorption coefficient (approximate).
73    ///
74    /// Uses the impedance at the face of the absorber layer (rigid backing).
75    /// α = 1 − |R|²
76    pub fn absorption_coefficient(&self, freq: f64) -> f64 {
77        let z_r = self.impedance_real(freq);
78        let z_i = self.impedance_imag(freq);
79        let ki = self.propagation_imag(freq);
80        let kd = ki * self.thickness;
81        let coth_kd = if kd > 50.0 { 1.0 } else { 1.0 / kd.tanh() };
82        let zs_r = z_r * coth_kd;
83        let zs_i = z_i * coth_kd;
84        let num_re = zs_r - 1.0;
85        let denom_re = zs_r + 1.0;
86        let r_sq = (num_re * num_re + zs_i * zs_i) / (denom_re * denom_re + zs_i * zs_i);
87        (1.0 - r_sq).clamp(0.0, 1.0)
88    }
89}
90/// Locally resonant acoustic metamaterial (LRAM) unit cell model.
91///
92/// A heavy mass m coated with a soft layer in a matrix medium.
93/// Near the local resonance frequency f_0 = (1/2π) * sqrt(k/m),
94/// the effective mass density becomes negative.
95///
96/// Reference: Liu et al., Science 289, 1734 (2000).
97#[allow(dead_code)]
98#[derive(Debug, Clone)]
99pub struct LocallyResonantUnit {
100    /// Mass of the heavy inclusion (kg).
101    pub mass: f64,
102    /// Stiffness of the soft coating (N/m).
103    pub stiffness: f64,
104    /// Volume of the unit cell (m³).
105    pub cell_volume: f64,
106    /// Background medium density (kg/m³).
107    pub rho_matrix: f64,
108}
109impl LocallyResonantUnit {
110    /// Create a new locally resonant unit cell.
111    pub fn new(mass: f64, stiffness: f64, cell_volume: f64, rho_matrix: f64) -> Self {
112        Self {
113            mass,
114            stiffness,
115            cell_volume,
116            rho_matrix,
117        }
118    }
119    /// Local resonance frequency (Hz).
120    pub fn resonance_frequency(&self) -> f64 {
121        (1.0 / (2.0 * PI)) * (self.stiffness / self.mass).sqrt()
122    }
123    /// Effective mass density at frequency f (Hz).
124    ///
125    /// rho_eff = rho_matrix + (mass / V_cell) * (omega_0² / (omega_0² - omega²))
126    pub fn effective_density(&self, freq: f64) -> f64 {
127        let omega_0 = 2.0 * PI * self.resonance_frequency();
128        let omega = 2.0 * PI * freq;
129        let omega_0_sq = omega_0 * omega_0;
130        let omega_sq = omega * omega;
131        let denom = omega_0_sq - omega_sq;
132        if denom.abs() < 1e-30 {
133            return f64::INFINITY;
134        }
135        let rho_incl = self.mass / self.cell_volume;
136        self.rho_matrix + rho_incl * omega_0_sq / denom
137    }
138    /// Whether the effective density is negative at frequency f.
139    pub fn is_negative_mass(&self, freq: f64) -> bool {
140        let rho_eff = self.effective_density(freq);
141        rho_eff < 0.0
142    }
143}
144/// Frequency-dependent absorption model using multiple relaxation mechanisms.
145///
146/// Models absorption in water including classical viscous losses and
147/// a magnesium sulfate (MgSO₄) relaxation term.
148///
149/// Reference: Francois & Garrison (1982), simplified.
150#[allow(dead_code)]
151#[derive(Debug, Clone)]
152pub struct WaterAbsorption {
153    /// Temperature \[°C\]
154    pub temperature_c: f64,
155    /// Salinity \[ppt (parts per thousand)\]
156    pub salinity_ppt: f64,
157    /// Depth \[m\]
158    pub depth_m: f64,
159}
160impl WaterAbsorption {
161    /// Create a new WaterAbsorption model.
162    #[allow(dead_code)]
163    pub fn new(temperature_c: f64, salinity_ppt: f64, depth_m: f64) -> Self {
164        Self {
165            temperature_c,
166            salinity_ppt,
167            depth_m,
168        }
169    }
170    /// Freshwater (zero salinity) at 20°C, sea surface.
171    #[allow(dead_code)]
172    pub fn freshwater() -> Self {
173        Self::new(20.0, 0.0, 0.0)
174    }
175    /// Seawater at typical ocean conditions (15°C, 35 ppt, 0 m).
176    #[allow(dead_code)]
177    pub fn seawater() -> Self {
178        Self::new(15.0, 35.0, 0.0)
179    }
180    /// Absorption coefficient \[dB/km\] at frequency f \[kHz\].
181    ///
182    /// Uses the Francois-Garrison formula (simplified).
183    #[allow(dead_code)]
184    pub fn absorption_db_per_km(&self, freq_khz: f64) -> f64 {
185        let t = self.temperature_c;
186        let s = self.salinity_ppt;
187        let d = self.depth_m / 1000.0;
188        let f = freq_khz;
189        let f2 = f * f;
190        let f1 = 0.78 * (s / 35.0).sqrt() * (t / 26.0).exp();
191        let a1 = 0.106 * f1 * f2 / (f1 * f1 + f2) * (-d / 6.0).exp();
192        let f2_rel = 42.0 * (t / 17.0).exp();
193        let a2 = 0.52 * (s / 35.0) * (1.0 + t / 40.0) * f2_rel * f2 / (f2_rel * f2_rel + f2)
194            * (-d / 6.0).exp();
195        let a3 = 4.9e-4 * f2 * (-t / 27.0 + d / 17.0).exp();
196        a1 + a2 + a3
197    }
198}
199/// Material properties relevant to acoustic wave propagation.
200#[allow(dead_code)]
201#[derive(Debug, Clone)]
202pub struct AcousticMaterial {
203    /// Mass density in kg/m³
204    pub density: f64,
205    /// Bulk modulus in Pa
206    pub bulk_modulus: f64,
207    /// Shear modulus in Pa (0 for fluids)
208    pub shear_modulus: f64,
209    /// Structural loss factor (dimensionless)
210    pub loss_factor: f64,
211}
212impl AcousticMaterial {
213    /// Create a new acoustic material.
214    #[allow(dead_code)]
215    pub fn new(density: f64, bulk_modulus: f64, shear_modulus: f64, loss_factor: f64) -> Self {
216        Self {
217            density,
218            bulk_modulus,
219            shear_modulus,
220            loss_factor,
221        }
222    }
223    /// Water at 20°C.
224    #[allow(dead_code)]
225    pub fn water() -> Self {
226        Self {
227            density: 998.2,
228            bulk_modulus: 2.18e9,
229            shear_modulus: 0.0,
230            loss_factor: 0.0,
231        }
232    }
233    /// Steel (structural).
234    #[allow(dead_code)]
235    pub fn steel() -> Self {
236        Self {
237            density: 7850.0,
238            bulk_modulus: 166.67e9,
239            shear_modulus: 80.0e9,
240            loss_factor: 0.001,
241        }
242    }
243    /// Air at 20°C, 1 atm.
244    #[allow(dead_code)]
245    pub fn air() -> Self {
246        Self {
247            density: 1.204,
248            bulk_modulus: 141.9e3,
249            shear_modulus: 0.0,
250            loss_factor: 0.0,
251        }
252    }
253    /// Concrete (normal weight).
254    #[allow(dead_code)]
255    pub fn concrete() -> Self {
256        Self {
257            density: 2300.0,
258            bulk_modulus: 13.33e9,
259            shear_modulus: 10.0e9,
260            loss_factor: 0.02,
261        }
262    }
263    /// Longitudinal (compressional) wave velocity: c_L = sqrt((K + 4G/3) / rho).
264    #[allow(dead_code)]
265    pub fn longitudinal_velocity(&self) -> f64 {
266        let m_modulus = self.bulk_modulus + 4.0 * self.shear_modulus / 3.0;
267        (m_modulus / self.density).sqrt()
268    }
269    /// Shear wave velocity: c_S = sqrt(G / rho). Returns 0 for fluids.
270    #[allow(dead_code)]
271    pub fn shear_velocity(&self) -> f64 {
272        if self.shear_modulus == 0.0 {
273            0.0
274        } else {
275            (self.shear_modulus / self.density).sqrt()
276        }
277    }
278    /// Characteristic acoustic impedance: Z = rho * c_L (Pa·s/m).
279    #[allow(dead_code)]
280    pub fn acoustic_impedance(&self) -> f64 {
281        self.density * self.longitudinal_velocity()
282    }
283    /// Compute the insertion loss (IL) of a barrier made of this material.
284    ///
285    /// The insertion loss is the reduction in sound power level (SPL) achieved
286    /// by inserting a barrier between a source and receiver. Using the
287    /// simplified mass-law / Maekawa formula for a thin barrier:
288    ///
289    /// IL = 20 · log10(1 + π · N · f · (m_s / (ρ_air · c_air)) )
290    ///
291    /// where N = Fresnel number characterising diffraction over the barrier,
292    /// and m_s = ρ · h is the surface mass density of the barrier panel.
293    ///
294    /// For thin panels the dominant mechanism is mass-law transmission loss,
295    /// so we combine barrier TL with Maekawa diffraction:
296    ///
297    /// IL ≈ TL_mass + ΔL_diffraction
298    ///
299    /// Here we use the practical formula:
300    ///
301    /// IL = 10 · log10( 1 + 20 · N )   \[dB\]  (Maekawa 1968, simplified)
302    ///
303    /// # Arguments
304    /// * `frequency`     - Frequency f \[Hz\]
305    /// * `thickness`     - Barrier panel thickness h \[m\]
306    /// * `fresnel_number` - Fresnel number N = 2·δ/λ where δ is path-length
307    ///   difference (dimensionless, N ≥ 0)
308    ///
309    /// # Returns
310    /// Insertion loss IL \[dB\]
311    #[allow(dead_code)]
312    pub fn compute_insertion_loss(
313        &self,
314        frequency: f64,
315        thickness: f64,
316        fresnel_number: f64,
317    ) -> f64 {
318        let m_s = self.density * thickness;
319        let rho_air = 1.204_f64;
320        let c_air = 343.0_f64;
321        let tl_mass = 20.0 * (1.0 + PI * frequency * m_s / (rho_air * c_air)).log10();
322        let il_diffraction = if fresnel_number > 0.0 {
323            10.0 * (1.0 + 20.0 * fresnel_number).log10()
324        } else {
325            0.0
326        };
327        tl_mass + il_diffraction
328    }
329    /// Compute the Noise Reduction Coefficient (NRC) for this material.
330    ///
331    /// The NRC is a single-number rating of sound absorption, defined as the
332    /// arithmetic mean of the Sabine absorption coefficients at 250, 500,
333    /// 1000, and 2000 Hz, rounded to the nearest 0.05.
334    ///
335    /// For a material characterised by its loss factor η, the statistical
336    /// absorption coefficient at each frequency band is estimated from the
337    /// relationship between loss factor and random-incidence absorption
338    /// coefficient:
339    ///
340    /// α_stat(f) ≈ η(f) · m_s · ω / (ρ_air · c_air)
341    ///
342    /// where m_s = ρ · h is the surface mass density and ω = 2πf. When η
343    /// is independent of frequency (as stored in `self.loss_factor`), the
344    /// coefficients scale with f, so we compute and average them.
345    ///
346    /// The result is clamped to \[0, 1\].
347    ///
348    /// # Arguments
349    /// * `thickness` - Panel thickness h \[m\]
350    ///
351    /// # Returns
352    /// NRC \[dimensionless, 0–1\]
353    #[allow(dead_code)]
354    pub fn compute_noise_reduction_coefficient(&self, thickness: f64) -> f64 {
355        let bands = [250.0_f64, 500.0, 1000.0, 2000.0];
356        let rho_air = 1.204_f64;
357        let c_air = 343.0_f64;
358        let m_s = self.density * thickness;
359        let sum: f64 = bands
360            .iter()
361            .map(|&f| {
362                let omega = 2.0 * PI * f;
363                let alpha = self.loss_factor * m_s * omega / (rho_air * c_air);
364                alpha.clamp(0.0, 1.0)
365            })
366            .sum();
367        let nrc_raw = sum / bands.len() as f64;
368        let nrc_rounded = (nrc_raw / 0.05).round() * 0.05;
369        nrc_rounded.clamp(0.0, 1.0)
370    }
371    /// Compute the transmission loss (TL) by the mass law.
372    ///
373    /// The mass law is the classical result for a limp, homogeneous panel:
374    ///
375    /// TL = 20 · log10(m_s · f) − 47   \[dB\]  (for normal incidence, SI units)
376    ///
377    /// or equivalently for random-incidence (field incidence):
378    ///
379    /// TL = 20 · log10(m_s · f) − 47.3  \[dB\]
380    ///
381    /// where m_s = ρ · h \[kg/m²\] is the surface mass density.
382    /// The constant 47 arises from 20·log10(2π / (ρ_air · c_air)).
383    ///
384    /// At double the frequency or double the mass, TL increases by ~6 dB
385    /// ("mass law doubling rule").
386    ///
387    /// # Arguments
388    /// * `frequency` - Frequency f \[Hz\]
389    /// * `thickness` - Panel thickness h \[m\]
390    ///
391    /// # Returns
392    /// Transmission loss TL \[dB\]
393    #[allow(dead_code)]
394    pub fn compute_transmission_loss_mass_law(&self, frequency: f64, thickness: f64) -> f64 {
395        let m_s = self.density * thickness;
396        let rho_air = 1.204_f64;
397        let c_air = 343.0_f64;
398        let arg = PI * m_s * frequency / (rho_air * c_air);
399        20.0 * arg.log10()
400    }
401}
402/// Rubber (natural rubber, acoustic grade).
403#[allow(dead_code)]
404impl AcousticMaterial {
405    /// Natural rubber (acoustic grade).
406    pub fn rubber() -> Self {
407        Self {
408            density: 1200.0,
409            bulk_modulus: 2.0e9,
410            shear_modulus: 0.5e6,
411            loss_factor: 0.1,
412        }
413    }
414    /// Aluminium alloy (2024-T3).
415    pub fn aluminium() -> Self {
416        Self {
417            density: 2780.0,
418            bulk_modulus: 73.1e9,
419            shear_modulus: 27.0e9,
420            loss_factor: 0.0002,
421        }
422    }
423    /// Glass (borosilicate).
424    pub fn glass() -> Self {
425        Self {
426            density: 2230.0,
427            bulk_modulus: 46.0e9,
428            shear_modulus: 26.0e9,
429            loss_factor: 0.0005,
430        }
431    }
432    /// Sound speed through the material \[m/s\] (alias for longitudinal_velocity).
433    #[allow(dead_code)]
434    pub fn sound_speed(&self) -> f64 {
435        self.longitudinal_velocity()
436    }
437    /// Specific acoustic impedance Z = ρ·c normalised to air impedance.
438    ///
439    /// Z_air = 413 Pa·s/m (20°C, 1 atm).
440    #[allow(dead_code)]
441    pub fn relative_impedance(&self) -> f64 {
442        self.acoustic_impedance() / 413.0
443    }
444}
445/// Ultrasonic non-destructive evaluation (NDE) signal model.
446///
447/// Simulates pulse-echo time-of-flight and amplitude for a pulse reflected
448/// from a planar reflector.
449#[derive(Debug, Clone)]
450#[allow(dead_code)]
451pub struct UltrasonicNDE {
452    /// Wave velocity in the test material \[m/s\]
453    pub velocity: f64,
454    /// Ultrasonic frequency \[Hz\]
455    pub frequency: f64,
456    /// Initial pulse amplitude \[Pa\] (peak pressure)
457    pub pulse_amplitude: f64,
458    /// Attenuation coefficient \[dB/m\]
459    pub attenuation_db_m: f64,
460}
461impl UltrasonicNDE {
462    /// Create a new `UltrasonicNDE` model.
463    #[allow(dead_code)]
464    pub fn new(velocity: f64, frequency: f64, pulse_amplitude: f64, attenuation_db_m: f64) -> Self {
465        Self {
466            velocity,
467            frequency,
468            pulse_amplitude,
469            attenuation_db_m,
470        }
471    }
472    /// Round-trip time-of-flight \[s\] for a reflector at depth `depth_m` \[m\].
473    #[allow(dead_code)]
474    pub fn time_of_flight(&self, depth_m: f64) -> f64 {
475        2.0 * depth_m / self.velocity
476    }
477    /// Echo amplitude \[Pa\] after round-trip path of `2·depth` \[m\].
478    ///
479    /// Accounts for beam spreading (1/r law) and material attenuation.
480    #[allow(dead_code)]
481    pub fn echo_amplitude(&self, depth_m: f64, reflection_coeff: f64) -> f64 {
482        if depth_m <= 0.0 {
483            return self.pulse_amplitude * reflection_coeff.abs();
484        }
485        let d = 2.0 * depth_m;
486        let geo = 1.0 / (1.0 + d);
487        let alpha_np = self.attenuation_db_m / 8.686;
488        let att = (-alpha_np * d).exp();
489        self.pulse_amplitude * reflection_coeff.abs() * geo * att
490    }
491    /// Wavelength in the material \[m\].
492    #[allow(dead_code)]
493    pub fn wavelength(&self) -> f64 {
494        self.velocity / self.frequency
495    }
496    /// Near-field (Fresnel) distance \[m\] for a circular transducer of radius `r` \[m\].
497    ///
498    /// `N = r² / λ`
499    #[allow(dead_code)]
500    pub fn near_field_distance(&self, transducer_radius: f64) -> f64 {
501        transducer_radius * transducer_radius / self.wavelength()
502    }
503    /// Lateral resolution at depth `d` \[m\] (Rayleigh criterion for focused beam):
504    ///
505    /// `x_r ≈ 1.22 · λ · d / (2·r)`
506    #[allow(dead_code)]
507    pub fn lateral_resolution(&self, depth_m: f64, transducer_radius: f64) -> f64 {
508        1.22 * self.wavelength() * depth_m / (2.0 * transducer_radius)
509    }
510}
511/// Simplified 1D phononic crystal band gap estimator.
512///
513/// For a 1D bi-material phononic crystal with alternating layers A and B:
514///
515/// Band gap frequency: f_gap = c_eff / (2 * a)  where a is the lattice constant.
516/// The gap width depends on the impedance mismatch.
517///
518/// This uses a simplified formula from Sigalas & Economou (1992):
519/// Δω/ω_0 ≈ (2/π) * |arcsin((Z_b - Z_a) / (Z_b + Z_a))|
520#[allow(dead_code)]
521#[derive(Debug, Clone)]
522pub struct PhononicCrystal1D {
523    /// Material A impedance (Pa·s/m).
524    pub z_a: f64,
525    /// Material B impedance (Pa·s/m).
526    pub z_b: f64,
527    /// Wave speed in material A (m/s).
528    pub c_a: f64,
529    /// Wave speed in material B (m/s).
530    pub c_b: f64,
531    /// Thickness of material A layer (m).
532    pub d_a: f64,
533    /// Thickness of material B layer (m).
534    pub d_b: f64,
535}
536impl PhononicCrystal1D {
537    /// Create a new 1D phononic crystal.
538    pub fn new(z_a: f64, z_b: f64, c_a: f64, c_b: f64, d_a: f64, d_b: f64) -> Self {
539        Self {
540            z_a,
541            z_b,
542            c_a,
543            c_b,
544            d_a,
545            d_b,
546        }
547    }
548    /// Lattice constant a = d_a + d_b.
549    pub fn lattice_constant(&self) -> f64 {
550        self.d_a + self.d_b
551    }
552    /// First band gap centre frequency (Hz).
553    ///
554    /// Based on the Bragg condition: 2a = λ → f_Bragg = c_eff / (2a).
555    pub fn bragg_frequency(&self) -> f64 {
556        let c_eff = (self.c_a * self.d_a + self.c_b * self.d_b) / self.lattice_constant();
557        c_eff / (2.0 * self.lattice_constant())
558    }
559    /// Estimated relative band gap width Δf/f_0.
560    pub fn relative_gap_width(&self) -> f64 {
561        let r = (self.z_b - self.z_a).abs() / (self.z_a + self.z_b);
562        (2.0 / PI) * r.asin()
563    }
564    /// Absolute band gap width (Hz).
565    pub fn gap_width_hz(&self) -> f64 {
566        self.relative_gap_width() * self.bragg_frequency()
567    }
568    /// Whether a given frequency lies in the first band gap.
569    pub fn is_in_gap(&self, freq: f64) -> bool {
570        let f0 = self.bragg_frequency();
571        let half_gap = 0.5 * self.gap_width_hz();
572        (freq - f0).abs() < half_gap
573    }
574}