From f4ccf69ba2c8405a6081fd36c27b8c11df9d51ee Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Thu, 21 Nov 2019 17:43:29 +0100 Subject: [PATCH 1/7] add point source with far-field directivity of circular piston --- sfs/fd/source.py | 62 ++++++++++++++++++++++++++++++++++++++++++++++++ sfs/util.py | 24 ++++++++++++++++++- 2 files changed, 85 insertions(+), 1 deletion(-) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 6a0b8602..604af72a 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -224,6 +224,68 @@ def point_dipole(omega, x0, n0, grid, *, c=None): return 1 / (4 * _np.pi) * (1j * k + 1 / r) * _np.inner(offset, n0) / \ _np.power(r, 2) * _np.exp(-1j * k * r) +def point_circular_piston(omega, x0, n0, grid, R, *, c=None): + r"""Point source with circular piston far-field directivity. + + Parameters + ---------- + omega : float + Frequency of source. + x0 : (3,) array_like + Position of source. + n0 : (3,) array_like + Normal vector (direction) of the circular piston. + grid : triple of array_like + The grid that is used for the sound field calculations. + See `sfs.util.xyz_grid()`. + R : float + Radius of circular piston. + c : float, optional + Speed of sound. + Returns + ------- + numpy.ndarray + Sound pressure at positions given by *grid*. + + Notes + ----- + .. math:: + + G(\x-\x_0,\w) = \frac{1}{2\pi} + \frac{J_1(\wc R \sin(\Theta))}{\wc R \sin(\Theta)} + \frac{\e{-\i\wc|\x-\x_0|}}{|\x-\x_0|} + with + + .. math:: + + \cos(\Theta) = + \frac{\langle \mathbf{x} - \mathbf{x}_0 | \mathbf{n}_0 \rangle} + {|\mathbf{x} - \mathbf{x}_0|} + + Examples + -------- + .. plot:: + :context: close-figs + + n0 = 0, -1, 0 + R = 1.0 + p = sfs.fd.source.point_circular_piston(omega, x0, n0, grid, R) + sfs.plot2d.amplitude(p * normalization_point, grid) + plt.title("Circular Piston Radiator at {} m".format(x0)) + + """ + k = _util.wavenumber(omega, c) + x0 = _util.asarray_1d(x0) + n0 = _util.asarray_1d(n0) + grid = _util.as_xyz_components(grid) + + offset = grid - x0 + r = _np.linalg.norm(offset) + cosphi = _np.inner(offset, n0) / r + sinphi = _np.sqrt(1 - _np.power(cosphi, 2)) + + return 1 / (2 * _np.pi) * _np.exp(-1j * k * r) * _util.jinc(k * R * sinphi) + def point_modal(omega, x0, grid, L, *, N=None, deltan=0, c=None): """Point source in a rectangular room using a modal room model. diff --git a/sfs/util.py b/sfs/util.py index b054346a..acd6ae6e 100644 --- a/sfs/util.py +++ b/sfs/util.py @@ -7,7 +7,8 @@ import collections import numpy as np from numpy.core.umath_tests import inner1d -from scipy.special import spherical_jn, spherical_yn +from numpy.core.numeric import where +from scipy.special import spherical_jn, spherical_yn, j1 from . import default @@ -555,6 +556,27 @@ def spherical_hn2(n, z): return spherical_jn(n, z) - 1j * spherical_yn(n, z) +def jinc(x): + r"""Circular counterpart of Sinc function a.k.a. Jinc function + + .. math:: + + \mathrm{Jinc}(x) = \frac{J_1(x)}{x} + + where :math:`J_1(\cdot)` is the Bessel function of first order, + and :math:`x` its real-value argument. + + Parameters + ---------- + x : array_like + Argument of the Jinc function. + + """ + + x = np.asanyarray(x) + y = where(x == 0, 1.0e-20, x) + return j1(y)/y + def source_selection_plane(n0, n): """Secondary source selection for a plane wave. From 43753b28f7a620f76adf6bee85f4b3ca5d376e60 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 22 Nov 2019 14:50:15 +0100 Subject: [PATCH 2/7] fix wrong distance attenuation --- sfs/fd/source.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 604af72a..6a575d7d 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -284,7 +284,8 @@ def point_circular_piston(omega, x0, n0, grid, R, *, c=None): cosphi = _np.inner(offset, n0) / r sinphi = _np.sqrt(1 - _np.power(cosphi, 2)) - return 1 / (2 * _np.pi) * _np.exp(-1j * k * r) * _util.jinc(k * R * sinphi) + return 1 / (2 * _np.pi) * _np.exp(-1j * k * r) / r \ + _util.jinc(k * R * sinphi) def point_modal(omega, x0, grid, L, *, N=None, deltan=0, c=None): From 8d57e7e5d6211dc6ec48e8a33c5241aca6eb375f Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 22 Nov 2019 15:29:21 +0100 Subject: [PATCH 3/7] fix wrong syntax --- sfs/fd/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 6a575d7d..0ed0fdbf 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -284,7 +284,7 @@ def point_circular_piston(omega, x0, n0, grid, R, *, c=None): cosphi = _np.inner(offset, n0) / r sinphi = _np.sqrt(1 - _np.power(cosphi, 2)) - return 1 / (2 * _np.pi) * _np.exp(-1j * k * r) / r \ + return 1 / (2 * _np.pi) * _np.exp(-1j * k * r) / r * \ _util.jinc(k * R * sinphi) From 50d31dc61eeddbd19accbaf97ca5324f78564193 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Thu, 19 Dec 2019 11:08:52 +0100 Subject: [PATCH 4/7] fix issues with codestyle --- sfs/fd/source.py | 5 +++-- sfs/util.py | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 0ed0fdbf..898f3858 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -224,6 +224,7 @@ def point_dipole(omega, x0, n0, grid, *, c=None): return 1 / (4 * _np.pi) * (1j * k + 1 / r) * _np.inner(offset, n0) / \ _np.power(r, 2) * _np.exp(-1j * k * r) + def point_circular_piston(omega, x0, n0, grid, R, *, c=None): r"""Point source with circular piston far-field directivity. @@ -251,10 +252,10 @@ def point_circular_piston(omega, x0, n0, grid, R, *, c=None): ----- .. math:: - G(\x-\x_0,\w) = \frac{1}{2\pi} + G(\x-\x_0,\w) = \frac{1}{2\pi} \frac{J_1(\wc R \sin(\Theta))}{\wc R \sin(\Theta)} \frac{\e{-\i\wc|\x-\x_0|}}{|\x-\x_0|} - with + with .. math:: diff --git a/sfs/util.py b/sfs/util.py index acd6ae6e..e0e663e7 100644 --- a/sfs/util.py +++ b/sfs/util.py @@ -577,6 +577,7 @@ def jinc(x): y = where(x == 0, 1.0e-20, x) return j1(y)/y + def source_selection_plane(n0, n): """Secondary source selection for a plane wave. From 19c4c4637f2a4c0b27c75266028a5dc39509ac35 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Thu, 19 Dec 2019 11:11:17 +0100 Subject: [PATCH 5/7] fix missing blank line in docstring --- sfs/fd/source.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 898f3858..4b395f5e 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -243,6 +243,7 @@ def point_circular_piston(omega, x0, n0, grid, R, *, c=None): Radius of circular piston. c : float, optional Speed of sound. + Returns ------- numpy.ndarray From 4e7503c13aca513d1de0e8854ac362b0bfcdaada Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 20 Dec 2019 12:32:21 +0100 Subject: [PATCH 6/7] fix docstring in jinc function --- sfs/util.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sfs/util.py b/sfs/util.py index e0e663e7..2f9933cb 100644 --- a/sfs/util.py +++ b/sfs/util.py @@ -557,14 +557,15 @@ def spherical_hn2(n, z): def jinc(x): - r"""Circular counterpart of Sinc function a.k.a. Jinc function + r"""Jinc function (circular counterpart of Sinc function). .. math:: \mathrm{Jinc}(x) = \frac{J_1(x)}{x} where :math:`J_1(\cdot)` is the Bessel function of first order, - and :math:`x` its real-value argument. + and :math:`x` its real-value argument. For reference implementation, see + https://docs.scipy.org/doc/numpy/reference/generated/numpy.sinc.html. Parameters ---------- @@ -572,7 +573,6 @@ def jinc(x): Argument of the Jinc function. """ - x = np.asanyarray(x) y = where(x == 0, 1.0e-20, x) return j1(y)/y From d9dd70f4dc5d0b36bc336459dbd59ba56f30f137 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 20 Dec 2019 13:24:14 +0100 Subject: [PATCH 7/7] explain usage of cosine --- sfs/fd/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sfs/fd/source.py b/sfs/fd/source.py index 4b395f5e..c8092411 100644 --- a/sfs/fd/source.py +++ b/sfs/fd/source.py @@ -256,7 +256,7 @@ def point_circular_piston(omega, x0, n0, grid, R, *, c=None): G(\x-\x_0,\w) = \frac{1}{2\pi} \frac{J_1(\wc R \sin(\Theta))}{\wc R \sin(\Theta)} \frac{\e{-\i\wc|\x-\x_0|}}{|\x-\x_0|} - with + with :math:`\Theta` more conveniently defined by its cosine .. math::