add 2.5D WFS reference contour examples

sfs/fd/wfs.py

@@ -24,6 +24,7 @@
     grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02)
 
     array = sfs.array.circular(N=32, R=R)
+    array2 = sfs.array.circular(N=64, R=R)
 
     def plot(d, selection, secondary_source):
         p = sfs.fd.synthesize(d, selection, array, secondary_source, grid=grid)
@@ -217,6 +218,112 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
         normalize_gain = 4 * np.pi * np.linalg.norm(xs)
         plot(normalize_gain * d, selection, secondary_source)
 
+    .. plot::
+        :context: close-figs
+
+        xs = -4, 0, 0  # point source on negative x-axis
+        xref_line = 0  # reference contour is a straight line
+
+        lmb = 1 / 4  # m
+        f_tmp = sfs.default.c / lmb  # Hz
+        omega_tmp = 2 * np.pi * f_tmp  # rad/s
+
+        # calc reference contour xref(x0):
+        x0_tmp = array2.x.T[np.newaxis, :]
+        xs_tmp = np.array(xs)[np.newaxis, :, np.newaxis]
+        x0xs = x0_tmp - xs_tmp
+        x0xs_length = np.linalg.norm(x0xs, axis=1)
+        x0xs_unit = x0xs / x0xs_length
+        n0_ = array2.n.T[np.newaxis, :]
+        xref = np.zeros_like(x0_tmp)
+        # cf. https://github.com/spatialaudio/wfs_chapter_hda/blob/master/python/wfs25d_circSSD.py#L91  # noqa
+        # code is valid for a virtual point source on negative x-axis:
+        for i in range(array2.x.shape[0]):
+            cosbeta = np.dot(-n0_[0, :, i], [-1, 0, 0])  # use outward SSD normal
+            tmp = x0xs_unit[0, :, i]
+            tmp *= -x0xs_length[0, i]
+            tmp *= xref_line + R * cosbeta
+            tmp /= xs_tmp[0, 0, 0] + R * cosbeta
+            xref[0, :, i] = x0_tmp[0, :, i] + tmp
+        xref = np.squeeze(xref).T
+
+        d, selection, secondary_source = sfs.fd.wfs.point_25d(
+            omega_tmp, array2.x, array2.n, xs, xref=xref)
+        normalize_gain = 4 * np.pi * np.linalg.norm(xs)
+        p = sfs.fd.synthesize(d * normalize_gain,
+                              selection,
+                              array2,
+                              secondary_source,
+                              grid=grid)
+
+        plt.rcParams['figure.figsize'] = 10, 8
+        ax_amp, ax_lvl = plt.subplot(2, 2, 1), plt.subplot(2, 2, 2)
+        sfs.plot2d.amplitude(p, grid, vmax=2, vmin=-2, ax=ax_amp,
+                             colorbar_kwargs={'label': 'lin'})
+        sfs.plot2d.level(p, grid, vmax=6, vmin=-6,
+                        cmap='seismic', ax=ax_lvl,
+                        colorbar_kwargs={'label': 'dB'})
+        sfs.plot2d.loudspeakers(array2.x, array2.n,
+                                selection * array2.a,
+                                size=0.125, ax=ax_amp)
+        sfs.plot2d.loudspeakers(array2.x, array2.n,
+                                selection * array2.a,
+                                size=0.125, ax=ax_lvl)
+        # plot xref(x0) contour:
+        ax_amp.plot(xref[:, 0][selection],
+                    xref[:, 1][selection], 'C5o', ms=4)
+        ax_lvl.plot(xref[:, 0][selection],
+                    xref[:, 1][selection], 'C5o', ms=4)
+        ax_amp.grid(True), ax_lvl.grid(True)
+
+    .. plot::
+        :context: close-figs
+
+        xs = -4, 0, 0  # point source on negative x-axis
+        xref_dist = 4  # radius for reference contour circle with origin xs
+
+        lmb = 1 / 4  # m
+        f_tmp = sfs.default.c / lmb  # Hz
+        omega_tmp = 2 * np.pi * f_tmp  # rad/s
+
+        # calc reference contour xref(x0):
+        x0_tmp = array2.x.T[np.newaxis, :]
+        xs_tmp = np.array(xs)[np.newaxis, :, np.newaxis]
+        x0xs = x0_tmp - xs_tmp
+        x0xs_length = np.linalg.norm(x0xs, axis=1)
+        x0xs_unit = x0xs / x0xs_length
+        xref = x0_tmp + (xref_dist - x0xs_length) * x0xs_unit
+        xref = np.squeeze(xref).T
+
+        d, selection, secondary_source = sfs.fd.wfs.point_25d(
+            omega_tmp, array2.x, array2.n, xs, xref=xref)
+        normalize_gain = 4 * np.pi * np.linalg.norm(xs)
+        p = sfs.fd.synthesize(d * normalize_gain,
+                              selection,
+                              array2,
+                              secondary_source,
+                              grid=grid)
+
+        plt.rcParams['figure.figsize'] = 10, 8
+        ax_amp, ax_lvl = plt.subplot(2, 2, 1), plt.subplot(2, 2, 2)
+        sfs.plot2d.amplitude(p, grid, vmax=2, vmin=-2, ax=ax_amp,
+                             colorbar_kwargs={'label': 'lin'})
+        sfs.plot2d.level(p, grid, vmax=6, vmin=-6,
+                        cmap='seismic', ax=ax_lvl,
+                        colorbar_kwargs={'label': 'dB'})
+        sfs.plot2d.loudspeakers(array2.x, array2.n,
+                                selection * array2.a,
+                                size=0.125, ax=ax_amp)
+        sfs.plot2d.loudspeakers(array2.x, array2.n,
+                                selection * array2.a,
+                                size=0.125, ax=ax_lvl)
+        # plot xref(x0) contour:
+        ax_amp.plot(xref[:, 0][selection],
+                    xref[:, 1][selection], 'C5o', ms=4)
+        ax_lvl.plot(xref[:, 0][selection],
+                    xref[:, 1][selection], 'C5o', ms=4)
+        ax_amp.grid(True), ax_lvl.grid(True)
+
     """
     x0 = _util.asarray_of_rows(x0)
     n0 = _util.asarray_of_rows(n0)