import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
c = 343.0 # Speed of sound, SI unit
def dB(x):
""" dB(x) = 20 log10( abs(x) + 5 eps )
The added small positive number is to prevent the nuisance error when `x` = 0
"""
return 20.0*np.log10(np.abs(x) + 5.0*np.finfo(np.float64).eps)
def intg_func(r, x, k):
""" Integrand function for the Rayleigh integral to compute the on-axis response
from a rigid oscillation circular piston.
Note that this function returns a complex number, and `scipy.integrate.dblquad`
works only with real numbers. The real and imaginary parts of the integral
need to be separately computed with `scipy.integrate.dblquad`.
"""
d = np.sqrt(r**2 + x**2)
return 2*np.pi* r * np.exp(-1.0j*k*d) / d
def press_onaxis(x, k):
""" Compute the on-axis response from a rigid oscillating circular piston in a
sinusoidal motion perpendicular to its axis using the Rayleigh integral method.
The piston motion is normalized to have a volume velocity of 1.
Because the `scipy` `dblquad` double integral function works only with real numbers,
to integrate a function which return complex values, the real and imaginary parts
of the integral need to be computed separately, and the results combined back
together to return the complex values.
`x` is the on-axis distance to the center of the piston.
`k` is acoustic wavenumber. k = 2*pi*f/c, where c is the speed of sound.
"""
intgl_re = quad(lambda r : np.real(intg_func(r, x, k)),
0.0, 0.5*piston_diameter)
intgl_im = quad(lambda r : np.imag(intg_func(r, x, k)),
0.0, 0.5*piston_diameter)
# Multiply by 1.0j to get the correct phase response since we aren't computing
# the pressure from piston surface velocity but from volume velocity.
# Not important if we are only interested in pressure magnitude.
return 1.0j*(intgl_re[0] + 1.0j*intgl_im[0])
def plot_press_vs_dist(ka, xmeas):
""" Plot a normalized on-axis sound pressure (in dB) vs distance curve for
`ka` (ka = product of acoustic wavenumber and piston radius) at
distances `xmeas`.
"""
f = (ka/(0.5*piston_diameter)) * c/(2*np.pi)
pmeas = np.array([press_onaxis(x, ka/(0.5*piston_diameter)) for x in xmeas])
label = '{:.0f} Hz, ka = {:.1f}'.format(f, ka)
ax.semilogx(xmeas, dB(pmeas), label=label)
piston_diameter = 0.115
xmeas = np.logspace(np.log10(0.001), np.log10(10.0), 201) # X-coordinates of the measurement points, m
fig, ax = plt.subplots(figsize=(8, 5))
for ka in [0.1, 0.5, 1.0, 5.0, 10.0]:
plot_press_vs_dist(ka, xmeas)
ax.grid(True, axis='both', which='both')
ax.set_title('On-Axis Sound Pressure vs Distance, Piston Diameter = {:.3f} m'.format(piston_diameter))
ax.legend(loc='upper right')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Normalized Sound Pressure, (dB, Arbitrary)')
ax.set_ylim(-60, 0)
fig.tight_layout()
plt.show()