The sound field¶

Each source's loudness at every cell it reaches, falling off with the shell and shaded by what stands between.

Requires

each cell's sound Power P (W; zero where nothing sounds), its Transmittance of sound t (1, 0..1: the share a cell lets through), the Reach r of a source (cells), the lattice's spacing (m).

Provides

the Intensity at each cell (W/m on a plane: the power over the shell's length; W/m^2 in space).

Law

from each source, I = P / (2 pi d) times the product of the transmittances of the cells strictly between the source and the cell, on the line from one to the other, for every cell within the reach; the sources' intensities add. On a plane the shell is the circle's length, 2 pi d; in space the sphere's area, 4 pi d^2.

Source

Kinsler et al., Fundamentals of Acoustics, 4th ed., 5.11 (the spherical wave's intensity); Kuttruff, Room Acoustics, ch. 1.

Proof

curve, invariant, limit.

No Trait form: waits on the middle: a source's rays are rays.fan's samples, folded with `@ before` as light.point_source's, which draft 5 does not write yet; until then a reference.

sound.field(power, transmittance, reach, spacing)

Each source's loudness at every cell it reaches: the power over the shell, shaded by what stands between (the point-source law).

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import library
from lusine import dmath as D
from lusine.kernels.sound.field import field

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import library from lusine import dmath as D from lusine.kernels.sound.field import field

The reference¶

The kernel's reference, in numpy, over rays.fan's samples, one ray to each cell within the reach, folded with @ before; its Trait form is not written yet.

In [2]:

library.reference(field)

library.reference(field)

Out[2]:

the reference, field.py

def field(power, transmittance, reach, spacing=1.0):
    """THE SOUND FIELD: each source's loudness at every cell it reaches, falling
    off with the shell and shaded by what stands between.

    Requires: each cell's sound `Power` P (W; zero where nothing sounds), its
    `Transmittance` of sound t (1, 0..1: the share a cell lets through), the
    `Reach` r of a source (cells), the lattice's spacing (m).
    Provides: the `Intensity` at each cell (W/m on a plane: the power over the
    shell's length; W/m^2 in space).
    Law: from each source, I = P / (2 pi d) times the product of the
    transmittances of the cells strictly between the source and the cell, on
    the line from one to the other, for every cell within the reach; the
    sources' intensities add. On a plane the shell is the circle's length,
    2 pi d; in space the sphere's area, 4 pi d^2.
    Source: Kinsler et al., Fundamentals of Acoustics, 4th ed., 5.11 (the
    spherical wave's intensity); Kuttruff, Room Acoustics, ch. 1.
    Proof: curve, invariant, limit."""
    P = np.asarray(power, dtype=np.float64)
    t = np.asarray(transmittance, dtype=np.float64)
    out = np.zeros_like(P)
    offsets = within(reach)
    for sy, sx in sources(P):
        for (dy, dx), (y, x), between in _received(P.shape, (sy, sx), t, reach, offsets):
            d = float(np.sqrt(dy * dy + dx * dx)) * float(spacing)
            out[y, x] += P[sy, sx] / (2.0 * D.PI * d) * between
    return out

Curve: the shell's falloff¶

In the open the intensity on a plane falls as 1/d: the power spread over the circle's length, 2 pi d.

In [3]:

n, c = 61, 30
P = np.zeros((n, n)); P[c, c] = 1.0
I = field(P, np.ones((n, n)), 25)
d = np.arange(1, 26)
got = np.array([I[c, c + k] for k in d])
assert np.allclose(got * 2 * D.PI * d, 1.0, rtol=0, atol=1e-12)
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.loglog(d, got, 'o', label='kernel')
ax.loglog(d, 1 / (2 * D.PI * d), label='P / (2 pi d)')
ax.set_xlabel('distance (cells)')
ax.set_ylabel('intensity (W/m)')
ax.legend()
plt.show()

n, c = 61, 30 P = np.zeros((n, n)); P[c, c] = 1.0 I = field(P, np.ones((n, n)), 25) d = np.arange(1, 26) got = np.array([I[c, c + k] for k in d]) assert np.allclose(got * 2 * D.PI * d, 1.0, rtol=0, atol=1e-12) fig, ax = plt.subplots(figsize=(7, 3.2)) ax.loglog(d, got, 'o', label='kernel') ax.loglog(d, 1 / (2 * D.PI * d), label='P / (2 pi d)') ax.set_xlabel('distance (cells)') ax.set_ylabel('intensity (W/m)') ax.legend() plt.show()

Invariant: fields add¶

The law is linear in the power: two sources' field is the sum of each one's, walls and all.

In [4]:

rng = np.random.default_rng(1)
t = rng.uniform(0.3, 1.0, (31, 31))
a = np.zeros((31, 31)); a[10, 10] = 1.0
b = np.zeros((31, 31)); b[20, 22] = 3.0
assert np.allclose(field(a + b, t, 12), field(a, t, 12) + field(b, t, 12), rtol=0, atol=1e-14)
print('superposition holds')

rng = np.random.default_rng(1) t = rng.uniform(0.3, 1.0, (31, 31)) a = np.zeros((31, 31)); a[10, 10] = 1.0 b = np.zeros((31, 31)); b[20, 22] = 3.0 assert np.allclose(field(a + b, t, 12), field(a, t, 12) + field(b, t, 12), rtol=0, atol=1e-14) print('superposition holds')

superposition holds

Limit: a wall that lets nothing through¶

Behind an opaque wall there is silence; a wall letting half through halves the field straight behind it (an oblique ray crosses a thin wall in more cells, its longer path through it, and is shaded more).

In [5]:

n, c = 41, 20
P = np.zeros((n, n)); P[c, c] = 1.0
open_ = field(P, np.ones((n, n)), 15)
for share in (0.0, 0.5):
    w = np.ones((n, n)); w[:, c + 4] = share
    I = field(P, w, 15)
    assert np.allclose(I[c, c + 5:], share * open_[c, c + 5:], rtol=0, atol=1e-15)
    if share == 0.0:
        assert float(I[:, c + 5:].sum()) == 0.0
fig, ax = plt.subplots(1, 2, figsize=(8, 3.6))
ax[0].imshow(np.log10(open_ + 1e-6)); ax[0].set_title('open')
w = np.ones((n, n)); w[5:36, c + 4] = 0.0
ax[1].imshow(np.log10(field(P, w, 15) + 1e-6)); ax[1].set_title('a wall')
plt.show()

n, c = 41, 20 P = np.zeros((n, n)); P[c, c] = 1.0 open_ = field(P, np.ones((n, n)), 15) for share in (0.0, 0.5): w = np.ones((n, n)); w[:, c + 4] = share I = field(P, w, 15) assert np.allclose(I[c, c + 5:], share * open_[c, c + 5:], rtol=0, atol=1e-15) if share == 0.0: assert float(I[:, c + 5:].sum()) == 0.0 fig, ax = plt.subplots(1, 2, figsize=(8, 3.6)) ax[0].imshow(np.log10(open_ + 1e-6)); ax[0].set_title('open') w = np.ones((n, n)); w[5:36, c + 4] = 0.0 ax[1].imshow(np.log10(field(P, w, 15) + 1e-6)); ax[1].set_title('a wall') plt.show()