import numpy as np
from scipy.integrate import dblquad
import matplotlib.pyplot as plt
c = 343.0 # Speed of sound, SI unit
rho = 1.204 # Density of air, 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 vfunc(v_in=1.0, f=1000.0):
""" Function to give the surface velocity of the piston.
"""
# Normalize the piston velocity to 1 kHz
return v_in / (f / 1000.0)
def intg_func_cart(y, z, pt_meas, k, vp):
""" The integrand function (the part inside the Rayleigh integral) when using
Cartesian coordinates.
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.linalg.norm(pt_meas - np.array([0.0, y, z]))
return vp * np.exp(-1.0j*k*d) / d
def press_cart(pt_meas, f):
""" Compute the Rayleigh integral when using Cartesian coordinates.
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.
This function assumes a rigid piston, and thus the piston velocity `vp` is constant
over its surface.
`pt_meas` are the coordinates at which the complex pressure is returned.
`f` is the piston surface oscillation frequency in Hz
"""
w = 2.0*np.pi*f
k = w/c
# The piston (cone) velocity is obtained by calling vfunc()
vp = vfunc(v_in, f)
intgl_re = dblquad(lambda y, z : np.real(intg_func_cart(y, z, pt_meas, k, vp)),
-0.5*height, 0.5*height, -0.5*width, 0.5*width)
intgl_im = dblquad(lambda y, z : np.imag(intg_func_cart(y, z, pt_meas, k, vp)),
-0.5*height, 0.5*height, -0.5*width, 0.5*width)
return 1.0j*rho*f*(intgl_re[0] + 1.0j*intgl_im[0])
f = 440.0 # signal frequency
v_in = np.sqrt(8.0)/100 # Piston velocity amplitude, m/s
width = 1.0 # Panel width, m
height = 1.0 # Panel height, m
xmeas = np.logspace(np.log10(0.01), np.log10(10.0), 151) # X-coordinates of the measurement points, m
pmeas = np.array([press_cart(np.array([x, 0.0, 0.0]), f) for x in xmeas])
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(xmeas, dB(pmeas) - dB(pmeas[0]))
ax.grid(True, axis='both', which='both')
ax.set_title('On-Axis Sound Pressure Magnitude with Distance')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Normalized Sound Pressure, (dB)')
plt.show()
. This is a some problem for the integration of a plain wave of sub system and more or less point sourse mains. I would like to try line array mains, but I don't have the resources (all of them)) to do so.
