Radiation patterns in acoustics describe important aspects of the acoustic behavior of a given radiator shape. Basically, an opening will function as a high-pass filter, with a reduction or amplification of the effective transmission and/or absorption of elements. Also contains an efficient Struve function.

Download C# source code (version 3)

Github project

## Introduction

I think I know what you are thinking: what is this about, and is it interesting to me? Well, probably more than you think, since they touch on many aspects of effective computational methods. If you are an acoustician, the answers you can bring out could give you a vital insight into the behavior of an acoustic wave as it hits an opening.

As the screenshot shows you, this is a normalized impedance that will give you a reactance loss, or, in fact, sometimes, a gain due to the energy that hits the surface compared to a normal incident wave. The impedance will also come with a phase displacement due to the oscillating mass, and here an end correction is often applied at low frequency. This is due to the fact that the oscillating mass is not only over the designated area but also in its surroundings.

For this project, I will rely on some mathematics functions that are found in Math.NET in the NuGet package. There are special functions that are missing from that package that I will supplement with my own.

## Background: Circular baffle

The mathematical development of this started off well over 100 years ago, with the help of some of the earliest adopters and innovators in mathematical physics. Hermann von Helmholtz showed that the sound pressure from a vibrating surface was given by the two integrals:

$ p(R,t) = \oint\limits_{S}\, \frac{e^{i(\omega t - kr)}}{4 \pi r} \left( \frac{\partial p}{\partial n} \right)_S \, \mathrm{d}s - \oint\limits_{S}\, p(S) \cdot \frac{\partial}{\partial n} \left(\frac{e^{i(\omega t - kr)}}{4 \pi r} \right) \, \mathrm{d}s$

The importance of sound mathematical reasoning by applying various conservation laws and logic shows how precise this kind of knowledge can be. Remember that in 1900, no electrical machinery excisted to do measuring, so most of these logics were derived theoretically with the help of the most basic experiments. They did some calculations and made logical conclusions from them.

Given the velocity, we can find the pressure that contributes to the first integral. The second integral requires that the pressure over the whole surface be known, and this cannot be known in general. But when placed in an infinite baffle, a heavy wall that extends indefinitely to each side and is adjacent to the vibrating surface, this integral vanishes. Along now comes one of the other great acousticians from this time period, Lord Rayleigh. He derives the following integral:

$p(R,t) = i \frac{\rho_0 c_0 k}{2 \pi R} \oint\limits_{S}\, u_n(S) \cdot e^{i (\omega t - kr)} \, dS$

For a circle, this will give

$ p(R, \theta ,t) = i \frac{\rho_0 c_0 k a^2}{2R} \hat{u} \cdot e^{i(\omega t - kR)} \left( \frac{2 \cdot J_1(ka \cdot \sin \theta)}{ka \cdot \sin \theta} \right) $

This can be directly used to calculate the radiation impedance:

$ 1 - \frac{2 \cdot J_1(2ka)}{2ka} + i \frac{2 \cdot H_1(2ka)}{2ka} $

The code is simply:

public static Complex CircularBaffle(double k_0, double a, double Z_0 = 1)
{
if (k_0 == 0)
return Complex.Zero;
double k0a = k_0 * a;
Complex i = new Complex(0, 1);
return Z_0*(1 - MathNet.Numerics.SpecialFunctions.BesselJ(1,2*k0a)/k0a + i * SpecialFunctions.StruveH(1,2*k0a)/k0a);
}

Next, we need a good implementation of the Struve function, called StruveH.

## Efficient Implementation of the Struve Function

I have come to the realization that this section could be a whole article by itself. It touches on a subject that I probably should be writing an article about: the asymptotic expansion. I decided to just do a good enough implementation with 9 significant digit accuracy from the range z = 0 to z = 20 000. It was also valid above this range, but that was not so interesting for the acoustical implementation. This search for accuracy also led me to implement tests to check that the functions did what they were supposed to do. I simply used Wolfram Alpha to generate the correct values that I used to check against.

Most implementations of the Struve function shared over the web seem to be way too complicated. Given that you already have access to a good implementation of the Neumann function, the implementation is indeed very simple. But unfortunately, it is not very accurate, and for lower orders of the Struve function, the worse it gets.

First, for the small input argument, the implementation with a Taylor series is more than good enough. The formula looks like the one below:

$ \textbf{H}_\alpha (z) = \sum_{m=0}^{\infty} \frac{(-1)^m}{\Gamma \left(m + \frac{3}{2} \right) \cdot \Gamma \left(m + \alpha + \frac{3}{2}\right) } \left( \frac{z}{2} \right)^{2m + \alpha + 1} $

The implementation, which I did:

private static Complex StruveHTaylor(double n, Complex z)
{
Complex TermResult = Complex.Zero;
List<Complex> SummationTerms = new List<Complex>();
if (z == Complex.Zero) { return Complex.Zero; }
double MaxIteration = 25d;
Complex Log_z2 = Complex.Log(z * 0.5d);
Complex iPi = new Complex(0, Math.PI);
Complex error;
double tol = 1e-10;
for (double m = 0; m < MaxIteration; m += 1d)
{
Complex LogarithmicTermEvaluation =
m * iPi
+ (2 * m + n + 1) * Log_z2
- MathNet.Numerics.SpecialFunctions.GammaLn(m + 1.5d)
- MathNet.Numerics.SpecialFunctions.GammaLn(m + n + 1.5d);
TermResult += Complex.Exp(LogarithmicTermEvaluation);
SummationTerms.Add(TermResult);
Complex Summation = EpsilonAlgorithm(SummationTerms.ToArray());
if (m > 0)
{
error = SummationTerms.Last() - Summation;
if (tol > Complex.Abs(error))
{
return Summation;
}
}
}
return EpsilonAlgorithm(SummationTerms.ToArray());
}

I will here use the formula from "The asymptotics of the Struve function H_ν(z) for large complex order and argument " by R. B. Paris. I start off with the work of R. B. Dingle:

$H_n (z) - Y_n (z) = \frac{2 \left( \frac{1}{2} z \right)^n}{\sqrt{\pi} \, \Gamma \left( n + \frac{1}{2} \right)} \cdot \int_0^\infty{e^{-z \cdot x} \cdot \left( 1 + x^2\right)^{n - \frac{1}{2}} \, dx}$

This formula is only valid for $|arg \, z| < \frac{\pi}{2}$. However, it could easily be extended by the use of the relation:

$H_n(z \cdot e^{\pi m i}) = e^{\pi m i (n+1) } \cdot H_n(z) \qquad m = \pm1, \pm2, \cdots$

This integral could be calculated with the Gaussian-Lagerre quadrature, since the integral is already of this form. I did this in Matlab but only got asymptotic convergence at very high values (from z = 10 000). So, instead, this integral is expanded with a reverse power series. All this is based on the work and articles by R. B. Dingle, and you can even take a look at some of his work in the book (see the references in the Paris article). The book is out of print but can be downloaded from Professor Sir Michael Victor Berry, FRS home page. The formula Dingle comes up with looks like this:

$H_n (z) - Y_n (z) = \frac{(0.5 \cdot z)^{n-1}}{\sqrt{\pi} \, \Gamma(n+0.5)}\left( 1 - \frac{1 \cdot (1- 2n)}{x^2} + \frac{1(1-2n)3(3-2n)}{x^4} - \frac{1(1-2n)3(3-3n)5(5-2n)}{x^6} + \cdots \right)$

Now, we exchange the integral with a reverse power series, where 10 of the coefficients are given in Paris article.

private static Complex[] C_k(Complex q)
{
Complex q2 = q * q;
Complex q3 = q2 * q;
Complex q4 = q3 * q;
Complex q5 = q4 * q;
Complex q6 = q5 * q;
Complex q7 = q6 * q;
Complex q8 = q7 * q;
Complex q9 = q8 * q;
Complex q10 = q9 * q;
return new Complex[]{
1d,
2d * q,
6d * q2 - 0.5d,
20d * q3 - 4d * q,
70d * q4 - 45d / 2d * q2 + 3d / 8d,
252d * q5 - 112d * q3 + 23d / 4d * q,
924d * q6 - 525d * q4 + 301d / 6d * q2 - 5d / 16d,
3432d * q7 - 2376d * q5 + 345d * q3 - 22d / 3d * q,
12870d * q8 - 21021d / 2d * q6 + 16665d / 8d * q4 - 1425d / 16d * q2 + 35d / 128d,
48620d * q9 - 45760 * q7 + 139139d / 12d * q5 - 1595d / 2d * q3 + 563d / 64d * q,
184756d * q10 - 196911d * q8 +61061d * q6 - 287289d / 48d * q4 + 133529d / 960d * q2 - 63d / 256d
};
}

Suppose you have a good enough (professional) implementation of the Neumann function. Then the Taylor series, together with the asymptotic expansion at $x \to \infty$ will give you a very straight-forward implementation that is accurate by 9 digits in the range 0–20 000 (and way above, but I did not test it as I did not need it).

public static Complex StruveHnNeumannYn(double n, Complex z)
{
Complex Result = Complex.Zero;
Complex ConstantTerm = (n - 1) * Complex.Log(0.5 * z)
- 0.5 * Complex.Log(Math.PI)
- MathNet.Numerics.SpecialFunctions.GammaLn(n + 0.5);
Complex[] Ck = C_k(n / z);
for (int k = 0; k < Ck.Length; k++)
{
Complex item = ConstantTerm + Complex.Log(Ck[k])
+ MathNet.Numerics.SpecialFunctions.GammaLn(k + 1)
- k * Complex.Log(z);
Result += Complex.Exp(item);
}
return Result;
}

Putting it all together:

public static Complex StruveH(double n, Complex z)
{
if (Complex.Abs(z) < 20)
return StruveHTaylor(n, z);
else
return StruveHnNeumannYn(n, z) + MathNet.Numerics.SpecialFunctions.BesselY(n, z);
}

This implementation will be valid for nearly any range that you could possibly want in practice. MathNet does have an implementation of the modified Struve function, which is a bit unfortunate since you can express the modified Struve function with the complex argument of the Struve function with some phase shift and negative imaginary number. Written out in a mathematical formula:

$ -i \cdot e^{-in \pi / 2} \textbf{H}_\alpha (i \cdot z) = \textbf{L}_\alpha (z) $

## Elliptic piston

A formula for the radiation impedance of an elliptic piston, also derived by Mechel, looks like this:

$ 1 - \frac{\pi}{2} k^2 \cdot a \cdot b \cdot \int_0^{\pi/2} \frac{J_1(B)}{B^3} \, d \phi + i \cdot \frac{\pi}{2} k^2 \cdot a \cdot b \cdot \int_0^{\pi/2} \frac{H_1(B)}{B^3} \, d \phi$

The integral is dependent en the "elliptic function" B:

$B = k \cdot a \cdot \sqrt{\cos^2(\phi) + \beta^2 \cdot \sin^2(\phi)}$

The general method for solving this is to expand the integrand into a power series and iterate. For the Series from the BesselJ1 funtion:

$Z_(J_1) = \beta \left( \frac{(k\cdot a)^2}{2} + \sum_{n=2}^{m}{ C_n(J_1) \cdot I_n(J_1) }\right)$

$C_1(J_1) = \frac{(k \cdot a)^2}{2}$

$C_n(J_1) = \frac{-(k \cdot a)^2}{n \cdot (n + 1)} \cdot C_{n-1}(J_1)$

$ I_0(J_1) = \frac{1}{\beta} \qquad I_1(J_1) = 1$

$I_n(J_1) = \frac{2n - 3}{2n-2} (1+\beta^2) \cdot I_{n-1} (J_1) + \frac{2n - 4}{2n-2} \beta^2 \cdot I_{n-2} (J_1)$

For the imaginary part containing the StruveH1 function

$Z(H_1) = \frac{4 \cdot k \cdot b}{\pi^2} \left( \frac{4}{3} \cdot I_0(H_1) - \frac{16}{45}(k \cdot a)^2 \cdot I_1(H_1) + \sum_{n=2}^{m}{ C_n(H_1) \cdot I_n(H_1) }\right)$

$C_1(H_1) = -\frac{16 \cdot (k \cdot a)^2}{45}$

$C_n(H_1) = \frac{-4 \cdot (k \cdot a)^2}{(2n +1 )\cdot (2n + 3)} \cdot C_{n-1}(H_1)$

$ I_0(H_1) = EllipticK(1 - \beta^2) \qquad I_1(H_1) = EllipticE(1-\beta^2)$

$I_n(H_1) = \frac{2n - 2}{2n-1} (1+\beta^2) \cdot I_{n-1} (H_1) + \frac{2n - 3}{2n-1} \beta^2 \cdot I_{n-2} (H_1)$

There is a problem with these methods, and that has to do with the rapidly growing powers of the alternating series, so I will eventually struggle if the values of x gets high enough. Like always, we need an asymptotic expansion when $x \to \infty $.

I know I have a very good algorithm for the circular baffle, so I will simply use that for the asymptotic expansion. I simply add the two different normalized impedances with the two different radiuses from the elliptic disk and take the average of them. This trick only works well enough when the natural wave number coincides with the average diameter of the elliptic disk, and this roughly happens when the normalized real part of the elliptic piston becomes 1. Ideally, we should check if the iterations become unruly and switch then, but that will have to be for a later time, and that is, in any case, not guaranteed to give you an accurate result.

public static Complex EllipticBaffle(double k, double a, double b, double Z_0 = 1)
{
if (k == 0)
return Complex.Zero;
double ka = k * a;
double ka2 = ka * ka;
double beta = Math.Min(a, b) / Math.Max(a, b);
double Limit = Math.Ceiling(4 * ka) + 4d;
double MaxLimit = 80;
if (Limit < MaxLimit)
{
double Z_r = beta * (EllipticSumZ_r(ka2, beta, Limit));
double Z_i = 4 / Math.PI / Math.PI * ka * beta * EllipticSumZ_i(ka2, beta, Limit);
return Z_0 * new Complex(Z_r, Z_i);
}
else
{
Complex LongSideBaffel = CircularBaffle(k, a);
Complex ShortSideBaffel = CircularBaffle(k, b);
return Z_0 * (LongSideBaffel + ShortSideBaffel) / 2;
}
}

##

## Field-Excited Rectangular radiator

The following derivation follows Mechel's approach with some of my own tweaking to the final implementation. Given a rectangular element placed in a hard baffle wall with an obliquely incident plane wave moving towards the surface, the term field excited refers to a vibration pattern set in motion by a wave propagating towards the opening.

The object here is a rectangle A with side lengths a and b on a hard baffle wall. The recatangle is excited by a plan wave with a polar incident angle of $\theta$ and an azimuthal angle of $\phi$. The average velocity pattern on the surface is:

$\begin{equation} \begin{aligned} &v(x,y) = V_0 \cdot e^{-i(k_x x+k_y y)} \\ & k_x = k_0 \cdot \sin(\theta) \cdot \cos(\phi) = k_0 \cdot \mu_x \\ & k_y = k_0 \cdot \sin(\theta) \cdot \sin(\phi) = k_0 \cdot \mu_y \end{aligned} \label{eq:GeneralVelocityPatternOnA} \end{equation}$

Following Mechel's calculations, we have the fact that $|v(x,y)| = \text{constant}$ on the surface $A$ and that the radiation impedance follows from the average field impedance. The Fourier transform of the velocity distribution leads to a form with a double integration:

$\begin{equation} \begin{aligned} V(k_1,k_2) &= \iint\limits_A{v(x_0 , y_0) e^{-i (k_1 x + k_2 y_0)}} \, dx_0 \, dy_0 \\ &= V_0 a b \cdot \frac{\sin((k_1+k_x)a/2)}{(k_1+k_x)a/2} \cdot \frac{\sin((k_2+k_y)b/2)}{(k_2+k_y)b/2} \end{aligned} \label{eq:FTofVelocityPattern} \end{equation}$

using $\alpha_1 = k_1/k_0$ and $\alpha_2 = k_2/k_0$:

$\begin{equation} \frac{Z_r}{Z_0} = \frac{k_0 a \cdot k_0 b}{4 \pi^2} \iint\limits_{-\infty}^{\infty}\left( \frac{\sin(\mu_x - \alpha_1)k_0 a/2}{(\mu_x - \alpha_1)k_0 a/2} \right)^2 \left( \frac{\sin(\mu_y - \alpha_2)k_0 b/2}{(\mu_y - \alpha_2)k_0 b/2} \right)^2 \frac{d \alpha_1 \, d\alpha_2}{\sqrt{1-\alpha_1^2 -\alpha_2^2}} \label{eq:FTofVelocityPatternSubstitute} \end{equation}$

This could be translated into:

$\begin{equation} \frac{Z_r}{Z_0} = \frac{2i}{\pi k_0 a \cdot k_0 b} \int_{x=0}^{k_0 a} dx \int_0^{k_0 b} (k_0 a-x)(k_0 b-y)\cos(\mu_x x) \cos(\mu_y y) \frac{e^{-i \sqrt{x^2 +y^2}}}{\sqrt{x^2 +y^2}} dy \label{eq:FTofVelocityPatternSubstitute2} \end{equation}$

The integral is not really suited for numerical integration as the limits are $k_0 \cdot a $ and $k_0 \cdot b$, which can be very large values indeed. Therefore, it is better to perform a substitution of variables and translate the double integral into:

$\begin{equation} \frac{Z_r}{Z_0} = \frac{2i}{\pi k_0 a \cdot k_0 b} \left[ \int_0^{\arctan(b/a)}{I(k_0 a/\cos(\varphi)) \, d \varphi} + \int_{\arctan(b/a)}^{\pi/2}{I(k_0 b/\sin(\varphi) \, d\varphi} \right] \label{eq:RadiationPatternWithIntermediateIntegral} \end{equation}$

The intermediate integral is given by:

$\begin{equation} I(R) = \int_0^R{ (U + V \cdot r + W \cdot r^2) \cos(\alpha r) \cos(\beta r) \cdot e^{-i r} \, dr} \label{eq:IntermediateIntegral} \end{equation}$

With the variables:

$\begin{equation} \begin{aligned} U &= k_0 a \cdot k_0 b \\ V &=-(k_0 a \sin(\phi) + k_0 b \cos(\phi)) \\ W &= \sin(\phi) \cos(\phi) = \sin(2 \cdot \phi) \\ \alpha &= \mu_x \cos(\phi) = k_0 \cdot \sin(\theta_i) \cdot \cos( \phi_i ) \cdot \cos(\phi) \\ \beta &= \mu_y \sin(\phi) = k_0 \cdot \sin(\theta_i) \cdot \sin( \phi_i ) \cdot \sin(\phi) \end{aligned} \end{equation}$

The intermediate integral could be integrated by substituting the cosine terms with the identity:

$\begin{equation} \cos(x) = \frac{e^{-i \cdot x} + e^{i \cdot x}}{2} \end{equation}$

Since every term is multiplied by $\cos(x)$ we can take the constant $1/4$ out of the integral, which can be a little simplified:

$\begin{equation} \begin{split} I(R) &=\frac{1}{4} \int_0^R (e^{ir (- \alpha - \beta -1)} +e^{ir (\alpha - \beta -1)}+e^{ir (- \alpha + \beta -1)}+e^{ir (\alpha + \beta -1)} ) \\ & \qquad \qquad \cdot (U + V \cdot r + W \cdot r^2) \, dr \end{split} \label{eq:IntermediateIntegral_rewritten} \end{equation}$

We now introduce a variable called $c_n$ which has the following 4 forms:

$\begin{equation} \begin{aligned} c_1 &= - \alpha - \beta -1 \\ c_2 &= \alpha - \beta -1 \\ c_3 &= - \alpha + \beta -1 \\ c_4 &= \alpha + \beta -1 \end{aligned} \end{equation}$

Given that $n=1,2,3,4$ the general term for $c_n$:

$\begin{equation} \begin{aligned} c_n &= (-1)^n \cdot \alpha + (-1)^{\left(\frac{n^2 + n}{2}\right)} \cdot \beta -1 \\ &= (-1)^n \cdot \alpha + i^{\left(n(n+1)\right)} \cdot \beta -1 \end{aligned} \label{eq:c_n_generalform} \end{equation}$

There is also the possibility to combine the original values and simplify the expression for $c_n$:

$\begin{equation} c_n = (-1)^{n}\cdot k_0 \cdot \sin(\theta) \cdot \cos{\left( \phi + i^{(n(n-1))} \cdot \varphi \right)} \label{eq:c_nOriginalExpressionSimplifyed} \end{equation}$

However, the implementation is much simpler in code with the use of the bitwise and operator:

double mod1 = (1 & n) == 0 ? 1 : -1;
double mod2 = (2 & n) == 0 ? 1 : -1;
Complex C_n = mod1 * alpha(j) + mod2 * beta(j) - 1;

Resulting in:

$\begin{equation} I(R,c_n) = \frac{1}{4}\sum_{n=1}^{4} \left(\int_0^R{e^{ir \cdot c_n } \cdot U \, dr} + \int_0^R{e^{ir \cdot c_n } \cdot V \cdot r \, dr} + \int_0^R{e^{ir \cdot c_n } \cdot W \cdot r^2 \, dr} \right) \label{eq:IntermediateIntegralWithSummationOverCn} \end{equation}$

Using standard integration and integration by parts, these are easily solved:

$\begin{equation} \int_0^R{e^{i c_n r} U \, dr} = - \frac{i U e^{i c_n R} }{c_n} + \frac{i U}{c_n} \label{eq:IntermediateIntegral_part1} \end{equation}$

$\begin{equation} \int_0^R{e^{i c_n r} V \cdot r \, dr} = V e^{i c_n R} \left( \frac{1}{c_n^2} - \frac{i R}{c_n}\right) - \frac{V}{c_n^2} \label{eq:IntermediateIntegral_part2} \end{equation}$

$\begin{equation} \int_0^R{e^{i c_n r} W \cdot r^2 \, dr} =W e^{i c_n R} \left( \frac{2 i }{c_n^3} + \frac{2R}{c_n^2} - \frac{i R^2}{c_n} \right) - \frac{2 i W}{c_n^3} \label{eq:IntermediateIntegral_part3} \end{equation}$

Sorting out and rearranging gives the following:

$\begin{equation} I(R) = \frac{1}{4 \cdot C_n^3} \sum_{n=1}^{4} e^{i \cdot C_n \cdot R} \cdot \left( - U \cdot i C_n^2 + V* \left(C_n - i \cdot C_n^2 \cdot R \right) + W \cdot \left( -i \cdot C_n^2 \cdot R^2 + 2 \cdot C_n \cdot R + 2 \cdot i \right) \right) + U \cdot i \cdot C_n^2 - V \cdot C_n - W \cdot 2 \cdot i \end{equation}$

The code that generates this integral can be implemented as such:

public static Complex I_R(Complex R, Complex C_n, Complex U, Complex V, Complex W)
{
Complex C_n2 = C_n * C_n;
Complex C_n3 = C_n2 * C_n;
Complex i = new Complex(0, 1);
Complex expR = Complex.Exp(i * C_n * R);
Complex result =
( expR * ( - U * i * C_n2
+ V * (C_n - i * C_n2 * R)
+ W * (-i * C_n2 * R * R + 2 * C_n * R + 2 * i)
)
+( U * i * C_n2
- V * C_n
- W * 2 * i )
) / (C_n3 * 4);
return result;
}

The two integrals can now easily be implemented by a Cauchy sum. Which is the same as the Riemann sum, only that a Cauchy sum uses finitely many elements while Riemann sets the number of elements to infinity.

for (double i = 0; i < arctan; i += deltaPhi)
{
for (int n = 0; n < 4; n++)
{
double mod1 = (1 & n) == 0 ? 1 : -1;
double mod2 = (2 & n) == 0 ? 1 : -1;
Complex C_n = mod1 * alpha(i) + mod2 * beta(i) - 1;
result += I_R(k0a / Math.Cos(i), C_n, U, V(i), W(i));
}
}
for (double j = arctan; j < Math.PI / 2; j += deltaPhi)
{
for (int n = 0; n < 4; n++)
{
double mod1 = (1 & n) == 0 ? 1 : -1;
double mod2 = (2 & n) == 0 ? 1 : -1;
Complex C_n = mod1 * alpha(j) + mod2 * beta(j) - 1;
result += I_R(k0b / Math.Sin(j), C_n, U, V(j), W(j));
}
}

If you are tempted to use some Newton-Coates-style summation, like Gaussian quadrature, you should remember that the requirement for using this is a realistically well-behaved function. There is no guarantee of that here, at least not in advance. That is why I opted for a more safe, straight-forward summation. This code also seems to be very stable even for massively large openings, like 100 times 100 meters.

## Rectangular Piston

There is a direct implementation for a Rectangular Piston as well, and it is included as a reference. But don't use it, instead use the field excited rectangular with an incoming wave which is normal incident theta = 0 and phi = 0. This is much more stable. The code just so we are crystal clear here:

List<Complex> RetangularPistonResult = new List<Complex>();
for (double k = 0; k <= UpperMode; k += deltaWaveNumber)
{
RetangularPistonResult.Add(Acoustics.Radiation.FieldExcited.Rectangular(k, 0, 0, a, b));
}

## Wide and Narrow strips

With the updated Struve function and the use of a rectangular piston with width 10 times longer than the height, these are now also accurate for a very large range.

## End corrections

The imaginary part of the radiation impedance really tells you about the mass that is oscillating in the opening. At lower frequencies, the mass is actually larger than the mass only contained in the opening, so we have something called end corrections. These are dependent on the shape of the opening and the available surrounding area that can oscillate, and they must be added on for a more realistic imaginary impedance.

## Further work

There are many more implementations of different shapes of pistons, like spheres, cylinders, doughnuts, etc. There is also the case where several small perforated openings work in tandem. But for those, I will leave it to you to investigate the references for further study.

## References

- Formulas of Acoustics, 2
^{nd} ed. by F. P. Mechel et. al. - Building Acoustics by Tor Erik Vigran

I also used some material for articles that are referenced by those books, so I did not include them here.

I would also like to thank the late F. P. Mechel for his kind help and for sending his own code used in Formulas of Acoustics.

## History

- 14
^{th} January, 2024: Initial version - 17
^{th} January, 2024: Updated explanations and code - 25
^{th} January, 2024: Updated Struve to 9-digit accuracy and added tests to mathematical functions.

This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.