Analysis

Expressing (3.13) in cylindrical coordinates for a cylindrical duct of radius $R$:

$$Z_{nm} = \frac{i\omega\rho}{2\pi S^2} \int\limits_{0}^{2... ...dr_0 \psi_m(r_0,\theta_0) \psi_n(r,\theta) \frac{e^{-ikh}}{h}$$ (3.14)

where

$$h = [r^2 + r_0^2 - 2rr_0 \cos(\theta-\theta_0)]^{\frac{1}{2}}.$$ (3.15)

As discussed in section 2.4.1, we will be treating cylindrically symmetric modes only. From equation (2.48), the mode profile on the surface is

$$\psi_{n} = \frac{J_0(\gamma_n r/R)}{J_0(\gamma_n)}$$ (3.16)

where $\gamma_n$ is the $n$th zero of the Bessel function $J_1$ and is tabulated in appendix A.

Now we will give $e^{-ikh}/h$ in terms of $\tau$, a dummy variable of integration. We will show that all the other variables of integration will then have an analytic solution. Sonine's integral from Watson [52] p.416 gives:

$$\frac{e^{-ikh}}{h} = k \int\limits_0^{\infty} \tau (\tau^2-1)^{-\frac{1}{2}} J_0(\tau k h) d \tau.$$ (3.17)

The integrand is imaginary when $\tau<1$ and real when $\tau>1$. Care must be taken when choosing the sign of $( \tau^2 - 1)^{-\frac{1}{2}}$ with the negative and imaginary interpretation taken here when $\tau<1$. Notice that we have been using the opposite sign convention from Zorumski [37] for the imaginary part throughout because we are assuming a time factor of $e^{i\omega t}$ rather than $e^{-i \omega t}$.

Neumann's addition formula [52] p.358 is

$$J_0(\tau k h) = \sum\limits_{q=-\infty}^{\infty} J_q(\tau k r) J_q(\tau k r_0) e^{iq(\theta-\theta_0)},$$ (3.18)

which can be substituted into (3.17) to give

$$\frac{e^{-ikh}}{h} = k \sum\limits_{q=-\infty}^{\infty}e^{... ...tau^2-1)^{-\frac{1}{2}} J_q(\tau k r) J_q(\tau k r_0) d \tau.$$ (3.19)

Now substituting (3.19) into (3.14) we get

$$Z_{nm} = \frac{i\rho c}{2\pi S^2} \int\limits_{0}^{2 \pi} d\theta... ...limits_{0}^{2 \pi} d\theta_0 \int\limits_{0}^{R} r_0 dr_0 \psi_m(r_0) \psi_n(r)$$

$$\times k^2 \sum\limits_{q=-\infty}^{\infty} e^{iq(\theta-\theta_0... ...s_0^\infty \tau (\tau^2-1)^{-\frac{1}{2}} J_q(\tau k r) J_q(\tau k r_0) d \tau.$$ (3.20)

Note that $\int_0^{2 \pi} d\theta e^{i q \theta}= 0$ unless $q=0$, in which case it is equal to $2 \pi$. Integrating by $\theta$ and $\theta_0$ then gives a factor of $(2\pi)^2$. When rearranged, the integral can be reduced to

$$Z_{nm} = \frac{i\rho c}{S} \int\limits_0^\infty \tau (\tau^2-1)^{-\frac{1}{2}} D_n(\tau) D_m(\tau) d \tau,$$ (3.21)

where

$$D_n(\tau) = \frac{k \sqrt{2}}{R} \int\limits_0^R r J_0(\tau k r) \psi_n(r) dr.$$ (3.22)

The integration in equation (3.22) can be found analytically (see equation (A.1) in appendix A):

$$D_n(\tau) = \frac{-\sqrt{2} \tau J_1(\tau k R)} {(\frac{\gamma_n}{kR})^2 - \tau^2}.$$ (3.23)

The four dimensional integral has now been reduced to a one dimensional integration, with the variable $h$ (and therefore the singularity mentioned at the end of section 3.3.1) eliminated. Noticing that in equation (3.21) the integral is real for $0<\tau<1$ and imaginary for $1<\tau<\infty$, we split the integral into real and imaginary parts with variables $\tau = \sin{\phi}$ and $\tau = \cosh{\xi}$ respectively.

$$Z_{nm} = \frac{i\rho c}{S} \int\limits_0^{\frac{\pi}{2}} \sin{\ph... ...sin^2{\phi}-1)^{-\frac{1}{2}} D_n(\sin{\phi}) D_m(\sin{\phi}) \cos{\phi} d \phi$$

$$+ \frac{i\rho c}{S} \int\limits_0^\infty \cosh{\xi} (\cosh^2{\xi}-1)^{-\frac{1}{2}} D_n(\cosh{\xi}) D_m(\cosh{\xi}) \sinh{\xi} d \xi.$$ (3.24)

Now $\cosh^2{\xi}-1 = \sinh^2{\xi}$ and $\sin^2{\phi}-1 = -\cos^2{\phi}$. Remembering that the negative imaginary interpretation should be taken for the resulting $(-\cos^2{\phi})^{-\frac{1}{2}}$ we get

$$Z_{nm} = \frac{\rho c}{S} \int\limits_0^{\frac{\pi}{2}} \sin{\phi} D_n(\sin{\phi}) D_m(\sin{\phi}) d \phi$$

$$+ \frac{i\rho c}{S} \int\limits_0^\infty \cosh{\xi} D_n(\cosh{\xi}) D_m(\cosh{\xi}) d \xi.$$ (3.25)

The first integral can be performed by numerical integration using Simpson's rule or an equivalent. In the second integral, however, the range extends to infinity. The integrand is an oscillatory function of $\xi$ whose amplitude of oscillation decays exponentially to $10^{-6}$ typically at around $\xi=10$. Numerical integration can then be performed from 0 and 10 without incurring any significant numerical errors.