Source code for sfs.mono.source

"""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  # inch

    x0 = 1.5, 1, 0
    f = 500  # Hz
    omega = 2 * np.pi * f

    grid = sfs.util.xyz_grid([-2, 3], [-1, 2], 0, spacing=0.02)

"""

import itertools
import numpy as np
from scipy import special
from .. import util


[docs]def point(omega, x0, n0, grid, c=None): """Point source. :: 1 e^(-j w/c |x-x0|) G(x-x0, w) = --- ----------------- 4pi |x-x0| Examples -------- .. plot:: :context: close-figs p_point = sfs.mono.source.point(omega, x0, None, grid) sfs.plot.soundfield(p_point, grid) plt.title("Point Source at {} m".format(x0)) Normalization ... multiply by :math:`4\pi` ... .. plot:: :context: close-figs p_point *= 4 * np.pi sfs.plot.soundfield(p_point, grid) plt.title("Point Source at {} m (normalized)".format(x0)) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) grid = util.asarray_of_arrays(grid) r = np.linalg.norm(grid - x0) return np.squeeze(1/(4*np.pi) * np.exp(-1j * k * r) / 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 numpy.ndarray The grid that is used for the sound field calculations. 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.asarray_of_arrays(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 + 1 / (ksquared - km) * p0 * p1 * p2 return np.squeeze(p)
[docs]def line(omega, x0, n0, grid, c=None): """Line source parallel to the z-axis. Note: third component of x0 is ignored. :: (2) G(x-x0, w) = -j/4 H0 (w/c |x-x0|) Examples -------- .. plot:: :context: close-figs p_line = sfs.mono.source.line(omega, x0, None, grid) sfs.plot.soundfield(p_line, grid) plt.title("Line Source at {} m".format(x0[:2])) Normalization ... .. plot:: :context: close-figs p_line *= np.exp(-1j*7*np.pi/4) / np.sqrt(1/(8*np.pi*omega/sfs.defs.c)) sfs.plot.soundfield(p_line, grid) plt.title("Line Source at {} m (normalized)".format(x0[:2])) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x0 = x0[:2] # ignore z-component grid = util.asarray_of_arrays(grid) r = np.linalg.norm(grid[:2] - x0) p = -1j/4 * special.hankel2(0, k * r) p = _duplicate_zdirection(p, grid) return np.squeeze(p)
[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. :: (2) G(x-x0, w) = jk/4 H1 (w/c |x-x0|) cos(phi) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x0 = x0[:2] # ignore z-component n0 = n0[:2] grid = util.asarray_of_arrays(grid) dx = grid[:2] - x0 r = np.linalg.norm(dx) p = 1j*k/4 * special.hankel2(1, k * r) * np.inner(dx, n0) / r p = _duplicate_zdirection(p, grid) return np.squeeze(p)
[docs]def plane(omega, x0, n0, grid, c=None): """Plane wave. :: G(x, w) = e^(-i w/c n x) Example ------- .. plot:: :context: close-figs direction = 45 # degree n0 = sfs.util.direction_vector(np.radians(direction)) p_plane = sfs.mono.source.plane(omega, x0, n0, grid) sfs.plot.soundfield(p_plane, grid); 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.asarray_of_arrays(grid) return np.squeeze(np.exp(-1j * k * np.inner(grid - x0, 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