Update wfs.py

- add two examples for xref(x0) reference contour into docstring

Author: fs446

sfs/fd/wfs.py

@@ -217,6 +217,113 @@ 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
+
+        array = sfs.array.circular(N=128, R=R)        
+        f = 343*4  # Hz
+        omega = 2 * np.pi * f  # rad/s
+
+        # example is hard coded for a point source on negative x-axis:
+        xs = -4, 0, 0
+        xref_line = 0  # reference contour is a straight line
+
+        # calc xref(x0):
+        x0_ = array.x.T[np.newaxis, :]
+        xs_ = np.array(xs)[np.newaxis, :, np.newaxis]
+        x0xs = x0_ - xs_
+        x0xs_length = np.linalg.norm(x0xs, axis=1)
+        x0xs_unit = x0xs / x0xs_length
+        n0_ = array.n.T[np.newaxis, :]
+        xref = np.zeros_like(x0_)
+        # cf. https://github.com/spatialaudio/wfs_chapter_hda/blob/master/python/wfs25d_circSSD.py#L91  # noqa
+        # hard coded for a point source on negative x-axis
+        for i in range(array.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_[0, 0, 0] + R * cosbeta
+            xref[0, :, i] = x0_[0, :, i] + tmp
+        xref = np.squeeze(xref).T
+
+        d, selection, secondary_source = sfs.fd.wfs.point_25d(
+            omega, array.x, array.n, xs, xref=xref)
+        normalize_gain = 4 * np.pi * np.linalg.norm(xs)
+        p = sfs.fd.synthesize(d * normalize_gain,
+                              selection,
+                              array,
+                              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(array.x, array.n,
+                                selection * array.a,
+                                size=0.15, ax=ax_amp)
+        sfs.plot2d.loudspeakers(array.x, array.n,
+                                selection * array.a,
+                                size=0.15, 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
+
+        f = 343*4  # Hz
+        omega = 2 * np.pi * f  # rad/s
+
+        # example is hard coded for a point source on negative x-axis:
+        xs = -4, 0, 0
+        xref_dist = 4  # radius for reference contour circle with origin xs
+
+        # calc xref(x0):
+        x0_ = array.x.T[np.newaxis, :]
+        xs_ = np.array(xs)[np.newaxis, :, np.newaxis]
+        x0xs = x0_ - xs_
+        x0xs_length = np.linalg.norm(x0xs, axis=1)
+        x0xs_unit = x0xs / x0xs_length
+        xref = x0_ + (xref_dist - x0xs_length) * x0xs_unit
+        xref = np.squeeze(xref).T
+
+        d, selection, secondary_source = sfs.fd.wfs.point_25d(
+            omega, array.x, array.n, xs, xref=xref)
+        normalize_gain = 4 * np.pi * np.linalg.norm(xs)
+        p = sfs.fd.synthesize(d * normalize_gain,
+                              selection,
+                              array,
+                              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(array.x, array.n,
+                                selection * array.a,
+                                size=0.15, ax=ax_amp)
+        sfs.plot2d.loudspeakers(array.x, array.n,
+                                selection * array.a,
+                                size=0.15, 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)