diff --git a/.github/images/silencer_expansion_chamber.svg b/.github/images/silencer_expansion_chamber.svg
new file mode 100644
index 00000000..d922ee6c
--- /dev/null
+++ b/.github/images/silencer_expansion_chamber.svg
@@ -0,0 +1,1278 @@
+
+
+
+
+
+
+ image/svg+xml
+
+
+ Matplotlib v3.11.0, https://matplotlib.org/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 31.5
+
+
+
+
+
+
+
+
+
+
+ 63
+
+
+
+
+
+
+
+
+
+
+ 125
+
+
+
+
+
+
+
+
+
+
+ 250
+
+
+
+
+
+
+
+
+
+
+ 500
+
+
+
+
+
+
+
+
+
+
+ 1k
+
+
+
+
+
+
+
+
+
+
+ 2k
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Frequency [Hz]
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 0.0
+
+
+
+
+
+
+
+
+
+
+ 2.5
+
+
+
+
+
+
+
+
+
+
+ 5.0
+
+
+
+
+
+
+
+
+
+
+ 7.5
+
+
+
+
+
+
+
+
+
+
+ 10.0
+
+
+
+
+
+
+
+
+
+
+ 12.5
+
+
+
+
+
+
+
+
+
+
+ 15.0
+
+
+
+
+
+
+
+
+
+
+ 17.5
+
+
+
+
+
+
+
+
+
+
+ 20.0
+
+
+
+ Transmission loss [dB]
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Expansion-chamber transmission loss (Bies Eq. 8.111)
+
+
+
+
+
+ Area ratio m = Sexp/Sduct
+
+
+
+
+ m = 2 → 1.9 dB
+
+
+
+
+ m = 4 → 6.5 dB
+
+
+
+
+ m = 8 → 12.2 dB
+
+
+
+
+ m = 16 → 18.1 dB
+
+
+
+
+
+
+
+
+
diff --git a/.github/images/silencer_expansion_chamber_dark.svg b/.github/images/silencer_expansion_chamber_dark.svg
new file mode 100644
index 00000000..0bb18367
--- /dev/null
+++ b/.github/images/silencer_expansion_chamber_dark.svg
@@ -0,0 +1,1278 @@
+
+
+
+
+
+
+ image/svg+xml
+
+
+ Matplotlib v3.11.0, https://matplotlib.org/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 31.5
+
+
+
+
+
+
+
+
+
+
+ 63
+
+
+
+
+
+
+
+
+
+
+ 125
+
+
+
+
+
+
+
+
+
+
+ 250
+
+
+
+
+
+
+
+
+
+
+ 500
+
+
+
+
+
+
+
+
+
+
+ 1k
+
+
+
+
+
+
+
+
+
+
+ 2k
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Frequency [Hz]
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 0.0
+
+
+
+
+
+
+
+
+
+
+ 2.5
+
+
+
+
+
+
+
+
+
+
+ 5.0
+
+
+
+
+
+
+
+
+
+
+ 7.5
+
+
+
+
+
+
+
+
+
+
+ 10.0
+
+
+
+
+
+
+
+
+
+
+ 12.5
+
+
+
+
+
+
+
+
+
+
+ 15.0
+
+
+
+
+
+
+
+
+
+
+ 17.5
+
+
+
+
+
+
+
+
+
+
+ 20.0
+
+
+
+ Transmission loss [dB]
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Expansion-chamber transmission loss (Bies Eq. 8.111)
+
+
+
+
+
+ Area ratio m = Sexp/Sduct
+
+
+
+
+ m = 2 → 1.9 dB
+
+
+
+
+ m = 4 → 6.5 dB
+
+
+
+
+ m = 8 → 12.2 dB
+
+
+
+
+ m = 16 → 18.1 dB
+
+
+
+
+
+
+
+
+
diff --git a/.github/images/silencer_expansion_chamber_es.svg b/.github/images/silencer_expansion_chamber_es.svg
new file mode 100644
index 00000000..147741c7
--- /dev/null
+++ b/.github/images/silencer_expansion_chamber_es.svg
@@ -0,0 +1,1278 @@
+
+
+
+
+
+
+ image/svg+xml
+
+
+ Matplotlib v3.11.0, https://matplotlib.org/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 31,5
+
+
+
+
+
+
+
+
+
+
+ 63
+
+
+
+
+
+
+
+
+
+
+ 125
+
+
+
+
+
+
+
+
+
+
+ 250
+
+
+
+
+
+
+
+
+
+
+ 500
+
+
+
+
+
+
+
+
+
+
+ 1k
+
+
+
+
+
+
+
+
+
+
+ 2k
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Frecuencia [Hz]
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 0
+
+
+
+
+
+
+
+
+
+
+ 2,5
+
+
+
+
+
+
+
+
+
+
+ 5
+
+
+
+
+
+
+
+
+
+
+ 7,5
+
+
+
+
+
+
+
+
+
+
+ 10
+
+
+
+
+
+
+
+
+
+
+ 12,5
+
+
+
+
+
+
+
+
+
+
+ 15
+
+
+
+
+
+
+
+
+
+
+ 17,5
+
+
+
+
+
+
+
+
+
+
+ 20
+
+
+
+ Pérdida por transmisión [dB]
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Pérdida de transmisión de cámara de expansión (Bies Ec. 8,111)
+
+
+
+
+
+ Relación de áreas m = Sexp/Sduct
+
+
+
+
+ m = 2 → 1,9 dB
+
+
+
+
+ m = 4 → 6,5 dB
+
+
+
+
+ m = 8 → 12,2 dB
+
+
+
+
+ m = 16 → 18,1 dB
+
+
+
+
+
+
+
+
+
diff --git a/.github/images/silencer_expansion_chamber_es_dark.svg b/.github/images/silencer_expansion_chamber_es_dark.svg
new file mode 100644
index 00000000..c0431285
--- /dev/null
+++ b/.github/images/silencer_expansion_chamber_es_dark.svg
@@ -0,0 +1,1278 @@
+
+
+
+
+
+
+ image/svg+xml
+
+
+ Matplotlib v3.11.0, https://matplotlib.org/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 31,5
+
+
+
+
+
+
+
+
+
+
+ 63
+
+
+
+
+
+
+
+
+
+
+ 125
+
+
+
+
+
+
+
+
+
+
+ 250
+
+
+
+
+
+
+
+
+
+
+ 500
+
+
+
+
+
+
+
+
+
+
+ 1k
+
+
+
+
+
+
+
+
+
+
+ 2k
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Frecuencia [Hz]
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 0
+
+
+
+
+
+
+
+
+
+
+ 2,5
+
+
+
+
+
+
+
+
+
+
+ 5
+
+
+
+
+
+
+
+
+
+
+ 7,5
+
+
+
+
+
+
+
+
+
+
+ 10
+
+
+
+
+
+
+
+
+
+
+ 12,5
+
+
+
+
+
+
+
+
+
+
+ 15
+
+
+
+
+
+
+
+
+
+
+ 17,5
+
+
+
+
+
+
+
+
+
+
+ 20
+
+
+
+ Pérdida por transmisión [dB]
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Pérdida de transmisión de cámara de expansión (Bies Ec. 8,111)
+
+
+
+
+
+ Relación de áreas m = Sexp/Sduct
+
+
+
+
+ m = 2 → 1,9 dB
+
+
+
+
+ m = 4 → 6,5 dB
+
+
+
+
+ m = 8 → 12,2 dB
+
+
+
+
+ m = 16 → 18,1 dB
+
+
+
+
+
+
+
+
+
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 434342c6..820dd566 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -67,6 +67,39 @@ and this project adheres to [Semantic Versioning](https://semver.org/).
`W = 0.5 |F|^2 Re{Y}` (Cremer, Heckl & Petersson 2005, *Structure-Borne Sound*
3e, Table 5.1), the theoretical companions of the measured ISO 7626 mobilities.
+- Industrial noise control: silencers, HVAC and machine enclosures, a new
+ `phonometry.noise_control` domain (Bies, Hansen & Howard, *Engineering Noise
+ Control* 5th ed.; Munjal, *Acoustics of Ducts and Mufflers*). `expansion_chamber`,
+ `helmholtz_resonator`, `quarter_wave_resonator` and `extended_tube_chamber`
+ build **reactive silencers** by the one-dimensional four-pole
+ (transmission-matrix) method, returning a `ReactiveSilencerResult` with the
+ transmission loss, the insertion loss for configurable source/radiation
+ impedances, the compound transfer matrix and `.plot()`; the expansion-chamber
+ transmission loss matches the closed form `10 lg[1 + (1/4)(m - 1/m)^2 sin^2 kL]`
+ (Bies Eq. (8.111)) exactly, and the four-pole machinery is available under
+ `phonometry.noise_control` (`duct_matrix`, `shunt_matrix`, `cascade`,
+ `transmission_loss`, `insertion_loss`, `helmholtz_impedance`,
+ `quarter_wave_impedance`). The
+ `noise_control.hvac` module adds the duct end reflection (`end_reflection_loss`,
+ ASHRAE Table 8.14), bend insertion loss (`elbow_insertion_loss`, Table 8.11),
+ the plenum transmission loss by Wells' method (`plenum_attenuation`) and the
+ flow-generated (self) noise of straight ducts and mitred bends
+ (`flow_noise_straight_duct`, `flow_noise_bend`, VDI 2081). `enclosure_insertion_loss`
+ gives the machine-enclosure insertion loss `IL = R - C` (Bies Eqs. (7.103),
+ (7.111)) from a **caller-supplied** panel transmission loss `R` (array or
+ callable; never predicted here) and the interior room constant reused from
+ `phonometry.room.room_constant`. The four-pole expansion chamber is
+ cross-checked against the independent 2D FDTD solver.
+- Radiation of a rigid circular piston in an infinite baffle
+ (`phonometry.electroacoustics.piston`, re-exported at the top level). `radiating_piston`
+ computes the radiation impedance `rho c S (R1 + j X1)` from the piston
+ resistance `R1 = 1 - 2 J1(2ka)/2ka` and reactance `X1 = 2 H1(2ka)/2ka`
+ (Beranek & Mellow, *Acoustics* 2nd ed., Eqs. (13.117), (13.118)), the
+ low-frequency radiation mass `8 rho a^3 / 3` (Eq. (4.151)), the far-field
+ directivity `2 J1(ka sin theta)/(ka sin theta)` and the directivity index,
+ returning a `RadiatingPistonResult` with `.plot()`; the building blocks
+ `piston_resistance`, `piston_reactance` and `piston_directivity` are exposed
+ directly.
- Objective speech intelligibility STOI and ESTOI
(`phonometry.hearing.objective_intelligibility.stoi`, also `phonometry.stoi`).
`stoi(clean, degraded, fs)` computes the short-time objective intelligibility
diff --git a/docs/CONFORMANCE.md b/docs/CONFORMANCE.md
index 836c2784..e535f374 100644
--- a/docs/CONFORMANCE.md
+++ b/docs/CONFORMANCE.md
@@ -15,7 +15,7 @@
## Numerical conformance report
-✅ **338/338 conformance checks pass** across 42 domains and 210 standards - filters class 1 - weightings within IEC 61672-1 class 1.
+✅ **353/353 conformance checks pass** across 44 domains and 224 standards - filters class 1 - weightings within IEC 61672-1 class 1.
Each row pins a standard clause to its expected normative value and the value the library computes. Every section below is collapsible and stays collapsed while all of its rows pass; a section with any failing row opens automatically.
@@ -719,3 +719,34 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli
+
+✅ Electroacoustics : 100% (6/6)
+
+| Standard | Quantity | Expected (norm) | Computed | Δ | Status |
+|:---|:---|:---|:---|:---|:---:|
+| Beranek & Mellow 2e Eq. (13.117) | Piston resistance R1(x) = 1 - 2 J1(x)/x at x = 2ka = 2 | 0.423275 (+/-0.00001) | 0.423275 | 0 | ✅ |
+| Beranek & Mellow 2e Eq. (13.118) | Piston reactance X1(x) = 2 H1(x)/x at x = 2ka = 2 | 0.646764 (+/-0.00001) | 0.646764 | 0 | ✅ |
+| Beranek & Mellow 2e Eq. (13.117) (low-frequency limit) | R1 -> (ka)^2/2 as ka -> 0 (x = 0.02, ka = 0.01) | 0.00005 (+/-0.01%) | 0.00005 | 0 | ✅ |
+| Beranek & Mellow 2e Eq. (4.151) | Radiation mass M = 8 rho a^3 / 3 (a = 0.1 m, rho = 1.206) | 0.003216 kg (+/-0 kg) | 0.003216 kg | 0 kg | ✅ |
+| Beranek & Mellow 2e Eq. (13.102), Table 14.1 | First directivity null at ka sin(theta) = 3.8317 (first zero of J1) | 0 (+/-0.000001) | 0 | 0 | ✅ |
+| Beranek & Mellow 2e §4.19 (half-space baffle) | Directivity index DI -> 10 lg 2 = 3.01 dB as ka -> 0 | 3.0103 dB (+/-0.001 dB) | 3.0103 dB | 0 dB | ✅ |
+
+
+
+
+✅ Industrial noise control : 100% (9/9)
+
+| Standard | Quantity | Expected (norm) | Computed | Δ | Status |
+|:---|:---|:---|:---|:---|:---:|
+| Bies 5e Eq. (8.111) | Expansion-chamber peak TL = 10 lg[1 + (1/4)(m - 1/m)^2], m = 4 at kL = pi/2 | 6.5472 dB (+/-0 dB) | 6.5472 dB | 0 dB | ✅ |
+| Bies 5e Eq. (8.111) | Expansion-chamber trough TL = 0 at kL = pi (chamber transparent) | 0 dB (+/-0 dB) | 0 dB | 0 dB | ✅ |
+| Bies 5e Eq. (8.44) / Example 8.1 | Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz | 56.6 Hz (+/-0.1 Hz) | 56.6 Hz | 0.003 Hz | ✅ |
+| Bies 5e Eq. (8.46) | Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V)) (S=1e-4, l_e=0.02, V=1e-3) | 122.067 Hz (+/-0 Hz) | 122.067 Hz | 0 Hz | ✅ |
+| Bies 5e Eq. (8.73) | Side-branch TL = 20 lg abs(1 + rho c/(2 Sd Zb)) (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | ✅ |
+| Bies 5e Eqs. (8.141)/(8.148) (four-pole insertion loss) | Insertion loss = transmission loss for the anechoic reference Zs=Zr=rho c/S | 6.2498 dB (= TL) | 6.2498 dB | 0 dB | ✅ |
+| Bies 5e Eq. (8.275) (Wells' plenum method) | Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2) | 12.8541 dB (+/-0 dB) | 12.8541 dB | 0 dB | ✅ |
+| Bies 5e Table 8.14 (ASHRAE end reflection, flush) | Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node) | 10 dB (+/-0 dB) | 10 dB | 0 dB | ✅ |
+| Bies 5e Eqs. (7.103), (7.111) (enclosure, fully absorbing limit) | Enclosure correction C -> 10 lg 0.3 = -5.23 dB as alpha_i -> 1 | -5.2288 dB (+/-0.001 dB) | -5.2288 dB | 0 dB | ✅ |
+
+
+
diff --git a/docs/README.md b/docs/README.md
index 6b1f7ab0..f75b38e1 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -23,6 +23,7 @@ Full documentation for phonometry. Also available as a website:
- [Objective intelligibility (STOI & ESTOI)](objective-intelligibility.md): the correlation-based intelligibility measures for time-frequency weighted noisy speech from a clean/degraded pair: STOI (Taal et al. 2011), the clipped per-band envelope correlation, and ESTOI (Jensen & Taal 2016), the row- and column-normalised spectral correlation that tracks modulated maskers
- [Electroacoustics: distortion & frequency response](electroacoustics.md): the IEC 60268-3 distortion set (THD, nth-order harmonic, THD+N and SINAD via AES17, SMPTE and CCIF intermodulation, dynamic intermodulation and weighted THD), the AES17 dynamic range and idle channel noise, and the Bendat & Piersol H1/H2 frequency-response estimators with the ordinary coherence γ²
- [Swept-sine distortion and phase utilities](swept-sine-distortion.md): harmonic separation from one exponential sweep (Farina 2000) with the synchronized swept-sine of Novak et al. 2015 for coherent harmonic phases, THD as a function of the excitation frequency, and minimum phase from |H| (real cepstrum), group delay and excess phase
+- [Industrial noise control: silencers, HVAC & enclosures](noise-control.md): reactive silencers by the four-pole transmission-matrix method (expansion chambers, Helmholtz and quarter-wave resonators, extended tubes) with transmission and insertion loss, the HVAC duct methods (end reflection, elbows, plenums and flow-generated noise) and machine-enclosure insertion loss from a supplied panel R and the interior room constant (Bies, Hansen & Howard)
- [Programme loudness & true peak](program-loudness.md): the ITU-R BS.1770-5 programme loudness (K-weighting, gated 400 ms blocks, channel weights including the Annex 3 positions) and the oversampled true-peak level in dBTP, with the EBU R 128 −23 LUFS practice, the Tech 3341 EBU Mode momentary/short-term/integrated meters and the Tech 3342 loudness range, validated against the official EBU test signals
- [Underwater acoustics: radiated noise & pile driving](underwater-acoustics.md): the ISO 18405 reference levels (SPL, SEL, peak re 1 µPa), the ISO 17208 ship radiated noise level and equivalent monopole source level via the Lloyd's-mirror correction, and the ISO 18406 single-strike, peak and cumulative pile-driving sound exposure
- [Underwater sound propagation](underwater-propagation.md): closed-form transmission loss (geometrical spreading plus volume absorption by Francois-Garrison, Ainslie-McColm or Thorp), the speed of sound in sea water (UNESCO/Chen-Millero, Del Grosso, Mackenzie) with the sound-speed profile, the passive/active sonar equation, seabed reflection loss (Rayleigh), the ocean ambient-noise spectrum (Wenz wind/thermal plus JOMOPANS-ECHO ship-traffic source levels) and numerical solvers (normal modes, ray tracing, parabolic equation)
diff --git a/docs/api-reference.md b/docs/api-reference.md
index dca8129d..57bf7922 100644
--- a/docs/api-reference.md
+++ b/docs/api-reference.md
@@ -356,6 +356,20 @@ cycle and warn on use; they are removed in 4.0.
| `synchronized_sweep_signal` | `function` | **Synchronized exponential sweep (Novak et al. 2015, Eqs. 47/49).**
+Mostrar el código de esta figura
+
+```python
+import matplotlib.pyplot as plt
+import numpy as np
+from phonometry import expansion_chamber
+
+freqs = np.linspace(20.0, 2000.0, 2000)
+
+# Una línea para una cámara: TL vs frecuencia (con la pérdida de inserción si
+# se dan las impedancias de fuente y radiación).
+expansion_chamber(freqs, 0.3, 0.04, 0.01,
+ source_impedance=4e4, radiation_impedance=5e3).plot()
+plt.show()
+
+# A mano: la familia de relaciones de áreas de la figura de concepto.
+fig, ax = plt.subplots()
+for m in (2.0, 4.0, 8.0, 16.0):
+ res = expansion_chamber(freqs, 0.3, m * 0.01, 0.01)
+ ax.plot(freqs, res.transmission_loss, label=f"m = {int(m)}")
+ax.set_xlabel("Frecuencia [Hz]"); ax.set_ylabel("Pérdida de transmisión [dB]")
+ax.legend()
+plt.show()
+```
+
+
+
+### Resonadores laterales y de tubo extendido
+
+Un **resonador de Helmholtz** (área de cuello $S_n$, longitud efectiva $l_e$,
+volumen de cavidad $V$) y un **tubo de cuarto de onda** cerrado (longitud $l_e$)
+cortocircuitan el conducto en su frecuencia de sintonía, dando un pico agudo de
+pérdida de transmisión: $f_0 = \tfrac{c}{2\pi}\sqrt{S_n/(l_e V)}$ (Bies Ec. (8.46))
+y $f = c/4l_e$ (Ec. (8.44)). Una **cámara de tubo extendido** aloja ramas de
+cuarto de onda dentro de una cámara de expansión para rellenar sus caídas.
+
+```python
+import numpy as np
+from phonometry import (
+ helmholtz_resonator, quarter_wave_resonator, extended_tube_chamber,
+)
+
+f = np.linspace(20.0, 600.0, 4000)
+
+hr = helmholtz_resonator(f, duct_area=0.01, neck_area=1e-4,
+ neck_length=0.02, cavity_volume=1e-3)
+print(round(float(hr.resonances[0]), 1)) # frecuencia de sintonía, Hz
+
+qw = quarter_wave_resonator(f, duct_area=0.01, length=1.516, branch_area=2e-3,
+ speed_of_sound=343.24)
+print(round(float(qw.resonances[0]), 1)) # 56.6 Hz (Bies Ejemplo 8.1)
+
+# Una extensión de entrada de L/4 rellena la primera caída de la cámara.
+et = extended_tube_chamber(f, length=0.4, chamber_area=0.04, pipe_area=0.01,
+ inlet_extension=0.1)
+```
+
+Cada uno devuelve un `ReactiveSilencerResult` con `transmission_loss`,
+`insertion_loss` (cuando se dan las impedancias), la matriz compuesta
+`transfer_matrix` y `.plot()`. Las configuraciones avanzadas encadenan elementos
+directamente con `duct_matrix`, `shunt_matrix`, `cascade`, `transmission_loss` e
+`insertion_loss`.
+
+## 2. Atenuación y ruido de flujo en conductos HVAC
+
+Una instalación de ventilación atenúa el ruido del ventilador en los codos, en
+el extremo abierto del conducto y en los plenums, y regenera ruido allí donde el
+flujo se perturba. `phonometry.noise_control.hvac` reúne los métodos del
+Capítulo 8 de Bies.
+
+```python
+from phonometry.noise_control import hvac
+
+bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0]
+
+# Reflexión de baja frecuencia en el extremo abierto (ASHRAE Tabla 8.14).
+er = hvac.end_reflection_loss(bands, diameter=0.30, termination="flush")
+
+# Pérdida de inserción de un codo cuadrado a 90 grados con revestimiento (Tabla 8.11).
+el = hvac.elbow_insertion_loss(bands, width=0.3, bend_type="square", lined=True)
+
+# TL de plenum por el método de Wells (forma cerrada).
+tl = hvac.plenum_attenuation(exit_area=0.1, line_of_sight=1.0,
+ wall_area=20.0, mean_absorption=0.2)
+print(round(tl, 1)) # dB
+
+# Ruido autogenerado (de flujo) de un conducto recto (VDI 2081).
+fn = hvac.flow_noise_straight_duct(bands, flow_velocity=10.0, area=0.04)
+```
+
+Los métodos de reflexión de extremo y de codos interpolan las tablas ASHRAE
+(pasan exactamente por los nodos tabulados); los métodos de plenum (Wells) y de
+ruido de flujo (VDI 2081) son formas cerradas. `end_reflection_loss`,
+`elbow_insertion_loss`, `flow_noise_straight_duct` y `flow_noise_bend` devuelven
+un `HvacSpectrumResult` (atenuación o nivel de potencia sonora regenerada) con
+`.plot()`; `plenum_attenuation` devuelve la pérdida de transmisión directamente.
+Los conductos rectangulares usan el diámetro equivalente $D = \sqrt{4S/\pi}$.
+
+:::note
+Bies 5.ª ed. da la reflexión de extremo solo como la consulta de la Tabla 8.14
+de ASHRAE; no hay forma cerrada en esa edición. Este módulo reproduce e
+interpola la tabla.
+:::
+
+## 3. Cerramientos de máquina
+
+Un cerramiento sellado reduce el ruido radiado en la pérdida de transmisión del
+panel $R$, menos una penalización $C$ por la acumulación reverberante en la
+cavidad pequeña y dura (Bies Ecs. (7.103), (7.111)):
+
+$$
+\mathrm{IL} = R - C,\qquad C = 10\log_{10}\!\left(0.3 + \frac{S_E}{R_i}\right),
+$$
+
+con el área externa $S_E$ y la **constante de sala** interior
+$R_i = S_i \alpha_i/(1-\alpha_i)$ (la misma
+[`room_constant`](/phonometry/es/guides/room-image-sources/) del campo estacionario
+de sala). Un interior duro desperdicia gran parte de $R$; revestirlo lleva $C$
+hacia su límite inferior $10\log_{10}0.3 = -5.2$ dB.
+
+**La pérdida de transmisión del panel $R$ la aportas tú** como un array por banda
+(medido, o predicho por un modelo de panel), o como una función de la frecuencia.
+Este módulo nunca predice $R$; combina un $R$ dado con la absorción interior.
+
+```python
+import numpy as np
+from phonometry import enclosure_insertion_loss
+
+bands = np.array([125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0])
+panel_R = np.array([18.0, 24.0, 30.0, 36.0, 42.0, 46.0]) # medido, dB
+
+enc = enclosure_insertion_loss(
+ panel_R, external_area=6.0, internal_area=5.0,
+ internal_absorption=0.3, frequencies=bands,
+)
+print(np.round(enc.insertion_loss, 1)) # IL neta = R - C por banda
+enc.plot() # panel R, corrección C e IL
+```
+
+`enclosure_insertion_loss` devuelve un `EnclosureResult` con la
+`panel_transmission_loss` del panel, la `correction` interior, la `insertion_loss`
+neta, la `room_constant` interior y `.plot()`.
+
+## Contraste con el solucionador FDTD
+
+La cámara de expansión de cuatro polos se contrasta con el
+[solucionador de ondas FDTD 2D](/phonometry/es/guides/fdtd-simulation/)
+independiente: un conducto de onda plana que se ensancha en una cámara y vuelve a
+estrecharse transmite mucho menos en el máximo de TL de cuatro polos ($kL = \pi/2$)
+que en la caída transparente ($kL = \pi$), y la razón de amplitudes medida
+reproduce la pérdida de transmisión máxima en forma cerrada con una fracción de
+decibelio.
+
+## Véase también
+
+- [Electroacústica](/phonometry/es/guides/electroacoustics/): el pistón radiante
+ (impedancia de radiación y directividad), el modelo de radiador complementario.
+- [Potencia sonora](/phonometry/es/guides/sound-power/): la `Lw` de fuente que
+ alimenta un conducto o un cerramiento.
+- [Fuentes imagen y campo estacionario](/phonometry/es/guides/room-image-sources/):
+ la `room_constant` reutilizada por la corrección interior del cerramiento.
+- [Simulación de ondas FDTD 2D](/phonometry/es/guides/fdtd-simulation/): el
+ solucionador independiente que contrasta la cámara de expansión.
+- [Informe de conformidad](https://github.com/jmrplens/phonometry/blob/main/docs/CONFORMANCE.md):
+ las formas cerradas y los anclajes resueltos con los que se validan estas
+ implementaciones.
+- Referencia de la API:
+ [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/),
+ [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) y
+ [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/).
diff --git a/site/src/content/docs/guides/electroacoustics.mdx b/site/src/content/docs/guides/electroacoustics.mdx
index 8424023e..2f808ea4 100644
--- a/site/src/content/docs/guides/electroacoustics.mdx
+++ b/site/src/content/docs/guides/electroacoustics.mdx
@@ -346,6 +346,48 @@ hide in that sentence:
"1 m" figure is then a referred quantity, not what a microphone placed at
1 m would read.
+## 6. Radiating piston: radiation impedance and directivity
+
+The rigid circular **piston in an infinite baffle** is the canonical radiator
+behind a loudspeaker cone, the open end of a duct and the radiation efficiency
+of any finite vibrating surface (Beranek & Mellow §4.19, §13.7). Its mechanical
+radiation impedance is $Z_r = \rho c S\,(R_1 + jX_1)$ with $S = \pi a^2$ and the
+dimensionless resistance and reactance functions (Eqs. (13.117), (13.118))
+
+$$
+R_1(x) = 1 - \frac{2 J_1(x)}{x}, \qquad X_1(x) = \frac{2 H_1(x)}{x},
+\qquad x = 2ka,
+$$
+
+where $J_1$ is the Bessel function and $H_1$ the Struve function of order one.
+At low frequency $R_1 \to (ka)^2/2$ and the reactance is mass-like with the
+**radiation mass** $M_r = 8\rho a^3/3$ (Eq. (4.151)); at high frequency
+$R_1 \to 1$, $X_1 \to 0$ and the piston radiates as if into an infinite tube.
+The far field follows the directivity $D(\theta) = 2 J_1(ka\sin\theta)/(ka\sin\theta)$,
+whose first null is at $ka\sin\theta = 3.8317$.
+
+```python
+import numpy as np
+from phonometry import radiating_piston
+
+res = radiating_piston(radius=0.1, frequencies=np.geomspace(20, 20000, 200),
+ angles=np.linspace(0.0, np.pi / 2, 91))
+print(round(res.radiation_mass, 4)) # 8 rho a^3 / 3, kg
+print(round(float(res.directivity_index[0]), 2)) # 3.01 dB half-space limit
+res.plot() # R1 and X1 vs ka
+```
+
+`radiating_piston` returns a `RadiatingPistonResult` with the normalized
+`resistance`/`reactance`, the mechanical `radiation_resistance`/`radiation_reactance`,
+the `radiation_mass`, the `directivity_index`, the far-field `directivity`
+pattern (when `angles` are given) and `.plot()`. The building blocks
+`piston_resistance`, `piston_reactance` and `piston_directivity` are also
+callable directly. The piston is the companion radiator model of the
+[industrial noise-control silencers](/phonometry/guides/noise-control/).
+
## See also
-- API reference: [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/) and [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/).
+- [Industrial noise control](/phonometry/guides/noise-control/): reactive
+ silencers, HVAC duct methods and machine enclosures, sharing the radiating
+ piston as the canonical radiator.
+- API reference: [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/), [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/) and [`electroacoustics.piston`](/phonometry/reference/api/electroacoustics/piston/).
diff --git a/site/src/content/docs/guides/noise-control.mdx b/site/src/content/docs/guides/noise-control.mdx
new file mode 100644
index 00000000..bf1701f3
--- /dev/null
+++ b/site/src/content/docs/guides/noise-control.mdx
@@ -0,0 +1,272 @@
+---
+title: "Industrial Noise Control: Silencers, HVAC and Enclosures"
+description: "Reactive silencers by the four-pole transmission-matrix method (expansion chambers, Helmholtz and quarter-wave resonators, extended tubes), HVAC duct attenuation and flow noise, and machine-enclosure insertion loss, from Bies, Hansen & Howard."
+references:
+ - type: book
+ authors: ["Bies, D. A.", "Hansen, C. H.", "Howard, C. Q."]
+ year: 2017
+ title: "Engineering noise control"
+ edition: "5th ed."
+ publisher: "CRC Press"
+ doi: "10.1201/9781351228152"
+ note: "The muffler four-pole method and expansion-chamber TL (§8.8-8.9), the HVAC duct methods (§8.11-8.17) and the machine-enclosure noise reduction (§7.4) of this page."
+ - type: book
+ authors: ["Munjal, M. L."]
+ year: 2014
+ title: "Acoustics of ducts and mufflers"
+ edition: "2nd ed."
+ publisher: "Wiley"
+ doi: "10.1002/9781118443767"
+ note: "The transfer-matrix formulation and the element matrices behind §1."
+ - type: book
+ authors: ["Vér, I. L.", "Beranek, L. L."]
+ year: 2006
+ title: "Noise and vibration control engineering: Principles and applications"
+ edition: "2nd ed."
+ publisher: "Wiley"
+ doi: "10.1002/9780470172568"
+ note: "The companion treatment of mufflers, ducts and enclosures cross-checked in this page."
+---
+
+import ThemeImage from '../../../components/ThemeImage.astro';
+
+Three passive measures dominate applied noise control: **silencers** in a duct,
+the passive attenuations and regenerated noise of an **HVAC** run, and a
+**machine enclosure**. `phonometry.noise_control` covers all three with the
+engineering theory of Bies, Hansen & Howard. The silencers are built on the
+one-dimensional four-pole (transfer-matrix) method, which chains reusable
+acoustic elements; the HVAC methods are the ASHRAE tables and closed forms for
+bends, end reflection, plenums and flow noise; and the enclosure model combines
+a user-supplied panel transmission loss with the reverberant build-up inside
+the enclosure. The radiating [piston](/phonometry/guides/electroacoustics/) of
+the electroacoustics domain is the companion radiator model.
+
+## 1. Reactive silencers (four-pole method)
+
+A reactive silencer attenuates by *reflecting* sound with impedance
+discontinuities. Each acoustic element is a 2x2 **transfer matrix** relating the
+pressure `p` and volume velocity `Su` at its two ends, and a compound silencer
+is the ordered matrix product of its elements (Bies §8.9). A straight duct of
+length `L` and area `S` is
+
+$$
+\begin{bmatrix} \cos kL & j\,\tfrac{\rho c}{S}\sin kL \\[2pt]
+j\,\tfrac{S}{\rho c}\sin kL & \cos kL \end{bmatrix},
+\qquad k = \omega/c,
+$$
+
+and a side branch of impedance $Z_b$ is the shunt
+$\left[\begin{smallmatrix} 1 & 0 \\ 1/Z_b & 1 \end{smallmatrix}\right]$. The
+**transmission loss** follows from the compound matrix $T$ (Bies Eq. (8.141);
+for equal inlet/outlet areas $S$, Eq. (8.148))
+
+$$
+\mathrm{TL} = 20\log_{10}\!\left(\tfrac{1}{2}\left|\,T_{11}
++ \tfrac{T_{12}}{Z_c} + Z_c\,T_{21} + T_{22}\right|\right),
+\qquad Z_c = \frac{\rho c}{S},
+$$
+
+and the **insertion loss** for a source impedance $Z_s$ and radiation impedance
+$Z_r$ is the extra attenuation over a direct connection.
+
+### Expansion chamber
+
+The simplest silencer, a chamber of area $S_\text{exp}$ and length $L$ between
+pipes of area $S_\text{duct}$, has the closed-form transmission loss (Bies
+Eq. (8.111)) with area ratio $m = S_\text{exp}/S_\text{duct}$
+
+$$
+\mathrm{TL} = 10\log_{10}\!\left[1 + \tfrac{1}{4}\left(m - \tfrac{1}{m}\right)^2
+\sin^2 kL\right],
+$$
+
+peaking at $10\log_{10}[1 + \tfrac14(m-1/m)^2]$ when $kL = \pi/2, 3\pi/2, \dots$
+and dropping to $0$ at $kL = n\pi$ where the chamber is a half-wavelength long
+and transparent. The four-pole product reproduces this exactly.
+
+```python
+import numpy as np
+from phonometry import expansion_chamber
+
+freqs = np.linspace(20.0, 2000.0, 2000)
+res = expansion_chamber(freqs, length=0.3, chamber_area=0.04, pipe_area=0.01)
+print(round(res.transmission_loss.max(), 2)) # 6.55 dB peak (m = 4)
+# The troughs at f = n c / 2L are exactly 0 dB (no dissipation).
+print(round(float(res.transmission_loss[np.argmin(res.transmission_loss)]), 6))
+```
+
+The figure below shows the transmission loss for four area ratios: a larger
+mismatch $m$ lifts every peak, but the troughs stay at 0 dB and the peaks stay
+at the same frequencies, set only by the chamber length.
+
+
+Show the code for this figure
+
+```python
+import matplotlib.pyplot as plt
+import numpy as np
+from phonometry import expansion_chamber
+
+freqs = np.linspace(20.0, 2000.0, 2000)
+
+# One line for one chamber: TL vs frequency (with the insertion loss too if
+# source/radiation impedances are given).
+expansion_chamber(freqs, 0.3, 0.04, 0.01,
+ source_impedance=4e4, radiation_impedance=5e3).plot()
+plt.show()
+
+# By hand: the family of area ratios of the concept figure.
+fig, ax = plt.subplots()
+for m in (2.0, 4.0, 8.0, 16.0):
+ res = expansion_chamber(freqs, 0.3, m * 0.01, 0.01)
+ ax.plot(freqs, res.transmission_loss, label=f"m = {int(m)}")
+ax.set_xlabel("Frequency [Hz]"); ax.set_ylabel("Transmission loss [dB]")
+ax.legend()
+plt.show()
+```
+
+
+
+### Side-branch and extended-tube resonators
+
+A **Helmholtz resonator** (neck area $S_n$, effective length $l_e$, cavity
+volume $V$) and a closed **quarter-wave tube** (length $l_e$) each short the
+duct at their tuning frequency, giving a sharp transmission-loss spike there:
+$f_0 = \tfrac{c}{2\pi}\sqrt{S_n/(l_e V)}$ (Bies Eq. (8.46)) and $f = c/4l_e$
+(Eq. (8.44)). An **extended-tube chamber** buries quarter-wave side branches in
+an expansion chamber to fill the plain chamber's troughs.
+
+```python
+import numpy as np
+from phonometry import (
+ helmholtz_resonator, quarter_wave_resonator, extended_tube_chamber,
+)
+
+f = np.linspace(20.0, 600.0, 4000)
+
+hr = helmholtz_resonator(f, duct_area=0.01, neck_area=1e-4,
+ neck_length=0.02, cavity_volume=1e-3)
+print(round(float(hr.resonances[0]), 1)) # tuning frequency, Hz
+
+qw = quarter_wave_resonator(f, duct_area=0.01, length=1.516, branch_area=2e-3,
+ speed_of_sound=343.24)
+print(round(float(qw.resonances[0]), 1)) # 56.6 Hz (Bies Example 8.1)
+
+# An inlet extension of L/4 fills the first expansion-chamber trough.
+et = extended_tube_chamber(f, length=0.4, chamber_area=0.04, pipe_area=0.01,
+ inlet_extension=0.1)
+```
+
+Each returns a `ReactiveSilencerResult` with `transmission_loss`,
+`insertion_loss` (when impedances are given), the compound `transfer_matrix`
+and `.plot()`. Advanced layouts chain elements directly with `duct_matrix`,
+`shunt_matrix`, `cascade`, `transmission_loss` and `insertion_loss`.
+
+## 2. HVAC duct attenuation and flow noise
+
+A ventilation run attenuates fan noise at bends, at the open duct end and in
+plenums, and regenerates noise wherever the airflow is disturbed. `phonometry.noise_control.hvac`
+gathers the Bies Chapter 8 methods.
+
+```python
+from phonometry.noise_control import hvac
+
+bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0]
+
+# Low-frequency reflection back up an open duct end (ASHRAE Table 8.14).
+er = hvac.end_reflection_loss(bands, diameter=0.30, termination="flush")
+
+# Insertion loss of a lined square 90-degree elbow (ASHRAE Table 8.11).
+el = hvac.elbow_insertion_loss(bands, width=0.3, bend_type="square", lined=True)
+
+# Plenum-chamber TL by Wells' method (closed form).
+tl = hvac.plenum_attenuation(exit_area=0.1, line_of_sight=1.0,
+ wall_area=20.0, mean_absorption=0.2)
+print(round(tl, 1)) # dB
+
+# Flow-generated (self) noise of a straight duct (VDI 2081).
+fn = hvac.flow_noise_straight_duct(bands, flow_velocity=10.0, area=0.04)
+```
+
+The end-reflection and elbow methods interpolate the ASHRAE tables (they pass
+exactly through the tabulated nodes); the plenum (Wells) and flow-noise (VDI
+2081) methods are closed forms. `end_reflection_loss`, `elbow_insertion_loss`,
+`flow_noise_straight_duct` and `flow_noise_bend` return an `HvacSpectrumResult`
+(attenuation or regenerated sound power level) with `.plot()`; `plenum_attenuation`
+returns the transmission loss directly. Rectangular ducts use the equivalent
+diameter $D = \sqrt{4S/\pi}$.
+
+:::note
+Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14
+look-up; there is no closed form in that edition. This module reproduces and
+interpolates the table.
+:::
+
+## 3. Machine enclosures
+
+A sealed enclosure reduces the radiated noise by its panel transmission loss
+$R$, minus a penalty $C$ for the reverberant build-up inside the small, hard
+cavity (Bies Eqs. (7.103), (7.111)):
+
+$$
+\mathrm{IL} = R - C,\qquad C = 10\log_{10}\!\left(0.3 + \frac{S_E}{R_i}\right),
+$$
+
+with the external area $S_E$ and the interior **room constant**
+$R_i = S_i \alpha_i/(1-\alpha_i)$ (the same
+[`room_constant`](/phonometry/guides/room-image-sources/) as the steady-state
+room field). A hard interior wastes much of the panel $R$; lining it drives $C$
+toward its floor $10\log_{10}0.3 = -5.2$ dB.
+
+**The panel transmission loss $R$ is supplied by you** as a per-band array
+(measured, or predicted by a panel model), or as a callable of frequency. This
+module never predicts $R$ itself; it combines a given $R$ with the interior
+absorption.
+
+```python
+import numpy as np
+from phonometry import enclosure_insertion_loss
+
+bands = np.array([125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0])
+panel_R = np.array([18.0, 24.0, 30.0, 36.0, 42.0, 46.0]) # measured, dB
+
+enc = enclosure_insertion_loss(
+ panel_R, external_area=6.0, internal_area=5.0,
+ internal_absorption=0.3, frequencies=bands,
+)
+print(np.round(enc.insertion_loss, 1)) # net IL = R - C per band
+enc.plot() # panel R, correction C and IL
+```
+
+`enclosure_insertion_loss` returns an `EnclosureResult` with the panel
+`panel_transmission_loss`, the interior `correction`, the net `insertion_loss`,
+the interior `room_constant`, and `.plot()`.
+
+## Cross-check against the FDTD solver
+
+The four-pole expansion chamber is cross-checked against the independent 2D
+[FDTD wave solver](/phonometry/guides/fdtd-simulation/): a plane-wave duct that
+widens into a chamber and narrows back transmits far less at the four-pole TL
+peak ($kL = \pi/2$) than at the transparent trough ($kL = \pi$), and the
+measured amplitude ratio reproduces the closed-form peak transmission loss to a
+fraction of a decibel.
+
+## See also
+
+- [Electroacoustics](/phonometry/guides/electroacoustics/): the radiating
+ piston (radiation impedance and directivity), the companion radiator model.
+- [Sound Power](/phonometry/guides/sound-power/): the source `Lw` that feeds a
+ duct or an enclosure.
+- [Room image sources and steady field](/phonometry/guides/room-image-sources/):
+ the `room_constant` reused by the enclosure interior correction.
+- [2D FDTD wave simulation](/phonometry/guides/fdtd-simulation/): the
+ independent solver that cross-checks the expansion chamber.
+- [Conformance report](https://github.com/jmrplens/phonometry/blob/main/docs/CONFORMANCE.md):
+ the closed forms and worked anchors these implementations are validated
+ against.
+- API reference:
+ [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/),
+ [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) and
+ [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/).
diff --git a/site/src/content/docs/reference/api/electroacoustics/piston.md b/site/src/content/docs/reference/api/electroacoustics/piston.md
new file mode 100644
index 00000000..a8a01296
--- /dev/null
+++ b/site/src/content/docs/reference/api/electroacoustics/piston.md
@@ -0,0 +1,209 @@
+---
+title: "electroacoustics.piston"
+description: "Public API of phonometry.electroacoustics.piston (auto-generated)."
+sidebar:
+ label: "piston"
+---
+
+> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand.
+
+Radiation of a rigid circular piston set in an infinite baffle.
+
+The baffled circular piston is the canonical acoustic radiator: a flat rigid
+disc of radius `a` vibrating with a uniform normal velocity in an otherwise
+rigid infinite plane. It is the model behind a loudspeaker cone in a large
+cabinet, the open end of a duct, and the reference source for the radiation
+efficiency of any finite vibrating surface, so its two results -- the
+**radiation impedance** the air presents to the piston and the **directivity**
+of the far field -- are the base of the electroacoustics domain (Beranek &
+Mellow, *Acoustics: Sound Fields, Transducers and Vibration* 2nd ed., §4.4;
+Bies, Hansen & Howard, *Engineering Noise Control* 5th ed.).
+
+**Radiation impedance.** The reaction force of the air on the piston is
+`F = Z_r u` with the mechanical radiation impedance
+
+ Z_r = rho c S ( R1(2ka) + j X1(2ka) ), S = pi a^2,
+
+where `k = omega / c` is the wavenumber, `rho c` the characteristic
+impedance of air and `S` the piston area. The dimensionless **piston
+resistance** and **reactance** functions are (Beranek & Mellow Eq. (4.30))
+
+ R1(x) = 1 - 2 J1(x) / x, X1(x) = 2 H1(x) / x,
+
+with `J1` the Bessel function of the first kind and `H1` the Struve
+function, both of order one, evaluated at `x = 2ka`.
+
+* **Low frequency** (`ka << 1`): `R1 -> (ka)^2 / 2` so the radiated power
+ rises as `f^2`, and `X1 -> (8 / 3 pi) ka`. The reactance is mass-like,
+ `X_r = rho c S X1 = omega M_r` with the **radiation (accreted) mass**
+
+ M_r = 8 rho a^3 / 3
+
+ (Beranek & Mellow Eq. (4.32)): the piston drags an extra `8 rho a^3 / 3` of
+ air, equivalent to a layer `8a / 3 pi` thick over its face.
+* **High frequency** (`ka >> 1`): `R1 -> 1` and `X1 -> 0`, so
+ `Z_r -> rho c S` -- the piston radiates as if into an infinite tube and the
+ air loads it purely resistively.
+
+**Directivity.** The far-field pressure of the baffled piston varies with the
+polar angle `theta` from the axis as (Beranek & Mellow Eq. (4.42))
+
+ D(theta) = 2 J1(ka sin theta) / (ka sin theta), D(0) = 1.
+
+The main lobe narrows as `ka` grows; its first null is at
+`ka sin theta = 3.8317` (the first zero of `J1`), which exists only once
+`ka > 3.8317`. The **directivity factor** `Q` (on-axis intensity over the
+intensity of a point source of equal power radiating into the full sphere) and
+the **directivity index** `DI = 10 log10 Q` follow from integrating
+`|D|^2` over the radiating hemisphere,
+
+ Q = 2 / integral_0^(pi/2) |D(theta)|^2 sin theta d theta,
+
+which tends to `Q = 2` (`DI = 3.01 dB`, the half-space baffle gain) at low
+`ka` and to `Q ~ (ka)^2` (`DI ~ 20 log10 ka`) at high `ka`.
+
+## piston_directivity
+
+```python
+piston_directivity(ka: ArrayLike, theta: ArrayLike) -> np.ndarray | float
+```
+
+Far-field directivity `D = 2 J1(ka sin theta) / (ka sin theta)`.
+
+The pressure amplitude of a baffled circular piston relative to its on-axis
+value (Beranek & Mellow Eq. (4.42)), normalized so `D(0) = 1`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `ka` | Wavenumber-radius product `ka` (scalar or array). |
+| `theta` | Polar angle from the axis, rad (scalar or array). Broadcast against `ka`. |
+
+**Returns:** `D` (float for scalar inputs, else an array).
+
+## piston_reactance
+
+```python
+piston_reactance(x: ArrayLike) -> np.ndarray | float
+```
+
+Piston reactance function `X1(x) = 2 H1(x) / x` (`H1` Struve order 1).
+
+The imaginary part of the normalized radiation impedance of a baffled
+circular piston (Beranek & Mellow Eq. (4.30)). It rises as
+`(8 / 3 pi) ka` (mass-like) at low `x = 2ka` and decays to 0 at high
+`x`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `x` | Argument `x = 2ka` (scalar or array), dimensionless. |
+
+**Returns:** `X1(x)` (float for scalar input, else an array).
+
+## piston_resistance
+
+```python
+piston_resistance(x: ArrayLike) -> np.ndarray | float
+```
+
+Piston resistance function `R1(x) = 1 - 2 J1(x) / x`.
+
+The real part of the normalized radiation impedance of a baffled circular
+piston, as a function of `x = 2ka` (Beranek & Mellow Eq. (4.30)). It
+rises as `x^2 / 8 = (ka)^2 / 2` at low `x` and tends to 1 at high `x`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `x` | Argument `x = 2ka` (scalar or array), dimensionless. |
+
+**Returns:** `R1(x)` (float for scalar input, else an array).
+
+## radiating_piston
+
+```python
+radiating_piston(
+ radius: float,
+ frequencies: ArrayLike,
+ *,
+ speed_of_sound: float = 343.0,
+ density: float = 1.206,
+ angles: ArrayLike | None = None,
+) -> RadiatingPistonResult
+```
+
+Radiation impedance and directivity of a rigid baffled circular piston.
+
+Evaluates the piston resistance `R1(2ka)` and reactance `X1(2ka)`, the
+mechanical radiation impedance `rho c S (R1 + j X1)`, the low-frequency
+radiation mass `8 rho a^3 / 3` and the directivity index over the given
+frequencies (Beranek & Mellow §4.4). Pass `angles` to also sample the
+far-field directivity pattern `D(theta)`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `radius` | Piston radius `a`, m. |
+| `frequencies` | Frequencies `f`, Hz (scalar or 1-D array), all > 0. |
+| `speed_of_sound` | Speed of sound `c`, m/s (default 343). |
+| `density` | Air density `rho`, kg/m3 (default 1.206). |
+| `angles` | Optional polar angles `theta` from the axis, rad, at which to sample the directivity pattern. |
+
+**Returns:** A [`RadiatingPistonResult`](/phonometry/reference/api/electroacoustics/piston/#radiatingpistonresult).
+
+## RadiatingPistonResult
+
+```python
+RadiatingPistonResult(
+ frequencies: np.ndarray,
+ ka: np.ndarray,
+ resistance: np.ndarray,
+ reactance: np.ndarray,
+ radiation_resistance: np.ndarray,
+ radiation_reactance: np.ndarray,
+ radiation_mass: float,
+ directivity_index: np.ndarray,
+ angles: np.ndarray | None,
+ directivity: np.ndarray | None,
+ radius: float,
+ speed_of_sound: float,
+ density: float,
+)
+```
+
+Radiation impedance and directivity of a baffled circular piston.
+
+**Attributes**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz. |
+| `ka` | Wavenumber-radius product `ka` at each frequency. |
+| `resistance` | Normalized piston resistance `R1(2ka)` (real part of `Z_r / (rho c S)`). |
+| `reactance` | Normalized piston reactance `X1(2ka)` (imaginary part of `Z_r / (rho c S)`). |
+| `radiation_resistance` | Mechanical radiation resistance `rho c S R1`, N s/m. |
+| `radiation_reactance` | Mechanical radiation reactance `rho c S X1`, N s/m. |
+| `radiation_mass` | Low-frequency accreted air mass `M_r = 8 rho a^3/3`, kg (a single value; the mass limit of `radiation_reactance / omega`). |
+| `directivity_index` | Directivity index `DI = 10 log10 Q`, dB. |
+| `angles` | Polar angles of `directivity`, rad, or `None` if not requested. |
+| `directivity` | Far-field directivity `D(theta)` as a `(n_freq, n_angle)` array, or `None` if `angles` was not given. |
+| `radius` | Piston radius `a`, m. |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+| `density` | Air density `rho`, kg/m3. |
+
+### RadiatingPistonResult.plot()
+
+```python
+RadiatingPistonResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes
+```
+
+Plot the normalized piston resistance and reactance against `ka`.
+
+Reproduces the classic Beranek & Mellow figure: `R1` rising to 1 and
+`X1` peaking then decaying, over the `ka` range of the result.
+Requires matplotlib (`pip install phonometry[plot]`).
diff --git a/site/src/content/docs/reference/api/index.md b/site/src/content/docs/reference/api/index.md
index 93a41e24..99040fe6 100644
--- a/site/src/content/docs/reference/api/index.md
+++ b/site/src/content/docs/reference/api/index.md
@@ -183,6 +183,15 @@ La referencia de la API se genera a partir de los docstrings del código (en ing
| [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/) | Distortion metrics for electroacoustic equipment (IEC 60268-3 / AES17). |
| [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/) | Frequency-response and coherence estimators (Bendat & Piersol). |
| [`electroacoustics.swept_sine`](/phonometry/reference/api/electroacoustics/swept-sine/) | Harmonic-distortion separation with exponential sweeps (Farina / Novak). |
+| [`electroacoustics.piston`](/phonometry/reference/api/electroacoustics/piston/) | Radiation of a rigid circular piston set in an infinite baffle. |
+
+## Industrial noise control
+
+| Module | Summary |
+| :--- | :--- |
+| [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/) | Reactive silencers by the four-pole (transmission-matrix) method. |
+| [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) | HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise. |
+| [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/) | Insertion loss of a close or free-standing machine enclosure. |
## Program loudness
diff --git a/site/src/content/docs/reference/api/noise_control/enclosures.md b/site/src/content/docs/reference/api/noise_control/enclosures.md
new file mode 100644
index 00000000..986ce7e9
--- /dev/null
+++ b/site/src/content/docs/reference/api/noise_control/enclosures.md
@@ -0,0 +1,106 @@
+---
+title: "noise_control.enclosures"
+description: "Public API of phonometry.noise_control.enclosures (auto-generated)."
+sidebar:
+ label: "enclosures"
+---
+
+> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand.
+
+Insertion loss of a close or free-standing machine enclosure.
+
+Wrapping a machine in a sealed enclosure reduces the radiated noise by the
+transmission loss of its panels, *minus* a penalty for the reverberant build-up
+inside the small, hard cavity. Bies, Hansen & Howard, *Engineering Noise
+Control* 5th ed., §7.4.2 (Eqs. (7.103), (7.111)) write the net reduction as
+
+ IL = R - C, C = 10 log10[ 0.3 + S_E (1 - alpha_i) / (S_i alpha_i) ],
+
+where `R` is the field-incidence transmission loss of the enclosure panels,
+`S_E` the external surface area, `S_i` the internal surface area (including
+the machine) and `alpha_i` the mean absorption of the enclosure interior. The
+reverberant term is exactly `S_E` over the interior **room constant**
+`R_i = S_i alpha_i / (1 - alpha_i)` ([`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant)), so
+
+ C = 10 log10( 0.3 + S_E / R_i ).
+
+A hard interior (`alpha_i` small) makes `C` large and wastes much of the
+panel `R`; lining the enclosure drives `C` toward its floor
+`10 log10 0.3 = -5.2 dB` (a fully absorbing interior, where `IL = R + 5.2`).
+Bies terms this net reduction the enclosure *noise reduction*; it is the
+insertion loss of the enclosure.
+
+**The panel transmission loss `R` is supplied by the caller** -- measured, or
+predicted by a panel model -- as a per-band array or a callable of frequency.
+This module never predicts `R` itself; it combines a given `R` with the
+interior absorption. The interior room constant reuses
+[`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant).
+
+## enclosure_insertion_loss
+
+```python
+enclosure_insertion_loss(
+ panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike],
+ external_area: float,
+ internal_area: float,
+ internal_absorption: ArrayLike,
+ *,
+ frequencies: ArrayLike | None = None,
+) -> EnclosureResult
+```
+
+Net insertion loss of a machine enclosure (Bies Eqs. (7.103), (7.111)).
+
+`IL = R - C` with `C = 10 log10(0.3 + S_E / R_i)` and the interior room
+constant `R_i = S_i alpha_i / (1 - alpha_i)`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `panel_transmission_loss` | Panel transmission loss `R` per band, dB. Either a per-band array (measured or predicted elsewhere) or a callable mapping a frequency array to per-band `R` (then `frequencies` is required). This function does not predict `R`. |
+| `external_area` | External enclosure surface area `S_E`, m2. |
+| `internal_area` | Internal surface area `S_i` (including the machine), m2. |
+| `internal_absorption` | Mean interior absorption `alpha_i` in `(0, 1)` (scalar or per-band). |
+| `frequencies` | Band centre frequencies, Hz; required when `panel_transmission_loss` is a callable, optional otherwise (used to label the result and the plot). |
+
+**Returns:** An [`EnclosureResult`](/phonometry/reference/api/noise_control/enclosures/#enclosureresult).
+
+## EnclosureResult
+
+```python
+EnclosureResult(
+ frequencies: np.ndarray | None,
+ panel_transmission_loss: np.ndarray,
+ correction: np.ndarray,
+ insertion_loss: np.ndarray,
+ external_area: float,
+ internal_area: float,
+ room_constant: np.ndarray,
+)
+```
+
+Insertion loss of a machine enclosure over frequency (Bies §7.4.2).
+
+**Attributes**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz, or `None` if the panel `R` was given as a bare per-band array with no frequency labels. |
+| `panel_transmission_loss` | The supplied panel transmission loss `R` per band, dB. |
+| `correction` | The interior-build-up correction `C` per band, dB. |
+| `insertion_loss` | The net enclosure insertion loss `IL = R - C`, dB. |
+| `external_area` | External enclosure surface area `S_E`, m2. |
+| `internal_area` | Internal surface area `S_i`, m2. |
+| `room_constant` | Interior room constant `R_i` per band, m2. |
+
+### EnclosureResult.plot()
+
+```python
+EnclosureResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes
+```
+
+Plot the panel `R`, correction `C` and net insertion loss.
+
+Requires matplotlib (`pip install phonometry[plot]`); returns the
+`Axes`.
diff --git a/site/src/content/docs/reference/api/noise_control/hvac.md b/site/src/content/docs/reference/api/noise_control/hvac.md
new file mode 100644
index 00000000..4b9123ee
--- /dev/null
+++ b/site/src/content/docs/reference/api/noise_control/hvac.md
@@ -0,0 +1,220 @@
+---
+title: "noise_control.hvac"
+description: "Public API of phonometry.noise_control.hvac (auto-generated)."
+sidebar:
+ label: "hvac"
+---
+
+> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand.
+
+HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise.
+
+A ventilation duct network attenuates fan noise through several mechanisms
+that add up along the path, and it *regenerates* noise wherever the airflow is
+disturbed. This module gathers the engineering methods of Bies, Hansen &
+Howard, *Engineering Noise Control* 5th ed., Chapter 8, for the passive
+attenuations -- **duct end reflection** (§8.13, Table 8.14), **bends/elbows**
+(§8.11, Table 8.11) and **plenum chambers** (§8.17, Wells' method) -- and for
+the **flow-generated (self) noise** of straight ducts and bends (§8.15).
+
+The end-reflection and elbow methods are empirical look-up tables (ASHRAE);
+they are interpolated over the duct size and, for the elbows, over the
+frequency-to-width ratio `W / lambda`. The plenum and flow-noise methods are
+closed forms evaluated directly.
+
+:::note
+Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14
+look-up (there is no closed form in this edition); this module reproduces
+that table and interpolates it. Rectangular ducts use the equivalent
+diameter `D = sqrt(4 S / pi)`.
+:::
+
+## elbow_insertion_loss
+
+```python
+elbow_insertion_loss(
+ frequencies: ArrayLike,
+ width: float,
+ *,
+ bend_type: str = 'square',
+ vanes: bool = False,
+ lined: bool = False,
+ speed_of_sound: float = 343.0,
+) -> HvacSpectrumResult
+```
+
+Duct bend/elbow insertion loss per bend (Bies Table 8.11, ASHRAE).
+
+Indexed by the frequency-to-width ratio `W / lambda` (`lambda = c / f`).
+Lined bends assume the lining extends at least three duct diameters up- and
+downstream. Round bends are treated as unlined with no vanes.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `width` | Duct width `W` in the plane of the bend, m. |
+| `bend_type` | `"square"` or `"round"`. |
+| `vanes` | Turning vanes fitted (square bends only). |
+| `lined` | Acoustically lined bend (square bends only). |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+
+**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the insertion loss, dB per bend.
+
+## end_reflection_loss
+
+```python
+end_reflection_loss(
+ frequencies: ArrayLike,
+ diameter: float,
+ *,
+ termination: str = 'flush',
+ speed_of_sound: float = 343.0,
+) -> HvacSpectrumResult
+```
+
+Duct end reflection loss (Bies Table 8.14, ASHRAE).
+
+The low-frequency reflection of sound back up a duct at its open
+termination into a room. Interpolated over `log` diameter and `log`
+frequency from Table 8.14; it passes exactly through the tabulated
+`(diameter, octave band)` nodes.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `diameter` | Duct internal diameter `D`, m (use `D = sqrt(4 S / pi)` for a rectangular duct of area `S`). |
+| `termination` | `"flush"` (duct flush with a wall/ceiling) or `"free"` (free space / suspended in the room). |
+| `speed_of_sound` | Speed of sound `c`, m/s (kept for signature symmetry; the table is indexed by frequency directly). |
+
+**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the reflection loss, dB.
+
+## flow_noise_bend
+
+```python
+flow_noise_bend(
+ frequencies: ArrayLike,
+ flow_velocity: float,
+ area: float,
+ height: float,
+ *,
+ density: float = 1.206,
+) -> HvacSpectrumResult
+```
+
+Flow-generated octave-band sound power of a mitred bend (Bies Eqs. (8.252), (8.254)).
+
+`L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103` with the
+stream power level `L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117`
+(Bies Eq. (8.252)) and the Strouhal number `N_s = f H / U` (`H` the duct
+height in the plane of the bend). The radiated sound power grows as the
+sixth power of the stream speed at low `N_s` (the inner-corner drag
+dipole) and the eighth power at high `N_s` (the outer-corner shear
+quadrupole); equivalently, the *efficiency* referenced to the stream power
+grows as `U^3` and `U^5` respectively.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Octave-band centre frequencies `f`, Hz (1-D array). |
+| `flow_velocity` | Mean flow speed `U`, m/s. |
+| `area` | Duct cross-sectional area `S`, m2. |
+| `height` | Duct height `H` in the plane of the bend, m. |
+| `density` | Air density `rho`, kg/m3. |
+
+**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the band sound power level, dB re 1e-12 W.
+
+## flow_noise_straight_duct
+
+```python
+flow_noise_straight_duct(
+ frequencies: ArrayLike,
+ flow_velocity: float,
+ area: float,
+) -> HvacSpectrumResult
+```
+
+Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)).
+
+`L_WB = 7 + 50 log10(U) + 10 log10(S) - 2 - 26 log10(1.14 + 0.02 f / U)`
+in dB re 1e-12 W (VDI 2081-1), for airflow speed `U` in a duct of area
+`S`.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Octave-band centre frequencies `f`, Hz (1-D array). |
+| `flow_velocity` | Mean flow speed `U`, m/s. |
+| `area` | Duct cross-sectional area `S`, m2. |
+
+**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the band sound power level, dB re 1e-12 W.
+
+## HvacSpectrumResult
+
+```python
+HvacSpectrumResult(
+ frequencies: np.ndarray,
+ values: np.ndarray,
+ quantity: str,
+ label: str,
+)
+```
+
+A per-frequency HVAC quantity (attenuation or regenerated power level).
+
+**Attributes**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz. |
+| `values` | The quantity per frequency (dB, or dB re 1e-12 W for a sound power level). |
+| `quantity` | What `values` holds (`"attenuation"` or `"sound_power_level"`). |
+| `label` | A short human label of the element. |
+
+### HvacSpectrumResult.plot()
+
+```python
+HvacSpectrumResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes
+```
+
+Plot the quantity against a continuous log-frequency axis.
+
+Requires matplotlib (`pip install phonometry[plot]`).
+
+## plenum_attenuation
+
+```python
+plenum_attenuation(
+ exit_area: float,
+ line_of_sight: float,
+ wall_area: float,
+ mean_absorption: ArrayLike,
+ *,
+ angle: float = 0.0,
+) -> np.ndarray | float
+```
+
+Plenum-chamber transmission loss by Wells' method (Bies Eq. (8.275)).
+
+`TL = -10 log10[ S_out ( cos(theta) / (pi r^2) + (1 - alpha) / (S_w alpha) ) ]`,
+where the reverberant term uses the plenum room constant
+`R = S_w alpha / (1 - alpha)` ([`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant)). The
+method holds above the inlet cut-on and when the plenum is large compared
+with the wavelength; it underpredicts the low-frequency loss by 5-10 dB.
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `exit_area` | Outlet-opening area `S_out`, m2. |
+| `line_of_sight` | Straight-line inlet-to-outlet distance `r`, m. |
+| `wall_area` | Total internal wall area `S_w`, m2. |
+| `mean_absorption` | Mean Sabine wall absorption `alpha` in `(0, 1)` (scalar or per-band). |
+| `angle` | Angle `theta` between the inlet axis and the line to the outlet, rad (default 0). |
+
+**Returns:** The transmission loss, dB (float for scalar absorption, else a per-band array).
diff --git a/site/src/content/docs/reference/api/noise_control/silencers.md b/site/src/content/docs/reference/api/noise_control/silencers.md
new file mode 100644
index 00000000..506c48d4
--- /dev/null
+++ b/site/src/content/docs/reference/api/noise_control/silencers.md
@@ -0,0 +1,244 @@
+---
+title: "noise_control.silencers"
+description: "Public API of phonometry.noise_control.silencers (auto-generated)."
+sidebar:
+ label: "silencers"
+---
+
+> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand.
+
+Reactive silencers by the four-pole (transmission-matrix) method.
+
+A reactive silencer controls noise by *reflecting* it back to the source with
+impedance discontinuities -- sudden area changes and side branches -- rather
+than by dissipating it in absorptive material. The one-dimensional plane-wave
+theory represents each acoustic element by a 2x2 **transfer (four-pole)
+matrix** relating the sound pressure `p` and volume velocity `S u` at its
+two ends, and a compound silencer is the ordered matrix product of its
+elements (Bies, Hansen & Howard, *Engineering Noise Control* 5th ed., §8.8-8.9;
+Munjal, *Acoustics of Ducts and Mufflers*).
+
+**Transfer matrix** (Bies Eq. (8.133)), state vector `[p, S u]` with the
+characteristic acoustic impedance `Z = rho c / S`. The plane-wave element for
+a straight duct of length `L` and area `S` is (Bies Eq. (8.143), no flow)
+
+ [[ cos(kL), j (rho c / S) sin(kL) ],
+ [ j (S / rho c) sin(kL), cos(kL) ]], k = omega / c,
+
+and a **side branch** of acoustic impedance `Z_b` is the shunt element
+(Bies Eq. (8.144))
+
+ [[ 1, 0 ],
+ [ 1 / Z_b, 1 ]].
+
+**Transmission loss** from the compound matrix `T` (Bies Eq. (8.141), no
+flow; reduces to Eq. (8.148) for equal inlet/outlet areas):
+
+ TL = 10 log10[ (Z1 / Zn) (1/4) | T11 + T12 / Zn + Z1 T21 + (Z1 / Zn) T22 |^2 ]
+
+with `Z1 = rho c / S_in` and `Zn = rho c / S_out`. `TL` is the intrinsic
+attenuation for an anechoic termination. The **insertion loss** for a source of
+internal impedance `Z_s` radiating into a termination impedance `Z_r` is
+the extra attenuation of inserting the silencer in place of a direct
+connection,
+
+ IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) / (Z_s + Z_r) |,
+
+which is `0` when the silencer reduces to a through connection (`T = I`)
+and equals the transmission loss for the anechoic reference
+`Z_s = rho c / S_in`, `Z_r = rho c / S_out`.
+
+**Simple expansion chamber.** A chamber of area `S_exp` and length `L`
+between pipes of area `S_duct` has the closed-form transmission loss (Bies
+Eq. (8.111)) with area ratio `m = S_exp / S_duct`
+
+ TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ],
+
+peaking at `10 log10[1 + (1/4)(m - 1/m)^2]` when `kL = pi/2, 3pi/2, ...` and
+dropping to `0` at `kL = n pi` (no dissipation). The four-pole product
+reproduces this exactly, and the machinery extends to side-branch (Helmholtz,
+quarter-wave) and extended-tube resonators that the closed form cannot cover.
+
+## expansion_chamber
+
+```python
+expansion_chamber(
+ frequencies: ArrayLike,
+ length: float,
+ chamber_area: float,
+ pipe_area: float,
+ *,
+ speed_of_sound: float = 343.0,
+ density: float = 1.206,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult
+```
+
+Simple expansion-chamber silencer (Bies Eq. (8.111) / four-pole).
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `length` | Chamber length `L`, m. |
+| `chamber_area` | Chamber cross-sectional area `S_exp`, m2. |
+| `pipe_area` | Inlet/outlet pipe area `S_duct`, m2. |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+| `density` | Air density `rho`, kg/m3. |
+| `source_impedance` | Optional source impedance `Z_s` for the insertion loss, Pa s/m3. |
+| `radiation_impedance` | Optional radiation impedance `Z_r` for the insertion loss, Pa s/m3. |
+
+**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult) (its `transmission_loss` equals the closed form `10 log10[1 + (1/4)(m - 1/m)^2 sin^2(kL)]`).
+
+## extended_tube_chamber
+
+```python
+extended_tube_chamber(
+ frequencies: ArrayLike,
+ length: float,
+ chamber_area: float,
+ pipe_area: float,
+ *,
+ inlet_extension: float = 0.0,
+ outlet_extension: float = 0.0,
+ speed_of_sound: float = 343.0,
+ density: float = 1.206,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult
+```
+
+Extended-inlet/outlet expansion chamber (Bies §8.9.7).
+
+The inlet and outlet pipes extend a distance into the chamber, forming
+annular quarter-wave side branches (of area `S_exp - S_duct` and lengths
+equal to the extensions, Bies Eq. (8.156)) at the two junctions. Tuning the
+extensions (classically `L/4` and `L/2`) places quarter-wave peaks that
+fill the `kL = n pi` troughs of the plain expansion chamber. With both
+extensions `0` the result reduces exactly to [`expansion_chamber`](/phonometry/reference/api/noise_control/silencers/#expansion_chamber).
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `length` | Chamber length `L`, m. |
+| `chamber_area` | Chamber cross-sectional area `S_exp`, m2. |
+| `pipe_area` | Inlet/outlet pipe area `S_duct`, m2. |
+| `inlet_extension` | Inlet pipe extension into the chamber `L_a`, m. |
+| `outlet_extension` | Outlet pipe extension into the chamber `L_b`, m. |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+| `density` | Air density `rho`, kg/m3. |
+| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. |
+| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. |
+
+**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult).
+
+## helmholtz_resonator
+
+```python
+helmholtz_resonator(
+ frequencies: ArrayLike,
+ duct_area: float,
+ neck_area: float,
+ neck_length: float,
+ cavity_volume: float,
+ *,
+ resistance: float = 0.0,
+ speed_of_sound: float = 343.0,
+ density: float = 1.206,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult
+```
+
+Side-branch Helmholtz resonator on a duct (Bies Eqs. (8.144), (8.152)).
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `duct_area` | Main-duct cross-sectional area `S_d`, m2. |
+| `neck_area` | Resonator neck area `S_neck`, m2. |
+| `neck_length` | Effective neck length `l_e`, m. |
+| `cavity_volume` | Cavity volume `V`, m3. |
+| `resistance` | Neck acoustic resistance `R`, Pa s/m3 (default 0). |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+| `density` | Air density `rho`, kg/m3. |
+| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. |
+| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. |
+
+**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult); `resonances` holds `f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))`.
+
+## quarter_wave_resonator
+
+```python
+quarter_wave_resonator(
+ frequencies: ArrayLike,
+ duct_area: float,
+ length: float,
+ branch_area: float,
+ *,
+ speed_of_sound: float = 343.0,
+ density: float = 1.206,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult
+```
+
+Closed quarter-wave side-branch tube on a duct (Bies Eqs. (8.144), (8.146)).
+
+**Parameters**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz (1-D array). |
+| `duct_area` | Main-duct cross-sectional area `S_d`, m2. |
+| `length` | Effective branch length `l_e`, m. |
+| `branch_area` | Branch tube area `S`, m2. |
+| `speed_of_sound` | Speed of sound `c`, m/s. |
+| `density` | Air density `rho`, kg/m3. |
+| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. |
+| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. |
+
+**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult); `resonances` holds the odd multiples of `f = c / (4 l_e)` within the frequency range.
+
+## ReactiveSilencerResult
+
+```python
+ReactiveSilencerResult(
+ frequencies: np.ndarray,
+ transmission_loss: np.ndarray,
+ insertion_loss: np.ndarray | None,
+ transfer_matrix: np.ndarray,
+ kind: str,
+ resonances: np.ndarray | None = None,
+)
+```
+
+Transmission and insertion loss of a reactive silencer over frequency.
+
+**Attributes**
+
+| Name | Description |
+| :--- | :--- |
+| `frequencies` | Frequencies `f`, Hz. |
+| `transmission_loss` | Transmission loss per frequency, dB. |
+| `insertion_loss` | Insertion loss per frequency, dB, or `None` when no source/radiation impedance was supplied. |
+| `transfer_matrix` | The compound `(n_freq, 2, 2)` four-pole matrix. |
+| `kind` | A short label of the device (e.g. `"expansion chamber"`). |
+| `resonances` | Notable resonance frequencies, Hz (e.g. the resonator tuning frequency), or `None`. |
+
+### ReactiveSilencerResult.plot()
+
+```python
+ReactiveSilencerResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes
+```
+
+Plot the transmission (and insertion) loss against frequency.
+
+Requires matplotlib (`pip install phonometry[plot]`); returns the
+`Axes`.
diff --git a/site/src/generated/api-sidebar.mjs b/site/src/generated/api-sidebar.mjs
index bba59195..c1025a42 100644
--- a/site/src/generated/api-sidebar.mjs
+++ b/site/src/generated/api-sidebar.mjs
@@ -179,6 +179,16 @@ export const apiSidebar = {
'reference/api/electroacoustics/distortion',
'reference/api/electroacoustics/frequency-response',
'reference/api/electroacoustics/swept-sine',
+ 'reference/api/electroacoustics/piston',
+ ],
+ },
+ {
+ label: 'Industrial noise control',
+ translations: { es: 'Control de ruido industrial' },
+ items: [
+ 'reference/api/noise_control/silencers',
+ 'reference/api/noise_control/hvac',
+ 'reference/api/noise_control/enclosures',
],
},
{
diff --git a/src/phonometry/__init__.py b/src/phonometry/__init__.py
index 0b30156b..bf8fa342 100644
--- a/src/phonometry/__init__.py
+++ b/src/phonometry/__init__.py
@@ -309,6 +309,13 @@
coherence,
transfer_function,
)
+from .electroacoustics.piston import (
+ RadiatingPistonResult,
+ piston_directivity,
+ piston_reactance,
+ piston_resistance,
+ radiating_piston,
+)
from .electroacoustics.swept_sine import (
SweptSineDistortionResult,
swept_sine_distortion,
@@ -918,6 +925,22 @@
SignalSource,
fdtd_simulation,
)
+from .noise_control.silencers import (
+ ReactiveSilencerResult,
+ expansion_chamber,
+ extended_tube_chamber,
+ helmholtz_resonator,
+ quarter_wave_resonator,
+)
+from .noise_control.hvac import (
+ HvacSpectrumResult,
+ elbow_insertion_loss,
+ end_reflection_loss,
+ flow_noise_bend,
+ flow_noise_straight_duct,
+ plenum_attenuation,
+)
+from .noise_control.enclosures import EnclosureResult, enclosure_insertion_loss
from ._version import __version__
# Public methods
@@ -993,6 +1016,11 @@
"swept_sine_distortion",
"synchronized_sweep_signal",
"SweptSineDistortionResult",
+ "radiating_piston",
+ "piston_resistance",
+ "piston_reactance",
+ "piston_directivity",
+ "RadiatingPistonResult",
"program_loudness",
"integrated_loudness",
"loudness_range",
@@ -1694,6 +1722,23 @@
"circular_aperture_transmission_coefficient",
"composite_transmission_loss",
"transmission_loss_from_coefficient",
+ # Industrial noise control (silencers, HVAC, enclosures). The low-level
+ # four-pole primitives (duct_matrix, shunt_matrix, cascade, transmission_loss,
+ # insertion_loss, ...) stay under phonometry.noise_control to keep the top
+ # level clean and avoid clashing with the underwater transmission_loss.
+ "expansion_chamber",
+ "helmholtz_resonator",
+ "quarter_wave_resonator",
+ "extended_tube_chamber",
+ "ReactiveSilencerResult",
+ "end_reflection_loss",
+ "elbow_insertion_loss",
+ "plenum_attenuation",
+ "flow_noise_straight_duct",
+ "flow_noise_bend",
+ "HvacSpectrumResult",
+ "enclosure_insertion_loss",
+ "EnclosureResult",
]
diff --git a/src/phonometry/_plot/electroacoustics.py b/src/phonometry/_plot/electroacoustics.py
index 304f85b3..c5ac2b23 100644
--- a/src/phonometry/_plot/electroacoustics.py
+++ b/src/phonometry/_plot/electroacoustics.py
@@ -22,6 +22,7 @@
from matplotlib.axes import Axes
from ..electroacoustics.distortion import HarmonicDistortionResult
from ..electroacoustics.frequency_response import FrequencyResponseResult
+ from ..electroacoustics.piston import RadiatingPistonResult
from ..electroacoustics.swept_sine import SweptSineDistortionResult
#: Shared frequency-axis label of the electroacoustics renderers.
@@ -200,3 +201,34 @@ def _thd_panel(axt: Axes) -> None:
axes[1].set_xlabel("Excitation frequency [Hz]")
format_frequency_axis(axes[1])
return axes
+
+
+def plot_piston_impedance(
+ result: "RadiatingPistonResult", ax: Axes | None = None, **kwargs: Any
+) -> Axes:
+ """Normalized radiation resistance and reactance of a baffled piston.
+
+ Draws the piston resistance ``R1(2ka)`` and reactance ``X1(2ka)`` against
+ ``ka`` on a logarithmic ``ka`` axis (the classic Beranek & Mellow figure):
+ ``R1`` rising to 1, ``X1`` peaking near ``ka ~ 0.5`` then decaying.
+
+ :param result: A :class:`~phonometry.electroacoustics.piston.RadiatingPistonResult`.
+ :param ax: Existing axes, or ``None`` to create a figure.
+ :param kwargs: Forwarded to the resistance ``Axes.plot`` (its primary curve).
+ :return: The axes.
+ """
+ ax = ax if ax is not None else _new_axes()
+ ka = np.asarray(result.ka, dtype=np.float64)
+ kwargs.setdefault("color", _C_PRIMARY)
+ kwargs.setdefault("label", r"$R_1$ (resistance)")
+ ax.semilogx(ka, np.asarray(result.resistance), lw=1.8, **kwargs)
+ ax.semilogx(ka, np.asarray(result.reactance), color=_C_SECONDARY, lw=1.8,
+ ls="--", label=r"$X_1$ (reactance)")
+ ax.axhline(1.0, color=_C_TERTIARY, ls=":", lw=1.0,
+ label=r"$R_1 \to 1$ ($ka \gg 1$)")
+ ax.set_xlabel(r"$ka$")
+ ax.set_ylabel(r"Normalized radiation impedance $Z_r / \rho c S$")
+ ax.set_title("Baffled circular piston radiation impedance")
+ ax.grid(True, which="both", alpha=0.3)
+ ax.legend(loc="best", fontsize="small")
+ return ax
diff --git a/src/phonometry/_plot/noise_control.py b/src/phonometry/_plot/noise_control.py
new file mode 100644
index 00000000..1179142a
--- /dev/null
+++ b/src/phonometry/_plot/noise_control.py
@@ -0,0 +1,134 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Plot renderers for the noise_control domain (lazy imports from result .plot())."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+
+from .common import (
+ _C_PRIMARY,
+ _C_REFERENCE,
+ _C_SECONDARY,
+ _C_TERTIARY,
+ _new_axes,
+ format_frequency_axis,
+)
+
+if TYPE_CHECKING:
+ from matplotlib.axes import Axes
+
+ from ..noise_control.enclosures import EnclosureResult
+ from ..noise_control.hvac import HvacSpectrumResult
+ from ..noise_control.silencers import ReactiveSilencerResult
+
+_FREQ_LABEL = "Frequency [Hz]"
+
+
+def plot_reactive_silencer(
+ result: "ReactiveSilencerResult", ax: "Axes | None" = None, **kwargs: Any
+) -> "Axes":
+ """Transmission (and insertion) loss of a reactive silencer over frequency.
+
+ :param result: A
+ :class:`~phonometry.noise_control.silencers.ReactiveSilencerResult`.
+ :param ax: Existing axes, or ``None`` to create a figure.
+ :param kwargs: Forwarded to the transmission-loss ``Axes.plot``.
+ :return: The axes.
+ """
+ ax = ax if ax is not None else _new_axes()
+ f = np.asarray(result.frequencies, dtype=np.float64)
+ kwargs.setdefault("color", _C_PRIMARY)
+ kwargs.setdefault("label", "Transmission loss")
+ kwargs.setdefault("lw", 1.8)
+ ax.plot(f, np.asarray(result.transmission_loss), **kwargs)
+ if result.insertion_loss is not None:
+ ax.plot(f, np.asarray(result.insertion_loss), color=_C_SECONDARY,
+ lw=1.4, ls="--", label="Insertion loss")
+ if result.resonances is not None:
+ labeled = False
+ for fr in np.atleast_1d(np.asarray(result.resonances)):
+ if f.min() <= fr <= f.max():
+ ax.axvline(fr, color=_C_TERTIARY, ls=":", lw=1.0,
+ label="Resonance" if not labeled else None)
+ labeled = True
+ ax.set_xlabel(_FREQ_LABEL)
+ ax.set_ylabel("Loss [dB]")
+ ax.set_title(f"Reactive silencer: {result.kind}")
+ ax.grid(True, which="both", alpha=0.3)
+ format_frequency_axis(ax)
+ ax.legend(loc="best", fontsize="small")
+ return ax
+
+
+def plot_hvac_spectrum(
+ result: "HvacSpectrumResult", ax: "Axes | None" = None, **kwargs: Any
+) -> "Axes":
+ """Per-frequency HVAC attenuation or regenerated sound power level.
+
+ :param result: A :class:`~phonometry.noise_control.hvac.HvacSpectrumResult`.
+ :param ax: Existing axes, or ``None`` to create a figure.
+ :param kwargs: Forwarded to ``Axes.plot``.
+ :return: The axes.
+ """
+ ax = ax if ax is not None else _new_axes()
+ f = np.asarray(result.frequencies, dtype=np.float64)
+ is_power = result.quantity == "sound_power_level"
+ kwargs.setdefault("color", _C_SECONDARY if is_power else _C_PRIMARY)
+ kwargs.setdefault("label", result.label)
+ kwargs.setdefault("lw", 1.8)
+ kwargs.setdefault("marker", "o")
+ kwargs.setdefault("ms", 3)
+ ax.plot(f, np.asarray(result.values), **kwargs)
+ ax.set_xlabel(_FREQ_LABEL)
+ ax.set_ylabel(
+ "Sound power level [dB re 1 pW]" if is_power else "Attenuation [dB]"
+ )
+ ax.set_title(result.label)
+ ax.grid(True, which="both", alpha=0.3)
+ format_frequency_axis(ax)
+ ax.legend(loc="best", fontsize="small")
+ return ax
+
+
+def plot_enclosure(
+ result: "EnclosureResult", ax: "Axes | None" = None, **kwargs: Any
+) -> "Axes":
+ """Panel R, interior correction C and net insertion loss of an enclosure.
+
+ :param result: An
+ :class:`~phonometry.noise_control.enclosures.EnclosureResult`.
+ :param ax: Existing axes, or ``None`` to create a figure.
+ :param kwargs: Forwarded to the insertion-loss ``Axes.plot``.
+ :return: The axes.
+ """
+ ax = ax if ax is not None else _new_axes()
+ n = np.asarray(result.insertion_loss).size
+ if result.frequencies is not None:
+ x = np.asarray(result.frequencies, dtype=np.float64)
+ continuous = True
+ else:
+ x = np.arange(n, dtype=np.float64)
+ continuous = False
+ ax.plot(x, np.asarray(result.panel_transmission_loss), color=_C_REFERENCE,
+ lw=1.3, ls="--", marker="s", ms=3, label="Panel R")
+ ax.plot(x, np.asarray(result.correction), color=_C_TERTIARY, lw=1.3,
+ ls=":", marker="^", ms=3, label="Interior correction C")
+ kwargs.setdefault("color", _C_PRIMARY)
+ kwargs.setdefault("label", "Insertion loss (R - C)")
+ kwargs.setdefault("lw", 1.9)
+ kwargs.setdefault("marker", "o")
+ kwargs.setdefault("ms", 3)
+ ax.plot(x, np.asarray(result.insertion_loss), **kwargs)
+ ax.set_ylabel("Level [dB]")
+ ax.set_title("Machine enclosure insertion loss")
+ ax.grid(True, which="both", alpha=0.3)
+ if continuous:
+ ax.set_xlabel(_FREQ_LABEL)
+ format_frequency_axis(ax)
+ else:
+ ax.set_xlabel("Band")
+ ax.set_xticks(x)
+ ax.legend(loc="best", fontsize="small")
+ return ax
diff --git a/src/phonometry/electroacoustics/__init__.py b/src/phonometry/electroacoustics/__init__.py
index 02659bd7..96093d87 100644
--- a/src/phonometry/electroacoustics/__init__.py
+++ b/src/phonometry/electroacoustics/__init__.py
@@ -21,6 +21,13 @@
weighted_thd,
)
from .frequency_response import FrequencyResponseResult, coherence, transfer_function
+from .piston import (
+ RadiatingPistonResult,
+ piston_directivity,
+ piston_reactance,
+ piston_resistance,
+ radiating_piston,
+)
from .swept_sine import (
SweptSineDistortionResult,
swept_sine_distortion,
@@ -31,6 +38,7 @@
"FrequencyResponseResult",
"HarmonicDistortionResult",
"ModulationDistortionResult",
+ "RadiatingPistonResult",
"SweptSineDistortionResult",
"coherence",
"difference_frequency_distortion",
@@ -41,6 +49,10 @@
"idle_channel_noise",
"itu_r_468_weighting",
"modulation_distortion",
+ "piston_directivity",
+ "piston_reactance",
+ "piston_resistance",
+ "radiating_piston",
"sinad",
"swept_sine_distortion",
"synchronized_sweep_signal",
diff --git a/src/phonometry/electroacoustics/piston.py b/src/phonometry/electroacoustics/piston.py
new file mode 100644
index 00000000..55e92a75
--- /dev/null
+++ b/src/phonometry/electroacoustics/piston.py
@@ -0,0 +1,291 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""
+Radiation of a rigid circular piston set in an infinite baffle.
+
+The baffled circular piston is the canonical acoustic radiator: a flat rigid
+disc of radius ``a`` vibrating with a uniform normal velocity in an otherwise
+rigid infinite plane. It is the model behind a loudspeaker cone in a large
+cabinet, the open end of a duct, and the reference source for the radiation
+efficiency of any finite vibrating surface, so its two results -- the
+**radiation impedance** the air presents to the piston and the **directivity**
+of the far field -- are the base of the electroacoustics domain (Beranek &
+Mellow, *Acoustics: Sound Fields, Transducers and Vibration* 2nd ed., §4.4;
+Bies, Hansen & Howard, *Engineering Noise Control* 5th ed.).
+
+**Radiation impedance.** The reaction force of the air on the piston is
+``F = Z_r u`` with the mechanical radiation impedance
+
+ Z_r = rho c S ( R1(2ka) + j X1(2ka) ), S = pi a^2,
+
+where ``k = omega / c`` is the wavenumber, ``rho c`` the characteristic
+impedance of air and ``S`` the piston area. The dimensionless **piston
+resistance** and **reactance** functions are (Beranek & Mellow Eq. (4.30))
+
+ R1(x) = 1 - 2 J1(x) / x, X1(x) = 2 H1(x) / x,
+
+with ``J1`` the Bessel function of the first kind and ``H1`` the Struve
+function, both of order one, evaluated at ``x = 2ka``.
+
+* **Low frequency** (``ka << 1``): ``R1 -> (ka)^2 / 2`` so the radiated power
+ rises as ``f^2``, and ``X1 -> (8 / 3 pi) ka``. The reactance is mass-like,
+ ``X_r = rho c S X1 = omega M_r`` with the **radiation (accreted) mass**
+
+ M_r = 8 rho a^3 / 3
+
+ (Beranek & Mellow Eq. (4.32)): the piston drags an extra ``8 rho a^3 / 3`` of
+ air, equivalent to a layer ``8a / 3 pi`` thick over its face.
+* **High frequency** (``ka >> 1``): ``R1 -> 1`` and ``X1 -> 0``, so
+ ``Z_r -> rho c S`` -- the piston radiates as if into an infinite tube and the
+ air loads it purely resistively.
+
+**Directivity.** The far-field pressure of the baffled piston varies with the
+polar angle ``theta`` from the axis as (Beranek & Mellow Eq. (4.42))
+
+ D(theta) = 2 J1(ka sin theta) / (ka sin theta), D(0) = 1.
+
+The main lobe narrows as ``ka`` grows; its first null is at
+``ka sin theta = 3.8317`` (the first zero of ``J1``), which exists only once
+``ka > 3.8317``. The **directivity factor** ``Q`` (on-axis intensity over the
+intensity of a point source of equal power radiating into the full sphere) and
+the **directivity index** ``DI = 10 log10 Q`` follow from integrating
+``|D|^2`` over the radiating hemisphere,
+
+ Q = 2 / integral_0^(pi/2) |D(theta)|^2 sin theta d theta,
+
+which tends to ``Q = 2`` (``DI = 3.01 dB``, the half-space baffle gain) at low
+``ka`` and to ``Q ~ (ka)^2`` (``DI ~ 20 log10 ka``) at high ``ka``.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+from scipy import special
+
+from .._internal.validation import require_positive
+
+if TYPE_CHECKING:
+ from matplotlib.axes import Axes
+
+#: First zero of the Bessel function ``J1``; the first directivity null of the
+#: baffled piston sits at ``ka sin theta = J1_FIRST_ZERO`` (Beranek & Mellow).
+J1_FIRST_ZERO = 3.8317059702075125
+
+#: Reference characteristic impedance factors of air at 20 degC, 101.325 kPa.
+_RHO_AIR = 1.206
+_C_AIR = 343.0
+
+
+def piston_resistance(x: ArrayLike) -> np.ndarray | float:
+ """Piston resistance function ``R1(x) = 1 - 2 J1(x) / x``.
+
+ The real part of the normalized radiation impedance of a baffled circular
+ piston, as a function of ``x = 2ka`` (Beranek & Mellow Eq. (4.30)). It
+ rises as ``x^2 / 8 = (ka)^2 / 2`` at low ``x`` and tends to 1 at high ``x``.
+
+ :param x: Argument ``x = 2ka`` (scalar or array), dimensionless.
+ :return: ``R1(x)`` (float for scalar input, else an array).
+ """
+ arr = np.asarray(x, dtype=np.float64)
+ if np.any(arr < 0.0) or not np.all(np.isfinite(arr)):
+ raise ValueError("'x' must be non-negative and finite.")
+ with np.errstate(invalid="ignore", divide="ignore"):
+ out = 1.0 - 2.0 * special.j1(arr) / arr
+ # J1(x)/x -> 1/2 as x -> 0, so R1 -> 0; the input is validated finite, so
+ # the 0/0 at x = 0 is the only non-finite point and is set to that limit.
+ out = np.where(np.isfinite(out), out, 0.0)
+ return out[()] if out.ndim == 0 else out
+
+
+def piston_reactance(x: ArrayLike) -> np.ndarray | float:
+ """Piston reactance function ``X1(x) = 2 H1(x) / x`` (``H1`` Struve order 1).
+
+ The imaginary part of the normalized radiation impedance of a baffled
+ circular piston (Beranek & Mellow Eq. (4.30)). It rises as
+ ``(8 / 3 pi) ka`` (mass-like) at low ``x = 2ka`` and decays to 0 at high
+ ``x``.
+
+ :param x: Argument ``x = 2ka`` (scalar or array), dimensionless.
+ :return: ``X1(x)`` (float for scalar input, else an array).
+ """
+ arr = np.asarray(x, dtype=np.float64)
+ if np.any(arr < 0.0) or not np.all(np.isfinite(arr)):
+ raise ValueError("'x' must be non-negative and finite.")
+ with np.errstate(invalid="ignore", divide="ignore"):
+ out = 2.0 * special.struve(1, arr) / arr
+ # H1(x)/x -> 0 as x -> 0 (H1 ~ 2 x^2 / 3 pi), so X1 -> 0; the 0/0 at x = 0
+ # is the only non-finite point (the input is validated finite).
+ out = np.where(np.isfinite(out), out, 0.0)
+ return out[()] if out.ndim == 0 else out
+
+
+def piston_directivity(ka: ArrayLike, theta: ArrayLike) -> np.ndarray | float:
+ """Far-field directivity ``D = 2 J1(ka sin theta) / (ka sin theta)``.
+
+ The pressure amplitude of a baffled circular piston relative to its on-axis
+ value (Beranek & Mellow Eq. (4.42)), normalized so ``D(0) = 1``.
+
+ :param ka: Wavenumber-radius product ``ka`` (scalar or array).
+ :param theta: Polar angle from the axis, rad (scalar or array). Broadcast
+ against ``ka``.
+ :return: ``D`` (float for scalar inputs, else an array).
+ """
+ ka_arr = np.asarray(ka, dtype=np.float64)
+ theta_arr = np.asarray(theta, dtype=np.float64)
+ if np.any(ka_arr < 0.0) or not np.all(np.isfinite(ka_arr)):
+ raise ValueError("'ka' must be non-negative and finite.")
+ if not np.all(np.isfinite(theta_arr)):
+ raise ValueError("'theta' must be finite.")
+ u = ka_arr * np.sin(theta_arr)
+ with np.errstate(invalid="ignore", divide="ignore"):
+ out = 2.0 * special.j1(u) / u
+ # 2 J1(u)/u -> 1 as u -> 0 (on-axis, or ka = 0): the 0/0 there is the only
+ # non-finite point and is set to that limit.
+ out = np.where(np.isfinite(out), out, 1.0)
+ return out[()] if out.ndim == 0 else out
+
+
+def _directivity_index(ka: NDArray[np.float64]) -> NDArray[np.float64]:
+ """Directivity index ``DI = 10 log10 Q`` per ``ka`` by hemisphere quadrature.
+
+ ``Q = 2 / integral_0^(pi/2) |D|^2 sin theta d theta`` (radiation into the
+ half-space in front of the baffle). Uses a fine trapezoid grid; ``ka -> 0``
+ gives the half-space baffle gain ``DI = 10 log10 2 = 3.01 dB`` and high
+ ``ka`` tends to ``Q = (ka)^2`` (``DI = 20 log10 ka``).
+ """
+ theta = np.linspace(0.0, 0.5 * np.pi, 2001)
+ sin_t = np.sin(theta)
+ d = np.asarray(
+ piston_directivity(ka[:, None], theta[None, :]), dtype=np.float64
+ )
+ integrand = d**2 * sin_t[None, :]
+ integral = np.trapezoid(integrand, theta, axis=-1)
+ return np.asarray(10.0 * np.log10(2.0 / integral), dtype=np.float64)
+
+
+@dataclass(frozen=True)
+class RadiatingPistonResult:
+ """Radiation impedance and directivity of a baffled circular piston.
+
+ :ivar frequencies: Frequencies ``f``, Hz.
+ :ivar ka: Wavenumber-radius product ``ka`` at each frequency.
+ :ivar resistance: Normalized piston resistance ``R1(2ka)`` (real part of
+ ``Z_r / (rho c S)``).
+ :ivar reactance: Normalized piston reactance ``X1(2ka)`` (imaginary part of
+ ``Z_r / (rho c S)``).
+ :ivar radiation_resistance: Mechanical radiation resistance
+ ``rho c S R1``, N s/m.
+ :ivar radiation_reactance: Mechanical radiation reactance
+ ``rho c S X1``, N s/m.
+ :ivar radiation_mass: Low-frequency accreted air mass ``M_r = 8 rho a^3/3``,
+ kg (a single value; the mass limit of ``radiation_reactance / omega``).
+ :ivar directivity_index: Directivity index ``DI = 10 log10 Q``, dB.
+ :ivar angles: Polar angles of ``directivity``, rad, or ``None`` if not
+ requested.
+ :ivar directivity: Far-field directivity ``D(theta)`` as a
+ ``(n_freq, n_angle)`` array, or ``None`` if ``angles`` was not given.
+ :ivar radius: Piston radius ``a``, m.
+ :ivar speed_of_sound: Speed of sound ``c``, m/s.
+ :ivar density: Air density ``rho``, kg/m3.
+ """
+
+ frequencies: np.ndarray
+ ka: np.ndarray
+ resistance: np.ndarray
+ reactance: np.ndarray
+ radiation_resistance: np.ndarray
+ radiation_reactance: np.ndarray
+ radiation_mass: float
+ directivity_index: np.ndarray
+ angles: np.ndarray | None
+ directivity: np.ndarray | None
+ radius: float
+ speed_of_sound: float
+ density: float
+
+ def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes":
+ """Plot the normalized piston resistance and reactance against ``ka``.
+
+ Reproduces the classic Beranek & Mellow figure: ``R1`` rising to 1 and
+ ``X1`` peaking then decaying, over the ``ka`` range of the result.
+ Requires matplotlib (``pip install phonometry[plot]``).
+ """
+ from .._plot.electroacoustics import plot_piston_impedance
+
+ return plot_piston_impedance(self, ax=ax, **kwargs)
+
+
+def radiating_piston(
+ radius: float,
+ frequencies: ArrayLike,
+ *,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+ angles: ArrayLike | None = None,
+) -> RadiatingPistonResult:
+ """Radiation impedance and directivity of a rigid baffled circular piston.
+
+ Evaluates the piston resistance ``R1(2ka)`` and reactance ``X1(2ka)``, the
+ mechanical radiation impedance ``rho c S (R1 + j X1)``, the low-frequency
+ radiation mass ``8 rho a^3 / 3`` and the directivity index over the given
+ frequencies (Beranek & Mellow §4.4). Pass ``angles`` to also sample the
+ far-field directivity pattern ``D(theta)``.
+
+ :param radius: Piston radius ``a``, m.
+ :param frequencies: Frequencies ``f``, Hz (scalar or 1-D array), all > 0.
+ :param speed_of_sound: Speed of sound ``c``, m/s (default 343).
+ :param density: Air density ``rho``, kg/m3 (default 1.206).
+ :param angles: Optional polar angles ``theta`` from the axis, rad, at which
+ to sample the directivity pattern.
+ :return: A :class:`RadiatingPistonResult`.
+ """
+ a = require_positive(radius, "radius")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64))
+ if f.ndim != 1 or f.size == 0:
+ raise ValueError("'frequencies' must be a non-empty 1-D array.")
+ if np.any(f <= 0.0) or not np.all(np.isfinite(f)):
+ raise ValueError("'frequencies' must be positive and finite.")
+
+ omega = 2.0 * np.pi * f
+ k = omega / c
+ ka = k * a
+ area = np.pi * a**2
+ r1 = np.asarray(piston_resistance(2.0 * ka), dtype=np.float64)
+ x1 = np.asarray(piston_reactance(2.0 * ka), dtype=np.float64)
+ rho_c_s = rho * c * area
+ radiation_mass = 8.0 * rho * a**3 / 3.0
+ di = _directivity_index(ka)
+
+ directivity: np.ndarray | None = None
+ angle_arr: np.ndarray | None = None
+ if angles is not None:
+ angle_arr = np.atleast_1d(np.asarray(angles, dtype=np.float64))
+ if angle_arr.ndim != 1 or angle_arr.size == 0:
+ raise ValueError("'angles' must be a non-empty 1-D array.")
+ if not np.all(np.isfinite(angle_arr)):
+ raise ValueError("'angles' must be finite.")
+ directivity = np.asarray(
+ piston_directivity(ka[:, None], angle_arr[None, :]),
+ dtype=np.float64,
+ )
+
+ return RadiatingPistonResult(
+ frequencies=f,
+ ka=ka,
+ resistance=r1,
+ reactance=x1,
+ radiation_resistance=rho_c_s * r1,
+ radiation_reactance=rho_c_s * x1,
+ radiation_mass=radiation_mass,
+ directivity_index=di,
+ angles=angle_arr,
+ directivity=directivity,
+ radius=a,
+ speed_of_sound=c,
+ density=rho,
+ )
diff --git a/src/phonometry/noise_control/__init__.py b/src/phonometry/noise_control/__init__.py
new file mode 100644
index 00000000..ba13487d
--- /dev/null
+++ b/src/phonometry/noise_control/__init__.py
@@ -0,0 +1,51 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""noise_control domain of phonometry (see module docstrings)."""
+
+from __future__ import annotations
+
+from .enclosures import EnclosureResult, enclosure_insertion_loss
+from .hvac import (
+ HvacSpectrumResult,
+ elbow_insertion_loss,
+ end_reflection_loss,
+ flow_noise_bend,
+ flow_noise_straight_duct,
+ plenum_attenuation,
+)
+from .silencers import (
+ ReactiveSilencerResult,
+ cascade,
+ duct_matrix,
+ expansion_chamber,
+ extended_tube_chamber,
+ helmholtz_impedance,
+ helmholtz_resonator,
+ insertion_loss,
+ quarter_wave_impedance,
+ quarter_wave_resonator,
+ shunt_matrix,
+ transmission_loss,
+)
+
+__all__ = [
+ "EnclosureResult",
+ "HvacSpectrumResult",
+ "ReactiveSilencerResult",
+ "cascade",
+ "duct_matrix",
+ "elbow_insertion_loss",
+ "enclosure_insertion_loss",
+ "end_reflection_loss",
+ "expansion_chamber",
+ "extended_tube_chamber",
+ "flow_noise_bend",
+ "flow_noise_straight_duct",
+ "helmholtz_impedance",
+ "helmholtz_resonator",
+ "insertion_loss",
+ "plenum_attenuation",
+ "quarter_wave_impedance",
+ "quarter_wave_resonator",
+ "shunt_matrix",
+ "transmission_loss",
+]
diff --git a/src/phonometry/noise_control/enclosures.py b/src/phonometry/noise_control/enclosures.py
new file mode 100644
index 00000000..574926e3
--- /dev/null
+++ b/src/phonometry/noise_control/enclosures.py
@@ -0,0 +1,173 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""
+Insertion loss of a close or free-standing machine enclosure.
+
+Wrapping a machine in a sealed enclosure reduces the radiated noise by the
+transmission loss of its panels, *minus* a penalty for the reverberant build-up
+inside the small, hard cavity. Bies, Hansen & Howard, *Engineering Noise
+Control* 5th ed., §7.4.2 (Eqs. (7.103), (7.111)) write the net reduction as
+
+ IL = R - C, C = 10 log10[ 0.3 + S_E (1 - alpha_i) / (S_i alpha_i) ],
+
+where ``R`` is the field-incidence transmission loss of the enclosure panels,
+``S_E`` the external surface area, ``S_i`` the internal surface area (including
+the machine) and ``alpha_i`` the mean absorption of the enclosure interior. The
+reverberant term is exactly ``S_E`` over the interior **room constant**
+``R_i = S_i alpha_i / (1 - alpha_i)`` (:func:`phonometry.room.room_constant`), so
+
+ C = 10 log10( 0.3 + S_E / R_i ).
+
+A hard interior (``alpha_i`` small) makes ``C`` large and wastes much of the
+panel ``R``; lining the enclosure drives ``C`` toward its floor
+``10 log10 0.3 = -5.2 dB`` (a fully absorbing interior, where ``IL = R + 5.2``).
+Bies terms this net reduction the enclosure *noise reduction*; it is the
+insertion loss of the enclosure.
+
+**The panel transmission loss ``R`` is supplied by the caller** -- measured, or
+predicted by a panel model -- as a per-band array or a callable of frequency.
+This module never predicts ``R`` itself; it combines a given ``R`` with the
+interior absorption. The interior room constant reuses
+:func:`phonometry.room.room_constant`.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Callable
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from .._internal.validation import require_positive
+from ..room.steady_field import room_constant
+
+if TYPE_CHECKING:
+ from matplotlib.axes import Axes
+
+
+@dataclass(frozen=True)
+class EnclosureResult:
+ """Insertion loss of a machine enclosure over frequency (Bies §7.4.2).
+
+ :ivar frequencies: Frequencies ``f``, Hz, or ``None`` if the panel ``R``
+ was given as a bare per-band array with no frequency labels.
+ :ivar panel_transmission_loss: The supplied panel transmission loss ``R``
+ per band, dB.
+ :ivar correction: The interior-build-up correction ``C`` per band, dB.
+ :ivar insertion_loss: The net enclosure insertion loss ``IL = R - C``, dB.
+ :ivar external_area: External enclosure surface area ``S_E``, m2.
+ :ivar internal_area: Internal surface area ``S_i``, m2.
+ :ivar room_constant: Interior room constant ``R_i`` per band, m2.
+ """
+
+ frequencies: np.ndarray | None
+ panel_transmission_loss: np.ndarray
+ correction: np.ndarray
+ insertion_loss: np.ndarray
+ external_area: float
+ internal_area: float
+ room_constant: np.ndarray
+
+ def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes":
+ """Plot the panel ``R``, correction ``C`` and net insertion loss.
+
+ Requires matplotlib (``pip install phonometry[plot]``); returns the
+ :class:`~matplotlib.axes.Axes`.
+ """
+ from .._plot.noise_control import plot_enclosure
+
+ return plot_enclosure(self, ax=ax, **kwargs)
+
+
+def _resolve_frequencies(
+ frequencies: ArrayLike | None,
+) -> NDArray[np.float64] | None:
+ """Validate an optional 1-D positive frequency grid (Hz)."""
+ if frequencies is None:
+ return None
+ freqs = np.atleast_1d(np.asarray(frequencies, dtype=np.float64))
+ if freqs.ndim != 1 or freqs.size == 0:
+ raise ValueError("'frequencies' must be a non-empty 1-D array.")
+ if np.any(freqs <= 0.0) or not np.all(np.isfinite(freqs)):
+ raise ValueError("'frequencies' must be positive and finite.")
+ return freqs
+
+
+def _resolve_panel_r(
+ panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike],
+ freqs: NDArray[np.float64] | None,
+) -> NDArray[np.float64]:
+ """Resolve the panel transmission loss ``R`` into a validated 1-D array."""
+ if callable(panel_transmission_loss):
+ if freqs is None:
+ raise ValueError(
+ "'frequencies' is required when 'panel_transmission_loss' "
+ "is a callable."
+ )
+ r = np.atleast_1d(np.asarray(panel_transmission_loss(freqs), dtype=np.float64))
+ else:
+ r = np.atleast_1d(np.asarray(panel_transmission_loss, dtype=np.float64))
+ if r.ndim != 1 or r.size == 0:
+ raise ValueError("'panel_transmission_loss' must be a non-empty 1-D array.")
+ if not np.all(np.isfinite(r)):
+ raise ValueError("'panel_transmission_loss' must be finite.")
+ return r
+
+
+def enclosure_insertion_loss(
+ panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike],
+ external_area: float,
+ internal_area: float,
+ internal_absorption: ArrayLike,
+ *,
+ frequencies: ArrayLike | None = None,
+) -> EnclosureResult:
+ """Net insertion loss of a machine enclosure (Bies Eqs. (7.103), (7.111)).
+
+ ``IL = R - C`` with ``C = 10 log10(0.3 + S_E / R_i)`` and the interior room
+ constant ``R_i = S_i alpha_i / (1 - alpha_i)``.
+
+ :param panel_transmission_loss: Panel transmission loss ``R`` per band, dB.
+ Either a per-band array (measured or predicted elsewhere) or a callable
+ mapping a frequency array to per-band ``R`` (then ``frequencies`` is
+ required). This function does not predict ``R``.
+ :param external_area: External enclosure surface area ``S_E``, m2.
+ :param internal_area: Internal surface area ``S_i`` (including the machine),
+ m2.
+ :param internal_absorption: Mean interior absorption ``alpha_i`` in
+ ``(0, 1)`` (scalar or per-band).
+ :param frequencies: Band centre frequencies, Hz; required when
+ ``panel_transmission_loss`` is a callable, optional otherwise (used to
+ label the result and the plot).
+ :return: An :class:`EnclosureResult`.
+ """
+ s_e = require_positive(external_area, "external_area")
+ s_i = require_positive(internal_area, "internal_area")
+
+ freqs = _resolve_frequencies(frequencies)
+ r = _resolve_panel_r(panel_transmission_loss, freqs)
+
+ alpha = np.asarray(internal_absorption, dtype=np.float64)
+ if alpha.ndim > 1:
+ raise ValueError("'internal_absorption' must be a scalar or a 1-D array.")
+ if np.any(alpha <= 0.0) or np.any(alpha >= 1.0) or not np.all(np.isfinite(alpha)):
+ raise ValueError("'internal_absorption' must lie strictly in (0, 1).")
+
+ r_i = np.atleast_1d(np.asarray(room_constant(s_i, alpha), dtype=np.float64))
+ r_i_b, r_b = np.broadcast_arrays(r_i, r)
+ if freqs is not None and freqs.shape != r_b.shape:
+ raise ValueError(
+ "'frequencies' must match the number of panel-R / absorption bands."
+ )
+ correction = 10.0 * np.log10(0.3 + s_e / r_i_b)
+ il = r_b - correction
+ return EnclosureResult(
+ frequencies=freqs,
+ panel_transmission_loss=np.array(r_b, dtype=np.float64),
+ correction=np.array(correction, dtype=np.float64),
+ insertion_loss=np.array(il, dtype=np.float64),
+ external_area=s_e,
+ internal_area=s_i,
+ room_constant=np.array(r_i_b, dtype=np.float64),
+ )
diff --git a/src/phonometry/noise_control/hvac.py b/src/phonometry/noise_control/hvac.py
new file mode 100644
index 00000000..9717a5c8
--- /dev/null
+++ b/src/phonometry/noise_control/hvac.py
@@ -0,0 +1,347 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""
+HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise.
+
+A ventilation duct network attenuates fan noise through several mechanisms
+that add up along the path, and it *regenerates* noise wherever the airflow is
+disturbed. This module gathers the engineering methods of Bies, Hansen &
+Howard, *Engineering Noise Control* 5th ed., Chapter 8, for the passive
+attenuations -- **duct end reflection** (§8.13, Table 8.14), **bends/elbows**
+(§8.11, Table 8.11) and **plenum chambers** (§8.17, Wells' method) -- and for
+the **flow-generated (self) noise** of straight ducts and bends (§8.15).
+
+The end-reflection and elbow methods are empirical look-up tables (ASHRAE);
+they are interpolated over the duct size and, for the elbows, over the
+frequency-to-width ratio ``W / lambda``. The plenum and flow-noise methods are
+closed forms evaluated directly.
+
+.. note::
+ Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14
+ look-up (there is no closed form in this edition); this module reproduces
+ that table and interpolates it. Rectangular ducts use the equivalent
+ diameter ``D = sqrt(4 S / pi)``.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from .._internal.validation import require_positive
+from ..room.steady_field import room_constant
+
+if TYPE_CHECKING:
+ from matplotlib.axes import Axes
+
+_C_AIR = 343.0
+
+# ---------------------------------------------------------------------------
+# Bies Table 8.14 -- duct end reflection loss (dB), ASHRAE.
+# Rows: internal diameter (mm). Columns: octave band centre (Hz).
+# Two termination conditions: "flush" (duct flush with a wall/ceiling) and
+# "free" (free space / suspended in the room).
+# ---------------------------------------------------------------------------
+_END_REFLECTION_BANDS: NDArray[np.float64] = np.array(
+ [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0]
+)
+_END_REFLECTION_DIAMETERS_MM: NDArray[np.float64] = np.array(
+ [150, 200, 250, 300, 400, 510, 610, 710, 810, 910, 1220, 1830], dtype=float
+)
+_END_REFLECTION_FLUSH: NDArray[np.float64] = np.array(
+ [
+ [18, 12, 7, 3, 1, 0],
+ [15, 10, 5, 2, 1, 0],
+ [14, 8, 4, 1, 0, 0],
+ [12, 7, 3, 1, 0, 0],
+ [10, 5, 2, 1, 0, 0],
+ [8, 4, 1, 0, 0, 0],
+ [7, 3, 1, 0, 0, 0],
+ [6, 2, 1, 0, 0, 0],
+ [5, 2, 1, 0, 0, 0],
+ [4, 2, 0, 0, 0, 0],
+ [3, 1, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0, 0],
+ ],
+ dtype=float,
+)
+_END_REFLECTION_FREE: NDArray[np.float64] = np.array(
+ [
+ [20, 14, 9, 5, 2, 1],
+ [18, 12, 7, 3, 1, 0],
+ [16, 11, 6, 2, 1, 0],
+ [14, 9, 5, 2, 1, 0],
+ [12, 7, 3, 1, 0, 0],
+ [10, 6, 2, 1, 0, 0],
+ [9, 5, 2, 1, 0, 0],
+ [8, 4, 1, 0, 0, 0],
+ [7, 3, 1, 0, 0, 0],
+ [6, 3, 1, 0, 0, 0],
+ [5, 2, 0, 0, 0, 0],
+ [3, 1, 0, 0, 0, 0],
+ ],
+ dtype=float,
+)
+
+# ---------------------------------------------------------------------------
+# Bies Table 8.11 -- elbow/bend insertion loss (dB per bend) vs W / lambda.
+# The five columns are the supported (bend_type, vanes, lined, round) cases;
+# each row is the value for the W/lambda band whose upper edge is the key.
+# ---------------------------------------------------------------------------
+_ELBOW_WL_UPPER: NDArray[np.float64] = np.array(
+ [0.14, 0.28, 0.55, 1.11, 2.22, np.inf]
+)
+_ELBOW_TABLE: dict[str, NDArray[np.float64]] = {
+ # square, no vanes, unlined
+ "square": np.array([0, 1, 5, 8, 4, 3], dtype=float),
+ # square, no vanes, lined
+ "square_lined": np.array([0, 1, 6, 11, 10, 10], dtype=float),
+ # square, with vanes, unlined
+ "square_vanes": np.array([0, 1, 4, 6, 4, 4], dtype=float),
+ # square, with vanes, lined
+ "square_vanes_lined": np.array([0, 1, 4, 7, 7, 7], dtype=float),
+ # round, no vanes, unlined
+ "round": np.array([0, 1, 2, 3, 3, 3], dtype=float),
+}
+
+
+def _frequencies(frequencies: ArrayLike) -> NDArray[np.float64]:
+ f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64))
+ if f.ndim != 1 or f.size == 0:
+ raise ValueError("'frequencies' must be a non-empty 1-D array.")
+ if np.any(f <= 0.0) or not np.all(np.isfinite(f)):
+ raise ValueError("'frequencies' must be positive and finite.")
+ return f
+
+
+@dataclass(frozen=True)
+class HvacSpectrumResult:
+ """A per-frequency HVAC quantity (attenuation or regenerated power level).
+
+ :ivar frequencies: Frequencies ``f``, Hz.
+ :ivar values: The quantity per frequency (dB, or dB re 1e-12 W for a
+ sound power level).
+ :ivar quantity: What ``values`` holds (``"attenuation"`` or
+ ``"sound_power_level"``).
+ :ivar label: A short human label of the element.
+ """
+
+ frequencies: np.ndarray
+ values: np.ndarray
+ quantity: str
+ label: str
+
+ def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes":
+ """Plot the quantity against a continuous log-frequency axis.
+
+ Requires matplotlib (``pip install phonometry[plot]``).
+ """
+ from .._plot.noise_control import plot_hvac_spectrum
+
+ return plot_hvac_spectrum(self, ax=ax, **kwargs)
+
+
+def end_reflection_loss(
+ frequencies: ArrayLike,
+ diameter: float,
+ *,
+ termination: str = "flush",
+ speed_of_sound: float = _C_AIR,
+) -> HvacSpectrumResult:
+ """Duct end reflection loss (Bies Table 8.14, ASHRAE).
+
+ The low-frequency reflection of sound back up a duct at its open
+ termination into a room. Interpolated over ``log`` diameter and ``log``
+ frequency from Table 8.14; it passes exactly through the tabulated
+ ``(diameter, octave band)`` nodes.
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param diameter: Duct internal diameter ``D``, m (use
+ ``D = sqrt(4 S / pi)`` for a rectangular duct of area ``S``).
+ :param termination: ``"flush"`` (duct flush with a wall/ceiling) or
+ ``"free"`` (free space / suspended in the room).
+ :param speed_of_sound: Speed of sound ``c``, m/s (kept for signature
+ symmetry; the table is indexed by frequency directly).
+ :return: A :class:`HvacSpectrumResult` of the reflection loss, dB.
+ """
+ f = _frequencies(frequencies)
+ d_mm = require_positive(diameter, "diameter") * 1000.0
+ if termination == "flush":
+ table = _END_REFLECTION_FLUSH
+ elif termination == "free":
+ table = _END_REFLECTION_FREE
+ else:
+ raise ValueError("'termination' must be 'flush' or 'free'.")
+ require_positive(speed_of_sound, "speed_of_sound")
+
+ log_d = np.log(_END_REFLECTION_DIAMETERS_MM)
+ log_f_band = np.log(_END_REFLECTION_BANDS)
+ # Interpolate the table in log-diameter for each band, then in log-freq.
+ per_band = np.array(
+ [np.interp(np.log(d_mm), log_d, table[:, j]) for j in range(table.shape[1])]
+ )
+ values = np.interp(np.log(f), log_f_band, per_band)
+ return HvacSpectrumResult(
+ frequencies=f, values=values, quantity="attenuation",
+ label=f"End reflection ({termination}, D = {diameter * 1000:.0f} mm)",
+ )
+
+
+def elbow_insertion_loss(
+ frequencies: ArrayLike,
+ width: float,
+ *,
+ bend_type: str = "square",
+ vanes: bool = False,
+ lined: bool = False,
+ speed_of_sound: float = _C_AIR,
+) -> HvacSpectrumResult:
+ """Duct bend/elbow insertion loss per bend (Bies Table 8.11, ASHRAE).
+
+ Indexed by the frequency-to-width ratio ``W / lambda`` (``lambda = c / f``).
+ Lined bends assume the lining extends at least three duct diameters up- and
+ downstream. Round bends are treated as unlined with no vanes.
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param width: Duct width ``W`` in the plane of the bend, m.
+ :param bend_type: ``"square"`` or ``"round"``.
+ :param vanes: Turning vanes fitted (square bends only).
+ :param lined: Acoustically lined bend (square bends only).
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :return: A :class:`HvacSpectrumResult` of the insertion loss, dB per bend.
+ """
+ f = _frequencies(frequencies)
+ w = require_positive(width, "width")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ if bend_type == "round":
+ if vanes or lined:
+ raise ValueError("round bends take neither vanes nor lining.")
+ key = "round"
+ elif bend_type == "square":
+ key = "square" + ("_vanes" if vanes else "") + ("_lined" if lined else "")
+ else:
+ raise ValueError("'bend_type' must be 'square' or 'round'.")
+ col = _ELBOW_TABLE[key]
+ wl = w * f / c
+ idx = np.searchsorted(_ELBOW_WL_UPPER, wl, side="left")
+ idx = np.clip(idx, 0, col.size - 1)
+ values = col[idx]
+ return HvacSpectrumResult(
+ frequencies=f, values=values, quantity="attenuation",
+ label=f"Elbow ({key.replace('_', ', ')}, W = {w * 1000:.0f} mm)",
+ )
+
+
+def plenum_attenuation(
+ exit_area: float,
+ line_of_sight: float,
+ wall_area: float,
+ mean_absorption: ArrayLike,
+ *,
+ angle: float = 0.0,
+) -> np.ndarray | float:
+ """Plenum-chamber transmission loss by Wells' method (Bies Eq. (8.275)).
+
+ ``TL = -10 log10[ S_out ( cos(theta) / (pi r^2) + (1 - alpha) / (S_w alpha) ) ]``,
+ where the reverberant term uses the plenum room constant
+ ``R = S_w alpha / (1 - alpha)`` (:func:`phonometry.room.room_constant`). The
+ method holds above the inlet cut-on and when the plenum is large compared
+ with the wavelength; it underpredicts the low-frequency loss by 5-10 dB.
+
+ :param exit_area: Outlet-opening area ``S_out``, m2.
+ :param line_of_sight: Straight-line inlet-to-outlet distance ``r``, m.
+ :param wall_area: Total internal wall area ``S_w``, m2.
+ :param mean_absorption: Mean Sabine wall absorption ``alpha`` in ``(0, 1)``
+ (scalar or per-band).
+ :param angle: Angle ``theta`` between the inlet axis and the line to the
+ outlet, rad (default 0).
+ :return: The transmission loss, dB (float for scalar absorption, else a
+ per-band array).
+ """
+ s_out = require_positive(exit_area, "exit_area")
+ r = require_positive(line_of_sight, "line_of_sight")
+ s_w = require_positive(wall_area, "wall_area")
+ alpha = np.asarray(mean_absorption, dtype=np.float64)
+ if np.any(alpha <= 0.0) or np.any(alpha >= 1.0) or not np.all(np.isfinite(alpha)):
+ raise ValueError("'mean_absorption' must lie strictly in (0, 1).")
+ r_const = np.asarray(room_constant(s_w, alpha), dtype=np.float64)
+ direct = np.cos(angle) / (np.pi * r**2)
+ reverberant = 1.0 / r_const
+ tl = -10.0 * np.log10(s_out * (direct + reverberant))
+ return float(tl) if tl.ndim == 0 else tl
+
+
+def flow_noise_straight_duct(
+ frequencies: ArrayLike,
+ flow_velocity: float,
+ area: float,
+) -> HvacSpectrumResult:
+ """Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)).
+
+ ``L_WB = 7 + 50 log10(U) + 10 log10(S) - 2 - 26 log10(1.14 + 0.02 f / U)``
+ in dB re 1e-12 W (VDI 2081-1), for airflow speed ``U`` in a duct of area
+ ``S``.
+
+ :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array).
+ :param flow_velocity: Mean flow speed ``U``, m/s.
+ :param area: Duct cross-sectional area ``S``, m2.
+ :return: A :class:`HvacSpectrumResult` of the band sound power level,
+ dB re 1e-12 W.
+ """
+ f = _frequencies(frequencies)
+ u = require_positive(flow_velocity, "flow_velocity")
+ s = require_positive(area, "area")
+ lw = (
+ 7.0
+ + 50.0 * np.log10(u)
+ + 10.0 * np.log10(s)
+ - 2.0
+ - 26.0 * np.log10(1.14 + 0.02 * f / u)
+ )
+ return HvacSpectrumResult(
+ frequencies=f, values=lw, quantity="sound_power_level",
+ label=f"Straight-duct flow noise (U = {u:.1f} m/s)",
+ )
+
+
+def flow_noise_bend(
+ frequencies: ArrayLike,
+ flow_velocity: float,
+ area: float,
+ height: float,
+ *,
+ density: float = 1.206,
+) -> HvacSpectrumResult:
+ """Flow-generated octave-band sound power of a mitred bend (Bies Eqs. (8.252), (8.254)).
+
+ ``L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103`` with the
+ stream power level ``L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117``
+ (Bies Eq. (8.252)) and the Strouhal number ``N_s = f H / U`` (``H`` the duct
+ height in the plane of the bend). The radiated sound power grows as the
+ sixth power of the stream speed at low ``N_s`` (the inner-corner drag
+ dipole) and the eighth power at high ``N_s`` (the outer-corner shear
+ quadrupole); equivalently, the *efficiency* referenced to the stream power
+ grows as ``U^3`` and ``U^5`` respectively.
+
+ :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array).
+ :param flow_velocity: Mean flow speed ``U``, m/s.
+ :param area: Duct cross-sectional area ``S``, m2.
+ :param height: Duct height ``H`` in the plane of the bend, m.
+ :param density: Air density ``rho``, kg/m3.
+ :return: A :class:`HvacSpectrumResult` of the band sound power level,
+ dB re 1e-12 W.
+ """
+ f = _frequencies(frequencies)
+ u = require_positive(flow_velocity, "flow_velocity")
+ s = require_positive(area, "area")
+ h = require_positive(height, "height")
+ rho = require_positive(density, "density")
+ lws = 30.0 * np.log10(u) + 10.0 * np.log10(s) + 10.0 * np.log10(rho) + 117.0
+ ns = f * h / u
+ lw = lws - 10.0 * np.log10(1.0 + 0.165 * ns**2) + 30.0 * np.log10(u) - 103.0
+ return HvacSpectrumResult(
+ frequencies=f, values=lw, quantity="sound_power_level",
+ label=f"Mitred-bend flow noise (U = {u:.1f} m/s)",
+ )
diff --git a/src/phonometry/noise_control/silencers.py b/src/phonometry/noise_control/silencers.py
new file mode 100644
index 00000000..5244c679
--- /dev/null
+++ b/src/phonometry/noise_control/silencers.py
@@ -0,0 +1,579 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""
+Reactive silencers by the four-pole (transmission-matrix) method.
+
+A reactive silencer controls noise by *reflecting* it back to the source with
+impedance discontinuities -- sudden area changes and side branches -- rather
+than by dissipating it in absorptive material. The one-dimensional plane-wave
+theory represents each acoustic element by a 2x2 **transfer (four-pole)
+matrix** relating the sound pressure ``p`` and volume velocity ``S u`` at its
+two ends, and a compound silencer is the ordered matrix product of its
+elements (Bies, Hansen & Howard, *Engineering Noise Control* 5th ed., §8.8-8.9;
+Munjal, *Acoustics of Ducts and Mufflers*).
+
+**Transfer matrix** (Bies Eq. (8.133)), state vector ``[p, S u]`` with the
+characteristic acoustic impedance ``Z = rho c / S``. The plane-wave element for
+a straight duct of length ``L`` and area ``S`` is (Bies Eq. (8.143), no flow)
+
+ [[ cos(kL), j (rho c / S) sin(kL) ],
+ [ j (S / rho c) sin(kL), cos(kL) ]], k = omega / c,
+
+and a **side branch** of acoustic impedance ``Z_b`` is the shunt element
+(Bies Eq. (8.144))
+
+ [[ 1, 0 ],
+ [ 1 / Z_b, 1 ]].
+
+**Transmission loss** from the compound matrix ``T`` (Bies Eq. (8.141), no
+flow; reduces to Eq. (8.148) for equal inlet/outlet areas):
+
+ TL = 10 log10[ (Z1 / Zn) (1/4) | T11 + T12 / Zn + Z1 T21 + (Z1 / Zn) T22 |^2 ]
+
+with ``Z1 = rho c / S_in`` and ``Zn = rho c / S_out``. ``TL`` is the intrinsic
+attenuation for an anechoic termination. The **insertion loss** for a source of
+internal impedance ``Z_s`` radiating into a termination impedance ``Z_r`` is
+the extra attenuation of inserting the silencer in place of a direct
+connection,
+
+ IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) / (Z_s + Z_r) |,
+
+which is ``0`` when the silencer reduces to a through connection (``T = I``)
+and equals the transmission loss for the anechoic reference
+``Z_s = rho c / S_in``, ``Z_r = rho c / S_out``.
+
+**Simple expansion chamber.** A chamber of area ``S_exp`` and length ``L``
+between pipes of area ``S_duct`` has the closed-form transmission loss (Bies
+Eq. (8.111)) with area ratio ``m = S_exp / S_duct``
+
+ TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ],
+
+peaking at ``10 log10[1 + (1/4)(m - 1/m)^2]`` when ``kL = pi/2, 3pi/2, ...`` and
+dropping to ``0`` at ``kL = n pi`` (no dissipation). The four-pole product
+reproduces this exactly, and the machinery extends to side-branch (Helmholtz,
+quarter-wave) and extended-tube resonators that the closed form cannot cover.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from .._internal.validation import require_non_negative, require_positive
+
+if TYPE_CHECKING:
+ from matplotlib.axes import Axes
+
+#: Reference air properties at 20 degC, 101.325 kPa.
+_C_AIR = 343.0
+_RHO_AIR = 1.206
+
+_Complex = NDArray[np.complex128]
+
+
+def _frequencies(frequencies: ArrayLike) -> NDArray[np.float64]:
+ """Validate a strictly positive, finite 1-D frequency grid (Hz)."""
+ f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64))
+ if f.ndim != 1 or f.size == 0:
+ raise ValueError("'frequencies' must be a non-empty 1-D array.")
+ if np.any(f <= 0.0) or not np.all(np.isfinite(f)):
+ raise ValueError("'frequencies' must be positive and finite.")
+ return f
+
+
+def duct_matrix(
+ frequencies: ArrayLike,
+ length: float,
+ area: float,
+ *,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+) -> _Complex:
+ """Four-pole matrix of a straight duct (Bies Eq. (8.143), no flow).
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param length: Duct length ``L``, m.
+ :param area: Cross-sectional area ``S``, m2.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :return: A ``(n_freq, 2, 2)`` complex transfer-matrix array.
+ """
+ f = _frequencies(frequencies)
+ length = require_non_negative(length, "length")
+ area = require_positive(area, "area")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ k = 2.0 * np.pi * f / c
+ z = rho * c / area
+ cos = np.cos(k * length)
+ sin = np.sin(k * length)
+ t = np.empty((f.size, 2, 2), dtype=np.complex128)
+ t[:, 0, 0] = cos
+ t[:, 0, 1] = 1j * z * sin
+ t[:, 1, 0] = 1j * sin / z
+ t[:, 1, 1] = cos
+ return t
+
+
+def shunt_matrix(branch_impedance: ArrayLike) -> _Complex:
+ """Four-pole matrix of a side branch of impedance ``Z_b`` (Bies Eq. (8.144)).
+
+ :param branch_impedance: Acoustic impedance ``Z_b`` of the branch,
+ Pa s/m3 (1-D complex array over frequency).
+ :return: A ``(n_freq, 2, 2)`` complex transfer-matrix array.
+ """
+ zb = np.atleast_1d(np.asarray(branch_impedance, dtype=np.complex128))
+ if zb.ndim != 1 or zb.size == 0:
+ raise ValueError("'branch_impedance' must be a non-empty 1-D array.")
+ t = np.zeros((zb.size, 2, 2), dtype=np.complex128)
+ t[:, 0, 0] = 1.0
+ # A lossless branch at exact resonance has Z_b -> 0 (it shorts the duct);
+ # 1/Z_b -> infinity gives the ideal infinite attenuation there.
+ with np.errstate(divide="ignore", invalid="ignore"):
+ t[:, 1, 0] = 1.0 / zb
+ t[:, 1, 1] = 1.0
+ return t
+
+
+def cascade(*matrices: _Complex) -> _Complex:
+ """Cascade element four-pole matrices from inlet to outlet.
+
+ The compound matrix is the ordered product ``T1 @ T2 @ ... @ Tn`` (the
+ state at the inlet equals the compound matrix times the state at the
+ outlet), broadcast over the frequency axis.
+
+ :param matrices: One or more ``(n_freq, 2, 2)`` arrays sharing ``n_freq``.
+ :return: The compound ``(n_freq, 2, 2)`` array.
+ """
+ if not matrices:
+ raise ValueError("cascade() needs at least one matrix.")
+ n = matrices[0].shape[0]
+ if any(m.shape[0] != n for m in matrices[1:]):
+ raise ValueError(
+ "cascade() matrices must share the same frequency grid (n_freq)."
+ )
+ total = matrices[0]
+ for m in matrices[1:]:
+ total = np.matmul(total, m)
+ return total
+
+
+def transmission_loss(
+ transfer_matrix: _Complex,
+ *,
+ inlet_area: float,
+ outlet_area: float,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+) -> NDArray[np.float64]:
+ """Transmission loss of a four-pole element (Bies Eq. (8.141), no flow).
+
+ :param transfer_matrix: A ``(n_freq, 2, 2)`` compound matrix.
+ :param inlet_area: Inlet pipe area ``S_in``, m2.
+ :param outlet_area: Outlet pipe area ``S_out``, m2.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :return: The transmission loss per frequency, dB.
+ """
+ s_in = require_positive(inlet_area, "inlet_area")
+ s_out = require_positive(outlet_area, "outlet_area")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ z1 = rho * c / s_in
+ zn = rho * c / s_out
+ t11 = transfer_matrix[:, 0, 0]
+ t12 = transfer_matrix[:, 0, 1]
+ t21 = transfer_matrix[:, 1, 0]
+ t22 = transfer_matrix[:, 1, 1]
+ term = t11 + t12 / zn + z1 * t21 + (z1 / zn) * t22
+ return np.asarray(
+ 10.0 * np.log10((z1 / zn) * 0.25 * np.abs(term) ** 2),
+ dtype=np.float64,
+ )
+
+
+def insertion_loss(
+ transfer_matrix: _Complex,
+ *,
+ source_impedance: ArrayLike,
+ radiation_impedance: ArrayLike,
+) -> NDArray[np.float64]:
+ """Insertion loss of a four-pole element for given end impedances.
+
+ The attenuation from inserting the element in place of a direct (zero
+ length) connection between a source of internal impedance ``Z_s`` and a
+ radiation (termination) impedance ``Z_r``:
+
+ IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) /
+ (Z_s + Z_r) |.
+
+ :param transfer_matrix: A ``(n_freq, 2, 2)`` compound matrix.
+ :param source_impedance: Source internal acoustic impedance ``Z_s``,
+ Pa s/m3 (scalar or per-frequency, real or complex).
+ :param radiation_impedance: Termination/radiation acoustic impedance
+ ``Z_r``, Pa s/m3 (scalar or per-frequency).
+ :return: The insertion loss per frequency, dB.
+ """
+ n = transfer_matrix.shape[0]
+ zs = np.broadcast_to(
+ np.asarray(source_impedance, dtype=np.complex128), (n,)
+ )
+ zr = np.broadcast_to(
+ np.asarray(radiation_impedance, dtype=np.complex128), (n,)
+ )
+ t11 = transfer_matrix[:, 0, 0]
+ t12 = transfer_matrix[:, 0, 1]
+ t21 = transfer_matrix[:, 1, 0]
+ t22 = transfer_matrix[:, 1, 1]
+ num = t11 * zr + t12 + zs * zr * t21 + zs * t22
+ return np.asarray(
+ 20.0 * np.log10(np.abs(num / (zs + zr))), dtype=np.float64
+ )
+
+
+def helmholtz_impedance(
+ frequencies: ArrayLike,
+ neck_area: float,
+ neck_length: float,
+ cavity_volume: float,
+ *,
+ resistance: float = 0.0,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+) -> _Complex:
+ """Acoustic impedance of a Helmholtz side branch (Bies Eq. (8.152)).
+
+ ``Z = R + j( rho omega l_e / S_neck - rho c^2 / (omega V) )`` with acoustic
+ mass ``rho l_e / S_neck`` and compliance ``V / (rho c^2)``; the resonance
+ ``f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))`` (Bies Eq. (8.46)) is where the
+ reactance vanishes and the branch shorts the duct.
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param neck_area: Neck cross-sectional area ``S_neck``, m2.
+ :param neck_length: Effective neck length ``l_e`` (with end corrections), m.
+ :param cavity_volume: Cavity volume ``V``, m3.
+ :param resistance: Acoustic resistance ``R``, Pa s/m3 (default 0, lossless).
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :return: The complex branch impedance per frequency, Pa s/m3.
+ """
+ f = _frequencies(frequencies)
+ s_neck = require_positive(neck_area, "neck_area")
+ le = require_positive(neck_length, "neck_length")
+ vol = require_positive(cavity_volume, "cavity_volume")
+ r = require_non_negative(resistance, "resistance")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ omega = 2.0 * np.pi * f
+ reactance = rho * omega * le / s_neck - rho * c**2 / (omega * vol)
+ return np.asarray(r + 1j * reactance, dtype=np.complex128)
+
+
+def quarter_wave_impedance(
+ frequencies: ArrayLike,
+ length: float,
+ area: float,
+ *,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+) -> _Complex:
+ """Acoustic impedance of a closed quarter-wave side branch (Bies Eq. (8.146)).
+
+ ``Z = -j (rho c / S) cot(k l_e)``; the reactance vanishes at
+ ``l_e = lambda / 4`` (``f = c / 4 l_e``), where the closed tube presents a
+ pressure node and shorts the duct.
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param length: Effective tube length ``l_e`` (with end correction), m.
+ :param area: Tube cross-sectional area ``S``, m2.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :return: The complex branch impedance per frequency, Pa s/m3.
+ """
+ f = _frequencies(frequencies)
+ le = require_positive(length, "length")
+ area = require_positive(area, "area")
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ k = 2.0 * np.pi * f / c
+ z = rho * c / area
+ # At the half-wave frequencies (k l_e = n pi) the closed tube is transparent
+ # (Z -> infinity); the divide-by-zero is expected and its shunt 1/Z is 0.
+ with np.errstate(divide="ignore", invalid="ignore"):
+ zb = -1j * z / np.tan(k * le)
+ return np.asarray(zb, dtype=np.complex128)
+
+
+@dataclass(frozen=True)
+class ReactiveSilencerResult:
+ """Transmission and insertion loss of a reactive silencer over frequency.
+
+ :ivar frequencies: Frequencies ``f``, Hz.
+ :ivar transmission_loss: Transmission loss per frequency, dB.
+ :ivar insertion_loss: Insertion loss per frequency, dB, or ``None`` when no
+ source/radiation impedance was supplied.
+ :ivar transfer_matrix: The compound ``(n_freq, 2, 2)`` four-pole matrix.
+ :ivar kind: A short label of the device (e.g. ``"expansion chamber"``).
+ :ivar resonances: Notable resonance frequencies, Hz (e.g. the resonator
+ tuning frequency), or ``None``.
+ """
+
+ frequencies: np.ndarray
+ transmission_loss: np.ndarray
+ insertion_loss: np.ndarray | None
+ transfer_matrix: np.ndarray
+ kind: str
+ resonances: np.ndarray | None = None
+
+ def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes":
+ """Plot the transmission (and insertion) loss against frequency.
+
+ Requires matplotlib (``pip install phonometry[plot]``); returns the
+ :class:`~matplotlib.axes.Axes`.
+ """
+ from .._plot.noise_control import plot_reactive_silencer
+
+ return plot_reactive_silencer(self, ax=ax, **kwargs)
+
+
+def _result(
+ f: NDArray[np.float64],
+ t: _Complex,
+ *,
+ inlet_area: float,
+ outlet_area: float,
+ c: float,
+ rho: float,
+ source_impedance: ArrayLike | None,
+ radiation_impedance: ArrayLike | None,
+ kind: str,
+ resonances: NDArray[np.float64] | None = None,
+) -> ReactiveSilencerResult:
+ """Assemble a :class:`ReactiveSilencerResult` from a compound matrix."""
+ tl = transmission_loss(
+ t, inlet_area=inlet_area, outlet_area=outlet_area,
+ speed_of_sound=c, density=rho,
+ )
+ il: NDArray[np.float64] | None = None
+ if source_impedance is not None and radiation_impedance is not None:
+ il = insertion_loss(
+ t, source_impedance=source_impedance,
+ radiation_impedance=radiation_impedance,
+ )
+ return ReactiveSilencerResult(
+ frequencies=f,
+ transmission_loss=tl,
+ insertion_loss=il,
+ transfer_matrix=t,
+ kind=kind,
+ resonances=resonances,
+ )
+
+
+def expansion_chamber(
+ frequencies: ArrayLike,
+ length: float,
+ chamber_area: float,
+ pipe_area: float,
+ *,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult:
+ """Simple expansion-chamber silencer (Bies Eq. (8.111) / four-pole).
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param length: Chamber length ``L``, m.
+ :param chamber_area: Chamber cross-sectional area ``S_exp``, m2.
+ :param pipe_area: Inlet/outlet pipe area ``S_duct``, m2.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :param source_impedance: Optional source impedance ``Z_s`` for the
+ insertion loss, Pa s/m3.
+ :param radiation_impedance: Optional radiation impedance ``Z_r`` for the
+ insertion loss, Pa s/m3.
+ :return: A :class:`ReactiveSilencerResult` (its ``transmission_loss``
+ equals the closed form ``10 log10[1 + (1/4)(m - 1/m)^2 sin^2(kL)]``).
+ """
+ f = _frequencies(frequencies)
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ require_positive(chamber_area, "chamber_area")
+ require_positive(pipe_area, "pipe_area")
+ t = duct_matrix(f, length, chamber_area, speed_of_sound=c, density=rho)
+ return _result(
+ f, t, inlet_area=pipe_area, outlet_area=pipe_area, c=c, rho=rho,
+ source_impedance=source_impedance,
+ radiation_impedance=radiation_impedance, kind="expansion chamber",
+ )
+
+
+def helmholtz_resonator(
+ frequencies: ArrayLike,
+ duct_area: float,
+ neck_area: float,
+ neck_length: float,
+ cavity_volume: float,
+ *,
+ resistance: float = 0.0,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult:
+ """Side-branch Helmholtz resonator on a duct (Bies Eqs. (8.144), (8.152)).
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param duct_area: Main-duct cross-sectional area ``S_d``, m2.
+ :param neck_area: Resonator neck area ``S_neck``, m2.
+ :param neck_length: Effective neck length ``l_e``, m.
+ :param cavity_volume: Cavity volume ``V``, m3.
+ :param resistance: Neck acoustic resistance ``R``, Pa s/m3 (default 0).
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3.
+ :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3.
+ :return: A :class:`ReactiveSilencerResult`; ``resonances`` holds
+ ``f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))``.
+ """
+ f = _frequencies(frequencies)
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ s_d = require_positive(duct_area, "duct_area")
+ zb = helmholtz_impedance(
+ f, neck_area, neck_length, cavity_volume, resistance=resistance,
+ speed_of_sound=c, density=rho,
+ )
+ t = shunt_matrix(zb)
+ f0 = c / (2.0 * np.pi) * np.sqrt(neck_area / (neck_length * cavity_volume))
+ return _result(
+ f, t, inlet_area=s_d, outlet_area=s_d, c=c, rho=rho,
+ source_impedance=source_impedance,
+ radiation_impedance=radiation_impedance,
+ kind="Helmholtz resonator",
+ resonances=np.array([f0], dtype=np.float64),
+ )
+
+
+def quarter_wave_resonator(
+ frequencies: ArrayLike,
+ duct_area: float,
+ length: float,
+ branch_area: float,
+ *,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult:
+ """Closed quarter-wave side-branch tube on a duct (Bies Eqs. (8.144), (8.146)).
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param duct_area: Main-duct cross-sectional area ``S_d``, m2.
+ :param length: Effective branch length ``l_e``, m.
+ :param branch_area: Branch tube area ``S``, m2.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3.
+ :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3.
+ :return: A :class:`ReactiveSilencerResult`; ``resonances`` holds the odd
+ multiples of ``f = c / (4 l_e)`` within the frequency range.
+ """
+ f = _frequencies(frequencies)
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ s_d = require_positive(duct_area, "duct_area")
+ le = require_positive(length, "length")
+ zb = quarter_wave_impedance(
+ f, le, branch_area, speed_of_sound=c, density=rho
+ )
+ t = shunt_matrix(zb)
+ f_fundamental = c / (4.0 * le)
+ odds = np.arange(1, int(2.0 * f.max() / f_fundamental) + 2, 2)
+ res = np.asarray(odds * f_fundamental, dtype=np.float64)
+ res = res[res <= f.max()]
+ if res.size == 0:
+ res = np.array([f_fundamental], dtype=np.float64)
+ return _result(
+ f, t, inlet_area=s_d, outlet_area=s_d, c=c, rho=rho,
+ source_impedance=source_impedance,
+ radiation_impedance=radiation_impedance,
+ kind="quarter-wave resonator",
+ resonances=res,
+ )
+
+
+def extended_tube_chamber(
+ frequencies: ArrayLike,
+ length: float,
+ chamber_area: float,
+ pipe_area: float,
+ *,
+ inlet_extension: float = 0.0,
+ outlet_extension: float = 0.0,
+ speed_of_sound: float = _C_AIR,
+ density: float = _RHO_AIR,
+ source_impedance: ArrayLike | None = None,
+ radiation_impedance: ArrayLike | None = None,
+) -> ReactiveSilencerResult:
+ """Extended-inlet/outlet expansion chamber (Bies §8.9.7).
+
+ The inlet and outlet pipes extend a distance into the chamber, forming
+ annular quarter-wave side branches (of area ``S_exp - S_duct`` and lengths
+ equal to the extensions, Bies Eq. (8.156)) at the two junctions. Tuning the
+ extensions (classically ``L/4`` and ``L/2``) places quarter-wave peaks that
+ fill the ``kL = n pi`` troughs of the plain expansion chamber. With both
+ extensions ``0`` the result reduces exactly to :func:`expansion_chamber`.
+
+ :param frequencies: Frequencies ``f``, Hz (1-D array).
+ :param length: Chamber length ``L``, m.
+ :param chamber_area: Chamber cross-sectional area ``S_exp``, m2.
+ :param pipe_area: Inlet/outlet pipe area ``S_duct``, m2.
+ :param inlet_extension: Inlet pipe extension into the chamber ``L_a``, m.
+ :param outlet_extension: Outlet pipe extension into the chamber ``L_b``, m.
+ :param speed_of_sound: Speed of sound ``c``, m/s.
+ :param density: Air density ``rho``, kg/m3.
+ :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3.
+ :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3.
+ :return: A :class:`ReactiveSilencerResult`.
+ """
+ f = _frequencies(frequencies)
+ c = require_positive(speed_of_sound, "speed_of_sound")
+ rho = require_positive(density, "density")
+ s_exp = require_positive(chamber_area, "chamber_area")
+ s_duct = require_positive(pipe_area, "pipe_area")
+ la = require_non_negative(inlet_extension, "inlet_extension")
+ lb = require_non_negative(outlet_extension, "outlet_extension")
+ if s_exp <= s_duct:
+ raise ValueError("'chamber_area' must exceed 'pipe_area'.")
+ if la + lb > length:
+ raise ValueError(
+ "the inlet and outlet extensions cannot together exceed the "
+ "chamber length (they would overlap inside the chamber)."
+ )
+ annulus = s_exp - s_duct
+
+ elements = []
+ resonances = []
+ if la > 0.0:
+ elements.append(shunt_matrix(
+ quarter_wave_impedance(f, la, annulus, speed_of_sound=c, density=rho)
+ ))
+ resonances.append(c / (4.0 * la))
+ elements.append(duct_matrix(f, length, s_exp, speed_of_sound=c, density=rho))
+ if lb > 0.0:
+ elements.append(shunt_matrix(
+ quarter_wave_impedance(f, lb, annulus, speed_of_sound=c, density=rho)
+ ))
+ resonances.append(c / (4.0 * lb))
+ t = cascade(*elements)
+ res = np.array(resonances, dtype=np.float64) if resonances else None
+ return _result(
+ f, t, inlet_area=s_duct, outlet_area=s_duct, c=c, rho=rho,
+ source_impedance=source_impedance,
+ radiation_impedance=radiation_impedance,
+ kind="extended-tube chamber", resonances=res,
+ )
diff --git a/tests/electroacoustics/test_piston.py b/tests/electroacoustics/test_piston.py
new file mode 100644
index 00000000..e2bcef6f
--- /dev/null
+++ b/tests/electroacoustics/test_piston.py
@@ -0,0 +1,127 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Tests for the baffled circular piston (Beranek & Mellow 2e, §4.19 / §13.7).
+
+Oracles are the exact closed forms (``R1 = 1 - 2 J1(x)/x``, ``X1 = 2 H1(x)/x``
+with ``x = 2ka``), their low- and high-frequency limits (Eqs. (13.117),
+(13.118), radiation mass Eq. (4.151)), the first directivity null at the first
+zero of ``J1`` (3.8317), and the half-space baffle directivity index of 3.01 dB
+at ``ka -> 0``. Point values were cross-checked against scipy in the source
+extraction.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+from scipy import special
+
+from phonometry import (
+ RadiatingPistonResult,
+ piston_directivity,
+ piston_reactance,
+ piston_resistance,
+ radiating_piston,
+)
+from phonometry.electroacoustics.piston import J1_FIRST_ZERO
+
+
+def test_resistance_reactance_closed_form() -> None:
+ x = np.array([0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0, 10.0])
+ assert np.allclose(piston_resistance(x), 1.0 - 2.0 * special.j1(x) / x)
+ assert np.allclose(piston_reactance(x), 2.0 * special.struve(1, x) / x)
+
+
+def test_digitized_reference_points() -> None:
+ # Points verified from Eqs. (13.117)/(13.118) in the source extraction.
+ ref = {
+ 0.1: (0.001249, 0.042413),
+ 1.0: (0.119899, 0.396915),
+ 3.0: (0.773961, 0.680073),
+ 5.0: (1.131032, 0.323125),
+ 10.0: (0.991305, 0.178366),
+ }
+ for x, (r1, x1) in ref.items():
+ assert piston_resistance(x) == pytest.approx(r1, abs=1e-5)
+ assert piston_reactance(x) == pytest.approx(x1, abs=1e-5)
+
+
+def test_low_frequency_limits() -> None:
+ # R1 -> (ka)^2/2 = x^2/8 and X1 -> (8/3pi) ka = (4/3pi) x as x -> 0.
+ x = 1e-3
+ assert piston_resistance(x) == pytest.approx(x**2 / 8.0, rel=1e-4)
+ assert piston_reactance(x) == pytest.approx(4.0 / (3.0 * math.pi) * x, rel=1e-4)
+
+
+def test_high_frequency_limit_resistive() -> None:
+ # R1 -> 1 and X1 -> 0 at large x: Z_r -> rho c S (purely resistive).
+ assert piston_resistance(80.0) == pytest.approx(1.0, abs=0.02)
+ assert abs(piston_reactance(80.0)) < 0.05
+
+
+def test_zero_argument_regular() -> None:
+ assert piston_resistance(0.0) == 0.0
+ assert piston_reactance(0.0) == 0.0
+ assert not np.any(np.isnan(piston_resistance(np.array([0.0, 1.0]))))
+
+
+def test_directivity_onaxis_and_first_null() -> None:
+ assert piston_directivity(5.0, 0.0) == pytest.approx(1.0)
+ ka = 6.0
+ theta_null = math.asin(J1_FIRST_ZERO / ka)
+ assert abs(piston_directivity(ka, theta_null)) < 1e-9
+ # No null exists when ka < first zero of J1.
+ theta = np.linspace(0.0, math.pi / 2, 200)
+ assert np.min(np.abs(piston_directivity(3.0, theta))) > 0.1
+
+
+def test_radiation_mass_coefficient() -> None:
+ a, rho = 0.12, 1.206
+ res = radiating_piston(a, [100.0], density=rho)
+ assert res.radiation_mass == pytest.approx(8.0 * rho * a**3 / 3.0)
+
+
+def test_radiation_reactance_tends_to_mass_times_omega() -> None:
+ # At low ka the mechanical reactance rho c S X1 = omega M_r.
+ a = 0.05
+ res = radiating_piston(a, [20.0], density=1.206, speed_of_sound=343.0)
+ omega = 2.0 * math.pi * 20.0
+ assert res.radiation_reactance[0] == pytest.approx(
+ omega * res.radiation_mass, rel=2e-3
+ )
+
+
+def test_directivity_index_half_space_limit() -> None:
+ # DI -> 10 log10(2) = 3.01 dB as ka -> 0 (radiation into the half-space).
+ res = radiating_piston(0.01, [1.0])
+ assert res.directivity_index[0] == pytest.approx(3.0103, abs=1e-3)
+
+
+def test_directivity_index_high_ka() -> None:
+ # DI -> 20 log10(ka) at high ka (Q = (ka)^2).
+ a, f, c = 1.0, 5000.0, 343.0
+ res = radiating_piston(a, [f], speed_of_sound=c)
+ ka = 2.0 * math.pi * f * a / c
+ assert res.directivity_index[0] == pytest.approx(20.0 * math.log10(ka), abs=0.3)
+
+
+def test_result_shapes_and_plot() -> None:
+ res = radiating_piston(0.1, [100.0, 1000.0], angles=[0.0, 0.3, 0.6])
+ assert isinstance(res, RadiatingPistonResult)
+ assert res.directivity is not None and res.directivity.shape == (2, 3)
+ assert res.ka.shape == (2,)
+ import matplotlib
+
+ matplotlib.use("Agg")
+ ax = res.plot()
+ assert ax.get_ylabel()
+
+
+def test_validation() -> None:
+ with pytest.raises(ValueError):
+ radiating_piston(-1.0, [100.0])
+ with pytest.raises(ValueError):
+ radiating_piston(0.1, [0.0])
+ with pytest.raises(ValueError):
+ piston_resistance(-1.0)
diff --git a/tests/noise_control/test_enclosures.py b/tests/noise_control/test_enclosures.py
new file mode 100644
index 00000000..439abaa3
--- /dev/null
+++ b/tests/noise_control/test_enclosures.py
@@ -0,0 +1,94 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Tests for machine-enclosure insertion loss (Bies §7.4.2).
+
+Oracles: the closed form ``IL = R - C`` with
+``C = 10 log10(0.3 + S_E / R_i)`` and ``R_i`` the interior room constant
+(Eqs. (7.103), (7.111)); the fully absorbing limit ``C -> 10 log10 0.3``
+(``IL = R + 5.23``); consistency of ``R_i`` with
+:func:`phonometry.room.room_constant`; and the panel ``R`` being supplied by
+the caller (array or callable), never predicted here.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from phonometry import enclosure_insertion_loss, room_constant
+from phonometry.noise_control.enclosures import EnclosureResult
+
+
+def test_closed_form() -> None:
+ r = np.array([20.0, 30.0, 40.0])
+ s_e, s_i, alpha = 6.0, 5.0, 0.3
+ res = enclosure_insertion_loss(r, s_e, s_i, alpha)
+ r_i = float(room_constant(s_i, alpha))
+ c = 10.0 * math.log10(0.3 + s_e / r_i)
+ assert np.allclose(res.correction, c)
+ assert np.allclose(res.insertion_loss, r - c)
+
+
+def test_room_constant_consistency() -> None:
+ res = enclosure_insertion_loss([30.0], 6.0, 5.0, np.array([0.3]))
+ assert res.room_constant[0] == pytest.approx(float(room_constant(5.0, 0.3)))
+
+
+def test_fully_absorbing_interior_floor() -> None:
+ # alpha -> 1: R_i -> inf, C -> 10 log10(0.3) = -5.23 dB, IL = R + 5.23.
+ res = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.999999)
+ assert res.correction[0] == pytest.approx(10.0 * math.log10(0.3), abs=1e-3)
+ assert res.insertion_loss[0] == pytest.approx(40.0 - 10.0 * math.log10(0.3), abs=1e-3)
+
+
+def test_harder_interior_wastes_more_panel() -> None:
+ live = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.05)
+ dead = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.5)
+ assert dead.insertion_loss[0] > live.insertion_loss[0]
+
+
+def test_per_band_absorption() -> None:
+ r = np.array([20.0, 25.0, 30.0, 35.0])
+ alpha = np.array([0.1, 0.2, 0.3, 0.4])
+ res = enclosure_insertion_loss(r, 6.0, 5.0, alpha)
+ assert res.insertion_loss.shape == (4,)
+ r_i = np.asarray(room_constant(5.0, alpha))
+ assert np.allclose(res.correction, 10.0 * np.log10(0.3 + 6.0 / r_i))
+
+
+def test_callable_panel_r() -> None:
+ freqs = np.array([125.0, 250.0, 500.0, 1000.0])
+
+ def mass_law(f: np.ndarray) -> np.ndarray:
+ return 20.0 * np.log10(f) - 10.0 # arbitrary user-supplied R(f)
+
+ res = enclosure_insertion_loss(mass_law, 6.0, 5.0, 0.3, frequencies=freqs)
+ assert np.allclose(res.panel_transmission_loss, mass_law(freqs))
+ assert res.frequencies is not None
+
+
+def test_callable_requires_frequencies() -> None:
+ with pytest.raises(ValueError, match="frequencies"):
+ enclosure_insertion_loss(lambda f: f, 6.0, 5.0, 0.3)
+
+
+def test_result_type_and_plot() -> None:
+ res = enclosure_insertion_loss(
+ [20.0, 30.0, 40.0], 6.0, 5.0, 0.3,
+ frequencies=[125.0, 250.0, 500.0],
+ )
+ assert isinstance(res, EnclosureResult)
+ import matplotlib
+
+ matplotlib.use("Agg")
+ assert res.plot().get_ylabel()
+ # No-frequency plot uses band indices.
+ assert enclosure_insertion_loss([20.0, 30.0], 6.0, 5.0, 0.3).plot() is not None
+
+
+def test_validation() -> None:
+ with pytest.raises(ValueError):
+ enclosure_insertion_loss([20.0], -1.0, 5.0, 0.3)
+ with pytest.raises(ValueError):
+ enclosure_insertion_loss([20.0], 6.0, 5.0, 1.0)
diff --git a/tests/noise_control/test_fdtd_crosscheck.py b/tests/noise_control/test_fdtd_crosscheck.py
new file mode 100644
index 00000000..aeacb03c
--- /dev/null
+++ b/tests/noise_control/test_fdtd_crosscheck.py
@@ -0,0 +1,96 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Independent FDTD cross-check of a single 2D expansion chamber.
+
+The four-pole expansion-chamber transmission loss (Bies Eq. (8.111)) predicts,
+for an area ratio ``m = S_exp / S_duct``, a ``TL`` peak of
+``10 log10[1 + (1/4)(m - 1/m)^2]`` at ``kL = pi/2`` (``f = c / 4L``) and a
+``TL = 0`` trough (fully transparent) at ``kL = pi`` (``f = c / 2L``).
+
+This test drives an *independent* solver -- the 2D FDTD engine of
+:mod:`phonometry.simulation` -- with a plane-wave duct that widens into a
+chamber and narrows back, and compares the transmitted amplitude at the peak
+and trough frequencies. The measured ``20 log10`` ratio of the transmitted
+amplitudes (trough over peak) must reproduce the closed-form peak ``TL``. The
+run is deterministic and small (a 16x200 grid, ~3000 steps, well under a
+second).
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from phonometry.noise_control import silencers as sl
+from phonometry.simulation.fdtd import CWSource, fdtd_simulation
+
+_C = 343.0
+_DX = 0.02
+_NY, _NX = 16, 200
+_DUCT_ROWS = (6, 7, 8, 9) # 4-cell duct; chamber is the full 16 rows -> m = 4
+_X1, _X2 = 80, 100 # chamber columns
+_L = (_X2 - _X1) * _DX # chamber length, 0.4 m
+_M = _NY / len(_DUCT_ROWS) # area (height) ratio = 4
+
+
+def _chamber_mask() -> np.ndarray:
+ mask = np.ones((_NY, _NX), dtype=bool)
+ for r in _DUCT_ROWS:
+ mask[r, :] = False # open the narrow inlet/outlet duct
+ mask[:, _X1:_X2] = False # open the full-height chamber
+ return mask
+
+
+def _transmitted_rms(frequency: float) -> float:
+ mask = _chamber_mask()
+ boundaries = {
+ "left": "absorbing", "right": "absorbing",
+ "top": "rigid", "bottom": "rigid",
+ }
+ src = CWSource(ix=25, iy=7, frequency=frequency, amplitude=1.0, ramp_cycles=6.0)
+ steps = 3000
+ dt = 0.6 * _DX / (_C * math.sqrt(2.0))
+ res = fdtd_simulation(
+ _C, _DX, steps * dt, sources=[src], shape=(_NY, _NX),
+ probes=[(170, 7)], boundaries=boundaries,
+ absorbing_layer_cells=10, obstacle_mask=mask,
+ )
+ p = res.pressures[0]
+ tail = p[int(p.size * 0.6):] # steady-state window
+ return float(np.sqrt(np.mean(tail**2)))
+
+
+def test_fdtd_matches_expansion_chamber_peak_tl() -> None:
+ f_peak = _C / (4.0 * _L) # kL = pi/2, four-pole TL maximum
+ f_trough = _C / (2.0 * _L) # kL = pi, four-pole TL = 0 (transparent)
+
+ down_peak = _transmitted_rms(f_peak)
+ down_trough = _transmitted_rms(f_trough)
+
+ # The chamber transmits more at the trough than at the peak.
+ assert down_trough > down_peak
+
+ # The measured attenuation difference reproduces the closed-form peak TL.
+ measured_delta = 20.0 * math.log10(down_trough / down_peak)
+ peak_tl = 10.0 * math.log10(1.0 + 0.25 * (_M - 1.0 / _M) ** 2)
+ assert peak_tl == pytest.approx(6.55, abs=0.05)
+ assert measured_delta == pytest.approx(peak_tl, abs=1.2)
+
+
+def test_fdtd_run_is_deterministic() -> None:
+ # Two independent runs of the (RNG-free) solver must be bit-identical.
+ f = _C / (4.0 * _L)
+ first = _transmitted_rms(f)
+ second = _transmitted_rms(f)
+ assert first == second
+
+
+def test_four_pole_peak_frequency_and_trough() -> None:
+ # The four-pole model itself: TL peak at f = c/4L, TL = 0 at f = c/2L.
+ f = np.array([_C / (4.0 * _L), _C / (2.0 * _L)])
+ res = sl.expansion_chamber(f, _L, _M * 0.01, 0.01)
+ assert res.transmission_loss[0] == pytest.approx(
+ 10.0 * math.log10(1.0 + 0.25 * (_M - 1.0 / _M) ** 2), abs=1e-6
+ )
+ assert res.transmission_loss[1] == pytest.approx(0.0, abs=1e-9)
diff --git a/tests/noise_control/test_hvac.py b/tests/noise_control/test_hvac.py
new file mode 100644
index 00000000..7c447ce2
--- /dev/null
+++ b/tests/noise_control/test_hvac.py
@@ -0,0 +1,129 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Tests for the HVAC duct methods (Bies Chapter 8).
+
+Oracles: exact points of the ASHRAE end-reflection table (Bies Table 8.14) and
+elbow table (Bies Table 8.11), Wells' plenum closed form (Eq. (8.275)) with its
+room-constant reverberant term and limits, and the VDI 2081 flow-noise formulas
+(Eqs. (8.251), (8.254)) evaluated directly.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from phonometry.noise_control import hvac
+
+
+def test_end_reflection_table_nodes_flush() -> None:
+ # Bies Table 8.14, flush: D = 200 mm -> [15, 10, 5, 2, 1, 0] dB.
+ bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0]
+ res = hvac.end_reflection_loss(bands, 0.200, termination="flush")
+ assert np.allclose(res.values, [15, 10, 5, 2, 1, 0], atol=1e-6)
+
+
+def test_end_reflection_table_nodes_free() -> None:
+ # Bies Table 8.14, free space: D = 150 mm -> [20, 14, 9, 5, 2, 1] dB.
+ bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0]
+ res = hvac.end_reflection_loss(bands, 0.150, termination="free")
+ assert np.allclose(res.values, [20, 14, 9, 5, 2, 1], atol=1e-6)
+
+
+def test_end_reflection_free_exceeds_flush() -> None:
+ res_flush = hvac.end_reflection_loss([125.0], 0.3, termination="flush")
+ res_free = hvac.end_reflection_loss([125.0], 0.3, termination="free")
+ assert res_free.values[0] > res_flush.values[0]
+
+
+def test_elbow_table_bands() -> None:
+ # Square, no vanes, unlined W = 0.3 m; W/lambda selects the table step.
+ c = 343.0
+ # f such that W/lambda = 0.4 (in 0.28-0.55) -> 5 dB
+ f = 0.4 * c / 0.3
+ res = hvac.elbow_insertion_loss([f], 0.3, bend_type="square")
+ assert res.values[0] == pytest.approx(5.0)
+ # W/lambda < 0.14 -> 0 dB
+ res0 = hvac.elbow_insertion_loss([0.1 * c / 0.3], 0.3)
+ assert res0.values[0] == 0.0
+
+
+def test_elbow_lined_beats_unlined_high_freq() -> None:
+ c = 343.0
+ f = 1.5 * c / 0.3 # W/lambda = 1.5 (1.11-2.22 band)
+ unlined = hvac.elbow_insertion_loss([f], 0.3, bend_type="square").values[0]
+ lined = hvac.elbow_insertion_loss([f], 0.3, bend_type="square", lined=True).values[0]
+ assert lined > unlined # 10 vs 4 dB
+
+
+def test_elbow_round_rejects_options() -> None:
+ with pytest.raises(ValueError):
+ hvac.elbow_insertion_loss([500.0], 0.3, bend_type="round", lined=True)
+
+
+def test_plenum_wells_closed_form() -> None:
+ # TL = -10 log10[S_out(cos theta/(pi r^2) + (1-alpha)/(S_w alpha))].
+ s_out, r, s_w, alpha = 0.1, 1.0, 20.0, 0.2
+ tl = hvac.plenum_attenuation(s_out, r, s_w, alpha)
+ direct = 1.0 / (math.pi * r**2)
+ reverb = (1.0 - alpha) / (s_w * alpha)
+ expected = -10.0 * math.log10(s_out * (direct + reverb))
+ assert tl == pytest.approx(expected)
+
+
+def test_plenum_more_absorption_more_loss() -> None:
+ a = hvac.plenum_attenuation(0.1, 1.0, 20.0, 0.1)
+ b = hvac.plenum_attenuation(0.1, 1.0, 20.0, 0.5)
+ assert b > a
+
+
+def test_plenum_per_band() -> None:
+ tl = hvac.plenum_attenuation(0.1, 1.0, 20.0, np.array([0.1, 0.3, 0.5]))
+ assert isinstance(tl, np.ndarray) and tl.shape == (3,)
+
+
+def test_flow_noise_straight_duct_formula() -> None:
+ # Bies Eq. (8.251), VDI 2081-1: note the -2 constant term.
+ f = np.array([250.0])
+ u, s = 10.0, 0.04
+ res = hvac.flow_noise_straight_duct(f, u, s)
+ expected = (
+ 7.0 + 50.0 * math.log10(u) + 10.0 * math.log10(s) - 2.0
+ - 26.0 * math.log10(1.14 + 0.02 * 250.0 / u)
+ )
+ assert res.values[0] == pytest.approx(expected)
+ # Pinned numeric anchor so the constant cannot silently drift.
+ assert res.values[0] == pytest.approx(35.4347, abs=1e-3)
+
+
+def test_flow_noise_scales_with_velocity() -> None:
+ # Regenerated noise rises steeply with flow speed (dominant 50 log10 U term,
+ # tempered by the frequency term which also depends on U).
+ f = np.array([250.0])
+ levels = [hvac.flow_noise_straight_duct(f, u, 0.04).values[0]
+ for u in (5.0, 10.0, 15.0)]
+ assert levels[0] < levels[1] < levels[2]
+ assert levels[2] - levels[0] > 20.0
+
+
+def test_flow_noise_bend_formula() -> None:
+ f = np.array([500.0])
+ u, s, h, rho = 12.0, 0.04, 0.2, 1.206
+ res = hvac.flow_noise_bend(f, u, s, h, density=rho)
+ lws = 30.0 * math.log10(u) + 10.0 * math.log10(s) + 10.0 * math.log10(rho) + 117.0
+ ns = 500.0 * h / u
+ expected = lws - 10.0 * math.log10(1.0 + 0.165 * ns**2) + 30.0 * math.log10(u) - 103.0
+ assert res.values[0] == pytest.approx(expected)
+
+
+def test_plot_and_validation() -> None:
+ import matplotlib
+
+ matplotlib.use("Agg")
+ hvac.end_reflection_loss([125.0, 250.0], 0.3).plot()
+ hvac.flow_noise_straight_duct([250.0, 500.0], 10.0, 0.04).plot()
+ with pytest.raises(ValueError):
+ hvac.end_reflection_loss([125.0], 0.3, termination="bad")
+ with pytest.raises(ValueError):
+ hvac.plenum_attenuation(0.1, 1.0, 20.0, 1.0)
diff --git a/tests/noise_control/test_silencers.py b/tests/noise_control/test_silencers.py
new file mode 100644
index 00000000..e98b5918
--- /dev/null
+++ b/tests/noise_control/test_silencers.py
@@ -0,0 +1,161 @@
+# Copyright (c) 2026. Jose M. Requena-Plens
+"""Tests for reactive silencers (Bies §8.8-8.9, four-pole method).
+
+Oracles: the simple-expansion-chamber closed form ``TL = 10 log10[1 + (1/4)
+(m - 1/m)^2 sin^2(kL)]`` (Bies Eq. (8.111)) with its published peak values
+(m = 2 -> 1.94, 4 -> 6.55, 8 -> 12.18, 16 -> 18.10 dB), the four-pole TL
+reducing to the side-branch closed form (Bies Eq. (8.73)), the quarter-wave
+tube tuning ``f = c/(4 l_e)`` (Example 8.1: 56.6 Hz), the Helmholtz resonance
+``f_0 = (c/2 pi) sqrt(S/(l_e V))``, matrix reciprocity, and limiting cases.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from phonometry.noise_control import silencers as sl
+
+
+def _k(f: np.ndarray, c: float = 343.0) -> np.ndarray:
+ return 2.0 * np.pi * f / c
+
+
+def test_expansion_chamber_matches_closed_form() -> None:
+ f = np.linspace(20.0, 2000.0, 2000)
+ for m in (2.0, 4.0, 8.0, 16.0):
+ s_duct, length = 0.01, 0.3
+ res = sl.expansion_chamber(f, length, m * s_duct, s_duct)
+ cf = 10.0 * np.log10(1.0 + 0.25 * (m - 1.0 / m) ** 2 * np.sin(_k(f) * length) ** 2)
+ assert np.allclose(res.transmission_loss, cf, atol=1e-9)
+
+
+@pytest.mark.parametrize(
+ ("m", "peak"), [(2.0, 1.94), (4.0, 6.55), (8.0, 12.18), (16.0, 18.10)]
+)
+def test_expansion_chamber_peaks(m: float, peak: float) -> None:
+ f = np.linspace(20.0, 2000.0, 6000)
+ res = sl.expansion_chamber(f, 0.3, m * 0.01, 0.01)
+ assert res.transmission_loss.max() == pytest.approx(peak, abs=0.02)
+
+
+def test_expansion_chamber_troughs_zero() -> None:
+ # TL = 0 at kL = n pi (chamber transparent, no dissipation).
+ c, length, s_duct = 343.0, 0.3, 0.01
+ f_trough = c / (2.0 * length) # kL = pi
+ res = sl.expansion_chamber(np.array([f_trough]), length, 0.04, s_duct)
+ assert res.transmission_loss[0] == pytest.approx(0.0, abs=1e-9)
+
+
+def test_side_branch_reduces_to_closed_form() -> None:
+ # A shunt branch Z_b between equal-area ducts: TL = 20 log10|1 + rho c/(2 Sd Zb)|.
+ f = np.linspace(50.0, 500.0, 200)
+ zb = sl.quarter_wave_impedance(f, 0.5, 0.002)
+ t = sl.shunt_matrix(zb)
+ tl = sl.transmission_loss(t, inlet_area=0.01, outlet_area=0.01)
+ closed = 20.0 * np.log10(np.abs(1.0 + 1.206 * 343.0 / (2.0 * 0.01 * zb)))
+ assert np.allclose(tl, closed, atol=1e-9)
+
+
+def test_quarter_wave_tuning_frequency() -> None:
+ # Bies Example 8.1: l_e = 1.516 m, c = 343.24 -> first peak at 56.6 Hz.
+ # Restrict below the third harmonic (169.8 Hz) so the fundamental is the
+ # only resonance in range (all odd harmonics peak equally for a QWT).
+ f = np.linspace(10.0, 120.0, 8000)
+ area = np.pi * 0.05**2 / 4.0
+ res = sl.quarter_wave_resonator(f, area, 1.516, area, speed_of_sound=343.24)
+ peak_f = f[np.argmax(res.transmission_loss)]
+ assert peak_f == pytest.approx(56.6, abs=0.2)
+ assert res.resonances[0] == pytest.approx(56.6, abs=0.1)
+
+
+def test_helmholtz_resonance_frequency() -> None:
+ c, s_neck, le, vol = 343.0, 1e-4, 0.02, 1e-3
+ res = sl.helmholtz_resonator(np.linspace(50.0, 400.0, 50), 0.01, s_neck, le, vol)
+ f0 = c / (2.0 * np.pi) * np.sqrt(s_neck / (le * vol))
+ assert res.resonances[0] == pytest.approx(f0)
+
+
+def test_helmholtz_peak_at_resonance() -> None:
+ # Lossless resonator: TL peaks sharply at f_0.
+ c, s_neck, le, vol = 343.0, 1e-4, 0.02, 1e-3
+ f0 = c / (2.0 * np.pi) * np.sqrt(s_neck / (le * vol))
+ f = np.linspace(0.6 * f0, 1.4 * f0, 4000)
+ res = sl.helmholtz_resonator(f, 0.01, s_neck, le, vol)
+ assert f[np.argmax(res.transmission_loss)] == pytest.approx(f0, rel=0.02)
+
+
+def test_insertion_loss_identity_is_zero() -> None:
+ f = np.linspace(50.0, 500.0, 50)
+ ident = sl.duct_matrix(f, 0.0, 0.01) # zero-length duct = identity
+ il = sl.insertion_loss(ident, source_impedance=400.0, radiation_impedance=50.0)
+ assert np.allclose(il, 0.0, atol=1e-12)
+
+
+def test_insertion_loss_equals_tl_for_anechoic_reference() -> None:
+ # The insertion loss with Z_s = rho c / S_in and Z_r = rho c / S_out equals
+ # the (anechoic) transmission loss; a positive, not negative, quantity.
+ c, rho, s = 343.0, 1.206, 0.01
+ z = rho * c / s
+ f = np.linspace(50.0, 1500.0, 400)
+ t = sl.expansion_chamber(f, 0.3, 0.04, s).transfer_matrix
+ tl = sl.transmission_loss(t, inlet_area=s, outlet_area=s)
+ il = sl.insertion_loss(t, source_impedance=z, radiation_impedance=z)
+ assert np.allclose(il, tl, atol=1e-9)
+ assert il.max() > 3.0 # a real silencer gives a positive insertion loss
+
+
+def test_transfer_matrix_reciprocity() -> None:
+ # Reciprocal passive elements have det(T) = 1.
+ f = np.linspace(50.0, 500.0, 20)
+ t = sl.expansion_chamber(f, 0.3, 0.04, 0.01).transfer_matrix
+ dets = np.linalg.det(t)
+ assert np.allclose(dets, 1.0, atol=1e-9)
+
+
+def test_cascade_order() -> None:
+ f = np.linspace(50.0, 500.0, 10)
+ a = sl.duct_matrix(f, 0.1, 0.01)
+ b = sl.duct_matrix(f, 0.2, 0.01)
+ # Two concatenated ducts equal one duct of the summed length.
+ both = sl.cascade(a, b)
+ whole = sl.duct_matrix(f, 0.3, 0.01)
+ assert np.allclose(both, whole, atol=1e-9)
+
+
+def test_extended_tube_reduces_to_expansion_chamber() -> None:
+ f = np.linspace(20.0, 2000.0, 500)
+ ext = sl.extended_tube_chamber(f, 0.3, 0.04, 0.01)
+ plain = sl.expansion_chamber(f, 0.3, 0.04, 0.01)
+ assert np.allclose(ext.transmission_loss, plain.transmission_loss, atol=1e-9)
+
+
+def test_extended_tube_fills_trough() -> None:
+ # An inlet extension tuned to L/4 raises the TL at the first chamber trough.
+ c, length = 343.0, 0.4
+ f = np.array([c / (2.0 * length)]) # first plain-chamber trough (TL = 0)
+ plain = sl.expansion_chamber(f, length, 0.04, 0.01)
+ tuned = sl.extended_tube_chamber(
+ f, length, 0.04, 0.01, inlet_extension=length / 4.0
+ )
+ assert plain.transmission_loss[0] == pytest.approx(0.0, abs=1e-9)
+ assert tuned.transmission_loss[0] > 3.0
+
+
+def test_insertion_loss_present_when_impedances_given() -> None:
+ f = np.linspace(50.0, 500.0, 20)
+ res = sl.expansion_chamber(
+ f, 0.3, 0.04, 0.01, source_impedance=4e4, radiation_impedance=5e3
+ )
+ assert res.insertion_loss is not None
+ plain = sl.expansion_chamber(f, 0.3, 0.04, 0.01)
+ assert plain.insertion_loss is None
+
+
+def test_validation() -> None:
+ with pytest.raises(ValueError):
+ sl.expansion_chamber([0.0], 0.3, 0.04, 0.01)
+ with pytest.raises(ValueError):
+ sl.expansion_chamber([100.0], 0.3, -0.04, 0.01)
+ with pytest.raises(ValueError, match="must exceed"):
+ sl.extended_tube_chamber([100.0], 0.3, 0.01, 0.02)
diff --git a/tests/test_package_architecture.py b/tests/test_package_architecture.py
index 8a9439e9..299b0a93 100644
--- a/tests/test_package_architecture.py
+++ b/tests/test_package_architecture.py
@@ -35,6 +35,9 @@
("building", "vibration"),
# double-wall cavity fill uses the porous equivalent-fluid model
("building", "materials"),
+ # HVAC plenum and machine enclosures reuse the room constant
+ # R = S*alpha/(1 - alpha) of the steady-state room field
+ ("noise_control", "room"),
}
