2.5D WFS Point Source Unified / Delft Version

sfs/mono/drivingfunction.py

@@ -156,7 +156,7 @@ def _wfs_point(omega, x0, n0, xs, c=None):
 
 
 def wfs_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
-    r"""Driving function for 2.5-dimensional WFS for a virtual point source.
+    r"""Driving function for 2.5-dimensional WFS of a virtual point source.
 
     Parameters
     ----------
@@ -169,9 +169,11 @@ def wfs_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
     xs : (3,) array_like
         Position of virtual point source.
     xref : (3,) array_like, optional
-        Reference point for synthesized sound field.
+        Reference point xref or contour xref(x0) for amplitude correct synthesis.
     c : float, optional
-        Speed of sound.
+        Speed of sound in m/s.
+    omalias : float, optional
+        Cut angular frequency for prefilter.
 
     Returns
     -------
@@ -186,38 +188,50 @@ def wfs_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
 
     Notes
     -----
+    Eq. (3.10), (3.11) Start, E.W. (1997): "Direct sound enhancement by wave
+    field synthesis", doctoral thesis, TU Delft
+
+    Eq. (2.137) from Schultz, F. (2016): https://doi.org/10.18453/rosdok_id00001765
+
     .. math::
 
-        D(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|}
-            \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}}
-            {|\x_0-\x_\text{s}|^\frac{3}{2}}
-            \e{-\i\wc |\x_0-\x_\text{s}|}
+        D(\x_0,\w) = \sqrt{8 \pi \, \i\wc}
+            \sqrt{\frac{|\x_\text{ref}-\x_0| \cdot
+            |\x_0-\x_\text{s}|}{|\x_\text{ref}-\x_0| + |\x_0-\x_\text{s}|}}
+            \scalarprod{\frac{\x_0-\x_\text{s}}{|\x_0-\x_\text{s}|}}{\n_0}
+            \frac{\e{-\i\wc |\x_0-\x_\text{s}|}}{4\pi\,|\x_0-\x_\text{s}|}
 
     Examples
     --------
     .. plot::
         :context: close-figs
 
-        d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point(
-            omega, array.x, array.n, xs)
-        plot(d, selection, secondary_source)
+        d, selection, secondary_source = (
+            sfs.mono.drivingfunction.wfs_25d_point_Unified_WIP(
+                omega, array.x, array.n, xs))
+        plot(4*np.pi*np.linalg.norm(xs) * d, selection, secondary_source)
 
     """
     x0 = util.asarray_of_rows(x0)
     n0 = util.asarray_of_rows(n0)
     xs = util.asarray_1d(xs)
     xref = util.asarray_1d(xref)
     k = util.wavenumber(omega, c)
+
     ds = x0 - xs
-    r = np.linalg.norm(ds, axis=1)
+    dr = xref - x0
+    s = np.linalg.norm(ds, axis=1)
+    r = np.linalg.norm(dr, axis=1)
+
     d = (
         wfs_25d_preeq(omega, omalias, c) *
-        np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) /
-        r ** (3 / 2) * np.exp(-1j * k * r))
+        np.sqrt(8 * np.pi) *
+        np.sqrt((r * s) / (r + s)) *
+        inner1d(n0, ds) / s *
+        np.exp(-1j * k * s) / (4 * np.pi * s))
     selection = util.source_selection_point(n0, x0, xs)
     return d, selection, secondary_source_point(omega, c)
 
-
 wfs_3d_point = _wfs_point