[FEATURE] add pytta.analytic with Green's functions, plane waves and a sphere

docs/pyttadoc/analytic.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+Analytic
+========
+
+.. automodule:: pytta.analytic
+   :members:

docs/pyttadoc/index.rst

@@ -57,11 +57,15 @@ Modules:
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.analytic:
+        analytic sound fields of point sources, plane waves and a sphere;
+
 .. toctree::
    classes/index
    utils/index
    apps/index
    functions
    generate
    rooms
+   analytic
    properties

pytta/__init__.py

@@ -72,6 +72,9 @@
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.analytic:
+        analytic sound fields of point sources, plane waves and a sphere;
+
 
 You can find usage examples in the toolbox's examples folder.
 
@@ -108,6 +111,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import analytic
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +134,7 @@
            'generate',
            'utils',
            'iso3741',
+           'analytic',
            'rooms',
 
            # Functions

pytta/analytic.py

@@ -0,0 +1,134 @@
+"""
+PyTTa analytic:
+---------------
+
+Analytic sound fields, as the Monopole tools of the ITA-Toolbox
+(applications/Analytic/Monopole): the free-field Green's function, a point
+source over an impedance plane, a plane wave, and a plane wave scattered by
+a sphere of given surface admittance.
+
+Available functions:
+
+    >>> pytta.analytic.green(f, source, points)
+    >>> pytta.analytic.green_impedance_plane(f, source, points, Z)
+    >>> pytta.analytic.plane_wave(f, direction, points)
+    >>> p_scat, p_inc = pytta.analytic.plane_wave_on_sphere(f, radius, points)
+
+Time dependence exp(j w t), so outgoing waves go as exp(-j k r). Points are
+(n, 3) cartesian coordinates in m; results are (frequencies, points).
+
+References
+----------
+
+Mechel, F. P. (2008). Formulas of Acoustics (2nd ed.). Springer.
+Section on the scattering of a plane wave by a sphere.
+
+Morse, P. M., & Ingard, K. U. (1968). Theoretical Acoustics. McGraw-Hill.
+"""
+
+import numpy as np
+from scipy.special import eval_legendre, spherical_jn, spherical_yn
+
+__all__ = ['green', 'green_impedance_plane', 'plane_wave', 'plane_wave_on_sphere']
+
+_C, _RHO = 343.2, 1.204
+
+
+def _grid(f, points):
+    f = np.atleast_1d(np.asarray(f, dtype=float))
+    P = np.atleast_2d(np.asarray(points, dtype=float))
+    return f, P, 2 * np.pi * f[:, None]
+
+
+def green(f, source, points, c=_C):
+    """
+    Free-field Green's function exp(-j k r) / (4 pi r) from source to each
+    point, as ita_analytic_green_function without a plane.
+    """
+    f, P, w = _grid(f, points)
+    r = np.linalg.norm(P - np.asarray(source, dtype=float), axis=1)[None, :]
+    return np.exp(-1j * w / c * r) / (4 * np.pi * r)
+
+
+def green_impedance_plane(f, source, points, Z, c=_C, rho=_RHO):
+    """
+    Point source over the plane z = 0 of surface impedance Z (Pa s/m, a
+    number or one per frequency), as ita_analytic_green_function with an
+    impedance plane: the free field plus the image source at (x0, y0, -z0)
+    weighted by the plane-wave reflection factor at the specular angle,
+    R = (Z cos(theta) - rho c) / (Z cos(theta) + rho c).
+
+    The image with the plane-wave factor is the usual approximation; it is
+    exact for a rigid (Z = inf) or a pressure-release (Z = 0) plane, and
+    leaves out the ground wave of the exact solution.
+    """
+    f, P, w = _grid(f, points)
+    s = np.asarray(source, dtype=float)
+    image = s * np.array([1, 1, -1])
+    r2 = np.linalg.norm(P - image, axis=1)[None, :]
+    cos_t = (P[:, 2] + s[2])[None, :] / r2
+    Z = np.broadcast_to(np.asarray(Z, dtype=complex).reshape(-1, 1) if np.ndim(Z) else Z,
+                        (f.size, 1))
+    with np.errstate(invalid='ignore', divide='ignore'):
+        R = np.where(np.isinf(Z), 1, (Z * cos_t - rho * c) / (Z * cos_t + rho * c))
+    return green(f, s, P, c) + R * np.exp(-1j * w / c * r2) / (4 * np.pi * r2)
+
+
+def plane_wave(f, direction, points, origin=(0, 0, 0), c=_C):
+    """Plane wave of unit amplitude, exp(-j k u . (r - origin)), u the direction."""
+    f, P, w = _grid(f, points)
+    u = np.asarray(direction, dtype=float)
+    u = u / np.linalg.norm(u)
+    return np.exp(-1j * w / c * ((P - np.asarray(origin, dtype=float)) @ u)[None, :])
+
+
+def plane_wave_on_sphere(f, radius, points, admittance=0., c=_C, max_order=80):
+    """
+    Plane wave exp(-j k z) scattered by a sphere of radius a at the origin,
+    as ita_analytic_plane_wave_on_sphere (series of Mechel, Formulas of
+    Acoustics): with the normalized surface admittance beta (rho c times the
+    admittance, 0 for a rigid sphere), the boundary condition
+    dp/dr = j k beta p at r = a gives
+
+        p_scat = sum_m (2m + 1) (-j)^m D_m P_m(cos theta) h_m(k r),
+        D_m = -(j_m'(ka) - j beta j_m(ka)) / (h_m'(ka) - j beta h_m(ka)),
+
+    with h_m = j_m - j y_m (outgoing for exp(j w t)).
+
+    Parameters
+    ----------
+    f : array_like
+        Frequencies (Hz).
+    radius : float
+    points : array_like
+        (n, 3) points on or outside the sphere.
+    admittance : complex or array_like
+        beta, a number or one per frequency. Default 0.
+    max_order : int
+        Highest order of the series. Default 80.
+
+    Returns
+    -------
+    (numpy.ndarray, numpy.ndarray)
+        Scattered and incident pressure, (frequencies, points).
+    """
+    f, P, w = _grid(f, points)
+    r = np.linalg.norm(P, axis=1)
+    if np.any(r < radius * (1 - 1e-12)):
+        raise ValueError('points must lie on or outside the sphere.')
+    cos_t = np.where(r > 0, P[:, 2] / np.where(r > 0, r, 1), 1)
+    k = w[:, 0] / c
+    beta = np.broadcast_to(np.asarray(admittance, dtype=complex), k.shape)[:, None]
+    m = np.arange(max_order + 1)[None, :]
+    ka = (k * radius)[:, None]
+    j, jd = spherical_jn(m, ka), spherical_jn(m, ka, derivative=True)
+    y, yd = spherical_yn(m, ka), spherical_yn(m, ka, derivative=True)
+    h, hd = j - 1j * y, jd - 1j * yd
+    D = -(jd - 1j * beta * j) / (hd - 1j * beta * h)          # (nf, orders)
+    coef = (2 * m + 1) * (-1j) ** m                            # (1, orders)
+    leg = np.array([eval_legendre(n, cos_t) for n in range(max_order + 1)])   # (orders, np)
+    kr = k[:, None, None] * r[None, None, :]                   # (nf, 1, np)
+    mm = np.arange(max_order + 1)[None, :, None]
+    hr = spherical_jn(mm, kr) - 1j * spherical_yn(mm, kr)      # (nf, orders, np)
+    p_scat = np.einsum('fm,mp,fmp->fp', D * coef, leg, hr)
+    return p_scat, plane_wave(f, (0, 0, 1), P, c=c)

tests/test_analytic.py

@@ -0,0 +1,76 @@
+"""
+Tests for pytta.analytic, after the Monopole tools of the ITA-Toolbox. Each
+test pins one rule of the field: its form, a boundary condition or a limit.
+"""
+import unittest
+
+import numpy as np
+
+from pytta import analytic as A
+
+C, RHO = 343.2, 1.204
+
+
+class TestAnalytic(unittest.TestCase):
+
+    def test_free_field_green_and_plane_wave(self):
+        # |G| = 1 / (4 pi r) with phase -k r, the same both ways; a plane
+        # wave has unit magnitude and phase -k u.r.
+        f = np.array([100., 1000.])
+        src, pts = (0.3, -0.2, 1.0), np.array([[1.0, 2.0, 0.5], [-2.0, 0.5, 3.0]])
+        G = A.green(f, src, pts)
+        r = np.linalg.norm(pts - src, axis=1)
+        k = 2 * np.pi * f[:, None] / C
+        e = np.abs(G - np.exp(-1j * k * r) / (4 * np.pi * r)).max()
+        back = A.green(f, pts[0], [src])[:, 0]
+        pw = A.plane_wave(f, (1, 1, 0), pts)
+        e_pw = np.abs(pw - np.exp(-1j * k * (pts @ np.array([1, 1, 0]) / np.sqrt(2)))).max()
+        print(f'\n  Green {e:.1e}, reciprocity {np.abs(back - G[:, 0]).max():.1e}, plane wave {e_pw:.1e}')
+        self.assertLess(e, 1e-15)
+        self.assertLess(np.abs(back - G[:, 0]).max(), 1e-15)
+        self.assertLess(e_pw, 1e-14)
+
+    def test_impedance_plane_limits(self):
+        # On a rigid plane the image doubles the free field; on a
+        # pressure-release plane the pressure is zero; straight above the
+        # source the image carries R = (Z - rho c) / (Z + rho c).
+        f = np.array([500.])
+        src = (0.0, 0.0, 1.0)
+        on = np.array([[2.0, 1.0, 0.0]])
+        rigid = A.green_impedance_plane(f, src, on, np.inf)
+        soft = A.green_impedance_plane(f, src, on, 0)
+        Z = 3 * RHO * C
+        above = np.array([[0.0, 0.0, 2.0]])
+        p = A.green_impedance_plane(f, src, above, Z)
+        R = (Z - RHO * C) / (Z + RHO * C)
+        expected = A.green(f, src, above) + R * A.green(f, (0, 0, -1.0), above)
+        print(f'\n  rigid / free {abs(rigid[0, 0] / A.green(f, src, on)[0, 0]):.6f}, '
+              f'soft {abs(soft[0, 0]):.1e}, normal incidence {abs(p - expected).max():.1e}')
+        self.assertAlmostEqual(abs(rigid[0, 0] / A.green(f, src, on)[0, 0]), 2, places=12)
+        self.assertLess(abs(soft[0, 0]), 1e-15)
+        self.assertLess(abs(p - expected).max(), 1e-15)
+
+    def test_sphere_boundary_condition(self):
+        # The total field meets dp/dr = j k beta p on the sphere: zero
+        # normal derivative for a rigid sphere, and the admittance condition
+        # for beta = 0.5 + 0.2j, checked by a one-sided second-order
+        # difference in r (points inside the sphere are rejected).
+        a, f = 0.1, np.array([300., 1500.])
+        th = np.linspace(0.1, np.pi - 0.1, 9)
+        u = np.column_stack([np.sin(th), np.zeros_like(th), np.cos(th)])
+        k = 2 * np.pi * f[:, None] / C
+        h = 1e-6
+        for beta in (0, 0.5 + 0.2j):
+            tot = [sum(A.plane_wave_on_sphere(f, a, (a + d) * u, beta, max_order=40)) for d in (0, h, 2 * h)]
+            dpdr = (-3 * tot[0] + 4 * tot[1] - tot[2]) / (2 * h)
+            e = np.abs(dpdr - 1j * k * beta * tot[0]).max() / np.abs(k * tot[0]).max()
+            print(f'\n  beta = {beta}: boundary condition met within {e:.1e}', end='')
+            self.assertLess(e, 1e-6)
+
+    def test_sphere_rejects_points_inside(self):
+        with self.assertRaises(ValueError):
+            A.plane_wave_on_sphere([100.], 0.1, [[0, 0, 0.05]])
+
+
+if __name__ == '__main__':
+    unittest.main()