"""Compute the sound field generated by a sound source.
.. plot::
:context: reset
import sfs
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 8, 4.5 # inch
x0 = 1.5, 1, 0
f = 500 # Hz
omega = 2 * np.pi * f
normalization_point = 4 * np.pi
normalization_line = \\
np.sqrt(8 * np.pi * omega / sfs.defs.c) * np.exp(1j * np.pi / 4)
grid = sfs.util.xyz_grid([-2, 3], [-1, 2], 0, spacing=0.02)
# Grid for vector fields:
vgrid = sfs.util.xyz_grid([-2, 3], [-1, 2], 0, spacing=0.1)
"""
import itertools
import numpy as np
from scipy import special
from .. import util
from .. import defs
[docs]def point(omega, x0, n0, grid, c=None):
"""Point source.
Notes
-----
::
1 e^(-j w/c |x-x0|)
G(x-x0, w) = --- -----------------
4pi |x-x0|
Examples
--------
.. plot::
:context: close-figs
p = sfs.mono.source.point(omega, x0, None, grid)
sfs.plot.soundfield(p, grid)
plt.title("Point Source at {} m".format(x0))
Normalization ...
.. plot::
:context: close-figs
sfs.plot.soundfield(p * normalization_point, grid,
colorbar_kwargs=dict(label="p / Pa"))
plt.title("Point Source at {} m (normalized)".format(x0))
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
grid = util.as_xyz_components(grid)
r = np.linalg.norm(grid - x0)
return 1 / (4*np.pi) * np.exp(-1j * k * r) / r
[docs]def point_velocity(omega, x0, n0, grid, c=None):
"""Velocity of a point source.
Returns
-------
XyzComponents
Particle velocity at positions given by `grid`.
See :class:`sfs.util.XyzComponents`.
Examples
--------
The particle velocity can be plotted on top of the sound pressure:
.. plot::
:context: close-figs
v = sfs.mono.source.point_velocity(omega, x0, None, vgrid)
sfs.plot.soundfield(p * normalization_point, grid)
sfs.plot.vectors(v * normalization_point, vgrid)
plt.title("Sound Pressure and Particle Velocity")
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
grid = util.as_xyz_components(grid)
offset = grid - x0
r = np.linalg.norm(offset)
v = point(omega, x0, n0, grid, c=c)
v *= (1+1j*k*r) / (defs.rho0 * defs.c * 1j*k*r)
return util.XyzComponents([v * o / r for o in offset])
[docs]def point_dipole(omega, x0, n0, grid, c=None):
"""Point source with dipole characteristics.
Parameters
----------
omega : float
Frequency of source.
x0 : (3,) array_like
Position of source.
n0 : (3,) array_like
Normal vector (direction) of dipole.
grid : triple of array_like
The grid that is used for the sound field calculations.
See :func:`sfs.util.xyz_grid`.
c : float, optional
Speed of sound.
Returns
-------
numpy.ndarray
Sound pressure at positions given by `grid`.
Notes
-----
::
d 1 / iw 1 \ (x-x0) n0
---- G(x-x0,w) = --- | ----- + ------- | ----------- e^(-i w/c |x-x0|)
d ns 4pi \ c |x-x0| / |x-x0|^2
Examples
--------
.. plot::
:context: close-figs
n0 = 0, 1, 0
p = sfs.mono.source.point_dipole(omega, x0, n0, grid)
sfs.plot.soundfield(p, grid)
plt.title("Dipole Point Source 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)
return 1 / (4*np.pi) * (1j * k + 1 / r) * np.inner(offset, n0) / \
np.power(r, 2) * np.exp(-1j * k * r)
[docs]def point_modal(omega, x0, n0, grid, L, N=None, deltan=0, c=None):
"""Point source in a rectangular room using a modal room model.
Parameters
----------
omega : float
Frequency of source.
x0 : (3,) array_like
Position of source.
n0 : (3,) array_like
Normal vector (direction) of source (only required for
compatibility).
grid : triple of array_like
The grid that is used for the sound field calculations.
See :func:`sfs.util.xyz_grid`.
L : (3,) array_like
Dimensionons of the rectangular room.
N : (3,) array_like or int, optional
Combination of modal orders in the three-spatial dimensions to
calculate the sound field for or maximum order for all
dimensions. If not given, the maximum modal order is
approximately determined and the sound field is computed up to
this maximum order.
deltan : float, optional
Absorption coefficient of the walls.
c : float, optional
Speed of sound.
Returns
-------
numpy.ndarray
Sound pressure at positions given by `grid`.
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
x, y, z = util.as_xyz_components(grid)
if N is None:
# determine maximum modal order per dimension
Nx = int(np.ceil(L[0]/np.pi * k))
Ny = int(np.ceil(L[1]/np.pi * k))
Nz = int(np.ceil(L[2]/np.pi * k))
mm = range(Nx)
nn = range(Ny)
ll = range(Nz)
elif np.isscalar(N):
# compute up to a given order
mm = range(N)
nn = range(N)
ll = range(N)
else:
# compute field for one order combination only
mm = [N[0]]
nn = [N[1]]
ll = [N[2]]
kmp0 = [((kx + 1j * deltan)**2, np.cos(kx * x) * np.cos(kx * x0[0]))
for kx in [m * np.pi / L[0] for m in mm]]
kmp1 = [((ky + 1j * deltan)**2, np.cos(ky * y) * np.cos(ky * x0[1]))
for ky in [n * np.pi / L[1] for n in nn]]
kmp2 = [((kz + 1j * deltan)**2, np.cos(kz * z) * np.cos(kz * x0[2]))
for kz in [l * np.pi / L[2] for l in ll]]
ksquared = k**2
p = 0
for (km0, p0), (km1, p1), (km2, p2) in itertools.product(kmp0, kmp1, kmp2):
km = km0 + km1 + km2
p = p + 8 / (ksquared - km) * p0 * p1 * p2
return p
[docs]def point_modal_velocity(omega, x0, n0, grid, L, N=None, deltan=0, c=None):
"""Velocity of point source in a rectangular room using a modal room model.
Parameters
----------
omega : float
Frequency of source.
x0 : (3,) array_like
Position of source.
n0 : (3,) array_like
Normal vector (direction) of source (only required for
compatibility).
grid : triple of array_like
The grid that is used for the sound field calculations.
See :func:`sfs.util.xyz_grid`.
L : (3,) array_like
Dimensionons of the rectangular room.
N : (3,) array_like or int, optional
Combination of modal orders in the three-spatial dimensions to
calculate the sound field for or maximum order for all
dimensions. If not given, the maximum modal order is
approximately determined and the sound field is computed up to
this maximum order.
deltan : float, optional
Absorption coefficient of the walls.
c : float, optional
Speed of sound.
Returns
-------
XyzComponents
Particle velocity at positions given by `grid`.
See :class:`sfs.util.XyzComponents`.
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
x, y, z = util.as_xyz_components(grid)
if N is None:
# determine maximum modal order per dimension
Nx = int(np.ceil(L[0]/np.pi * k))
Ny = int(np.ceil(L[1]/np.pi * k))
Nz = int(np.ceil(L[2]/np.pi * k))
mm = range(Nx)
nn = range(Ny)
ll = range(Nz)
elif np.isscalar(N):
# compute up to a given order
mm = range(N)
nn = range(N)
ll = range(N)
else:
# compute field for one order combination only
mm = [N[0]]
nn = [N[1]]
ll = [N[2]]
kmp0 = [((kx + 1j * deltan)**2, np.sin(kx * x) * np.cos(kx * x0[0]))
for kx in [m * np.pi / L[0] for m in mm]]
kmp1 = [((ky + 1j * deltan)**2, np.sin(ky * y) * np.cos(ky * x0[1]))
for ky in [n * np.pi / L[1] for n in nn]]
kmp2 = [((kz + 1j * deltan)**2, np.sin(kz * z) * np.cos(kz * x0[2]))
for kz in [l * np.pi / L[2] for l in ll]]
ksquared = k**2
vx = 0+0j
vy = 0+0j
vz = 0+0j
for (km0, p0), (km1, p1), (km2, p2) in itertools.product(kmp0, kmp1, kmp2):
km = km0 + km1 + km2
vx = vx - 8*1j / (ksquared - km) * p0
vy = vy - 8*1j / (ksquared - km) * p1
vz = vz - 8*1j / (ksquared - km) * p2
return util.XyzComponents([vx, vy, vz])
[docs]def line(omega, x0, n0, grid, c=None):
"""Line source parallel to the z-axis.
Note: third component of x0 is ignored.
Notes
-----
::
(2)
G(x-x0, w) = -j/4 H0 (w/c |x-x0|)
Examples
--------
.. plot::
:context: close-figs
p = sfs.mono.source.line(omega, x0, None, grid)
sfs.plot.soundfield(p, grid)
plt.title("Line Source at {} m".format(x0[:2]))
Normalization ...
.. plot::
:context: close-figs
sfs.plot.soundfield(p * normalization_line, grid,
colorbar_kwargs=dict(label="p / Pa"))
plt.title("Line Source at {} m (normalized)".format(x0[:2]))
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)[:2] # ignore z-component
grid = util.as_xyz_components(grid)
r = np.linalg.norm(grid[:2] - x0)
p = -1j/4 * _hankel2_0(k * r)
return _duplicate_zdirection(p, grid)
[docs]def line_velocity(omega, x0, n0, grid, c=None):
"""Velocity of line source parallel to the z-axis.
Returns
-------
XyzComponents
Particle velocity at positions given by `grid`.
See :class:`sfs.util.XyzComponents`.
Examples
--------
The particle velocity can be plotted on top of the sound pressure:
.. plot::
:context: close-figs
v = sfs.mono.source.line_velocity(omega, x0, None, vgrid)
sfs.plot.soundfield(p * normalization_line, grid)
sfs.plot.vectors(v * normalization_line, vgrid)
plt.title("Sound Pressure and Particle Velocity")
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)[:2] # ignore z-component
grid = util.as_xyz_components(grid)
offset = grid[:2] - x0
r = np.linalg.norm(offset)
v = -1/(4*defs.c*defs.rho0) * special.hankel2(1, k * r)
v = [v * o / r for o in offset]
assert v[0].shape == v[1].shape
if len(grid) > 2:
v.append(np.zeros_like(v[0]))
return util.XyzComponents([_duplicate_zdirection(vi, grid) for vi in v])
[docs]def line_dipole(omega, x0, n0, grid, c=None):
"""Line source with dipole characteristics parallel to the z-axis.
Note: third component of x0 is ignored.
Notes
-----
::
(2)
G(x-x0, w) = jk/4 H1 (w/c |x-x0|) cos(phi)
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)[:2] # ignore z-components
n0 = util.asarray_1d(n0)[:2]
grid = util.as_xyz_components(grid)
dx = grid[:2] - x0
r = np.linalg.norm(dx)
p = 1j*k/4 * special.hankel2(1, k * r) * np.inner(dx, n0) / r
return _duplicate_zdirection(p, grid)
[docs]def line_dirichlet_edge(omega, x0, grid, alpha=3/2*np.pi, Nc=None, c=None):
""" Sound field of an line source scattered at an edge with Dirichlet
boundary conditions.
[Michael Möser, Technische Akustik, 2012, Springer, eq.(10.18/19)]
Parameters
----------
omega : float
Angular frequency.
x0 : (3,) array_like
Position of line source.
XyzComponents
A grid that is used for calculations of the sound field.
alpha : float, optional
Outer angle of edge.
Nc : int, optional
Number of elements for series expansion of driving function. Estimated
if not given.
c : float, optional
Speed of sound
Returns
-------
() numpy.ndarray
Complex pressure at grid positions.
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
phi_s = np.arctan2(x0[1], x0[0])
if phi_s < 0:
phi_s = phi_s + 2*np.pi
r_s = np.linalg.norm(x0)
grid = util.XyzComponents(grid)
r = np.linalg.norm(grid[:2])
phi = np.arctan2(grid[1], grid[0])
phi = np.where(phi < 0, phi+2*np.pi, phi)
if Nc is None:
Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi)
epsilon = np.ones(Nc) # weights for series expansion
epsilon[0] = 2
p = np.zeros((grid[0].shape[1], grid[1].shape[0]), dtype=complex)
idxr = (r <= r_s)
idxa = (phi <= alpha)
for m in np.arange(Nc):
nu = m*np.pi/alpha
f = 1/epsilon[m] * np.sin(nu*phi_s) * np.sin(nu*phi)
p[idxr & idxa] = p[idxr & idxa] + f[idxr & idxa] * \
special.jn(nu, k*r[idxr & idxa]) * special.hankel2(nu, k*r_s)
p[~idxr & idxa] = p[~idxr & idxa] + f[~idxr & idxa] * \
special.jn(nu, k*r_s) * special.hankel2(nu, k*r[~idxr & idxa])
p = p * -1j*np.pi/alpha
pl = line(omega, x0, None, grid, c=c)
p[~idxa] = pl[~idxa]
return p
[docs]def plane(omega, x0, n0, grid, c=None):
"""Plane wave.
Notes
-----
::
G(x, w) = e^(-i w/c n x)
Examples
--------
.. plot::
:context: close-figs
direction = 45 # degree
n0 = sfs.util.direction_vector(np.radians(direction))
p = sfs.mono.source.plane(omega, x0, n0, grid)
sfs.plot.soundfield(p, grid, colorbar_kwargs=dict(label="p / Pa"))
plt.title("Plane wave with direction {} degree".format(direction))
"""
k = util.wavenumber(omega, c)
x0 = util.asarray_1d(x0)
n0 = util.asarray_1d(n0)
grid = util.as_xyz_components(grid)
return np.exp(-1j * k * np.inner(grid - x0, n0))
[docs]def plane_velocity(omega, x0, n0, grid, c=None):
"""Velocity of a plane wave.
Notes
-----
::
V(x, w) = 1/(rho c) e^(-i w/c n x) n
Returns
-------
XyzComponents
Particle velocity at positions given by `grid`.
See :class:`sfs.util.XyzComponents`.
Examples
--------
The particle velocity can be plotted on top of the sound pressure:
.. plot::
:context: close-figs
v = sfs.mono.source.plane_velocity(omega, x0, n0, vgrid)
sfs.plot.soundfield(p, grid)
sfs.plot.vectors(v, vgrid)
plt.title("Sound Pressure and Particle Velocity")
"""
v = plane(omega, x0, n0, grid, c=c) / (defs.rho0 * defs.c)
return util.XyzComponents([v * n for n in n0])
def _duplicate_zdirection(p, grid):
"""If necessary, duplicate field in z-direction."""
gridshape = np.broadcast(*grid).shape
if len(gridshape) > 2:
return np.tile(p, [1, 1, gridshape[2]])
else:
return p
def _hankel2_0(x):
"""Wrapper for Hankel function of the second type using fast versions
of the Bessel functions of first/second kind in scipy"""
return special.j0(x)-1j*special.y0(x)