text:
\section{Introduction}
	The spectral method is a highly precise numerical technique similar to the finite difference method (FDM) and the finite element method (FEM) \cite{Boyd2001}. Compared with the FDM and FEM, however, the spectral method is a global discrete approach, meaning its basis functions are defined in the entire solution domain, whereas the basis functions of the FDM and FEM are defined on the grid points and in the elements, respectively; hence, the spectral method is more accurate than either the FDM or the FEM \cite{Jshen2011}. Strictly speaking, the spectral method is a weighted residual method, but the basis functions selected by the spectral method---unlike those of the general weighted residual method---are sets of orthogonal polynomials, and the orthogonality of these basis functions ensures that the spectral method converges quickly \cite{Canuto2006}. The research conducted by Orszag and Gottilieb illustrated that the numerical error of the spectral method decreases exponentially with an increase in the truncation order \cite{Orszag1972,Gottlieb1977}. Consequently, benefiting from its high precision, the spectral method has been widely used in many numerical simulations of scientific and engineering problems \cite{Canuto1988,Guoby1998,Tuhw2021c}.

	Accurately and efficiently identifying and locating underwater targets have always been the core pursuit of computational ocean acoustics \cite{Etter2018}. The most critical step in achieving these objectives is accurately calculating the sound field \cite{Jensen2011}, for which the spectral method, being a highly precise numerical technique, has recently been employed. The resulting spectral method-based algorithms and programs have shown excellent performance in the two-dimensional simulation of underwater acoustic propagation.

	In 1993, Dzieciuch wrote a program to calculate the ideal fluid waveguide of the Munk sound speed profile using the Chebyshev--Tau spectral method \cite{aw}. However, this program could calculate only a single-layer water body in which the seawater density does not vary with depth and the seabed is a pressure-release boundary, and thus, it was not capable of solving higher-complexity acoustic propagation problems; nevertheless, this program was the first instance that the spectral method had been used to solve the underwater acoustic propagation problem, and it demonstrated remarkable accuracy. More recently, in 2016, Evans proposed the Legendre--Galerkin spectral method to solve for the normal modes of underwater acoustic propagation in two-layer media and developed the program rimLG \cite{rimLG}. Because rimLG also requires the seabed to be a pressure-release boundary and applies Gaussian integration in each parameter layer, its computational speed becomes quite slow when the acoustic profile data are dense, but it exhibits extremely high accuracy. Moreover, rimLG can simulate underwater acoustic waveguides with attenuation in the sediment, and the diversity of sound waveguides that can be synthesized is much greater than the capability of Dzieciuch's program. Subsequently, in 2019, Sabatini proposed a multidomain spectral collocation method that can accurately calculate the normal modes in the open ocean \cite{Sabatini2019}. This method can accurately compute the half-space boundary of the seafloor, although the size of the matrix eigenvalue problem is doubled.

	Since 2020, we have performed a series of studies to design highly precise numerical algorithms based on the spectral method for solving underwater acoustic propagation problems in a marine environment. In March 2020, we proposed an algorithm based on the Chebyshev--Tau spectral method to solve the underwater acoustic propagation problem in two-layer media \cite{Tuhw2020a,Tuhw2020b} and developed a corresponding simulation software named NM-CT \cite{NM-CT}. Numerical simulations showed that the accuracy of NM-CT is comparable to that of rimLG, and the numerical errors are all on the order of $10^{-14}$, which is much smaller than the order of $10^{-6}$ of the classic FDM; in addition, the computational speed is much faster than that of rimLG. Then, in June 2021, we proposed an algorithm for simulating underwater acoustic propagation in multilayer media based on the Legendre collocation spectral method \cite{Tuhw2021d,MultiLC}. Our model adopts an absorbing layer to simulate the acoustic half-space, which overcomes the shortcomings of Sabatini's algorithm, namely, matrices that are twice as large and a slow computational speed. Additionally, we conducted research on solving parabolic equation models with the spectral method \cite{Tuhw2021a,Tuhw2021b,SMPE} and subsequently proposed a spectral method-based algorithm for solving acoustic propagation in a two-dimensional, range-dependent marine environment \cite{Tuhw2022a,Tuhw2022b}.

	To date, the spectral method has shown remarkable success and potential in designing highly precise numerical algorithms for simulating underwater acoustic propagation. Nevertheless, the existing numerical algorithms for simulating three-dimensional underwater sound fields still suffer from low computational accuracies and slow computational speeds \cite{Lee1995,Linyt2019,Liuw2021}. Accordingly, the development of a highly precise three-dimensional underwater acoustic propagation model based on the spectral method is a research topic of considerable importance.

	In this paper, we develop a highly precise quasi-three-dimensional numerical algorithm for simulating underwater acoustic propagation using the integral transformation technique \cite{Luowy2016}, a stair-step approximation \cite{Evans1983}, and the Chebyshev--Tau spectral method.

	\section{Mathematical Modeling}
	\subsection{Quasi-three-dimensional marine environment}
	\label{section2.1}
	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale = 0.8]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(14,-4.5)--(12,-4.5)--(6.5,-7.5)--(4,-7.5)--(2,-5)--cycle;
			\fill[orange,opacity=0.7] (2,-5)--(4,-7.5)--(6.5,-7.5)--(12,-4.5)--(14,-4.5)--(14,-8.5)--(4,-8.5)--(2,-6)--cycle;
			\fill[brown,opacity=0.7] (2,-6)--(4,-8.5)--(14,-8.5)--(14,-9.5)--(4,-9.5)--(2,-7)--cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-9) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(5,-3.75) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-7);
			\draw[very thick](4,-2.5)--(4,-9.5);
			\draw[very thick](2,-5)--(4,-7.5);
			\draw[very thick](12,0)--(12,-2);
			\draw[very thick](14,-2.5)--(14,-9.5);
			\draw[very thick](12,-2)--(14,-4.5);

			\draw[very thick](2,-5)--(4.5,-5);
			\draw[very thick](4,-7.5)--(6.5,-7.5);
			\draw[very thick](4.5,-5)--(6.5,-7.5);
			\draw[very thick](6.49,-7.5)--(12,-4.5);
			\draw[very thick](4.5,-5)--(10.02,-2);
			\draw[very thick](12,-4.5)--(14,-4.5);
			\draw[very thick](10,-2)--(12,-2);
			\draw[very thick](10,-2)--(12,-4.5);

			\draw[dashed, very thick](2,-6)--(4,-8.5);
			\draw[dashed, very thick](4,-8.5)--(14,-8.5);
			\filldraw [red] (8,-2.7) circle [radius=2.5pt];
			\node at (2.7,-1.3){$o$};
			\node at (9,-5){Water};
			\node at (9,-7.5){Sediment};
			\node at (9,-9){Acoustic half--space};
		\end{tikzpicture}
		\caption{Schematic diagram of the quasi-three-dimensional marine environment.}
		\label{Figure1}
	\end{figure}

	A continental slope is a ubiquitous offshore seabed terrain, and the typical structure of a continental slope is shown in Fig.~\ref{Figure1}. Terrain and ocean acoustic parameters are range-dependent variables in the $xoz$--plane; in other words, parameters such as the sound speed, density and attenuation, as well as terrain, vary across the plane. In contrast, in the $yoz$--plane, the acoustic parameters are $y$-independent. Thus, this scenario is not a real three-dimensional problem but rather a quasi-three-dimensional problem. To accurately simulate this problem physically, we can always place the sound source in the plane of $y=0$; thus, the governing equation of the problem can be written as \cite{Jensen2011}:
	\begin{equation}
		\label{eq.1}
		\begin{aligned}
			&\rho(x,z) \frac{\partial}{\partial x}\left[\frac{1}{\rho(x,z)} \frac{\partial p}{\partial x}\right]+\frac{\partial^{2} p}{\partial y^{2}}+\\
			&\rho(x,z) \frac{\partial}{\partial z}\left[\frac{1}{\rho(x,z)} \frac{\partial p}{\partial z}\right]+k^{2}(x,z) p=-\delta\left(x-x_{s}\right) \delta(y) \delta\left(z-z_{s}\right)
		\end{aligned}
	\end{equation}
	where the harmonic point source is located at a horizontal range $x=x_s$ and a depth $z=z_s$, $k(x,z)=2\pi f/c(x,z)$ is the wavenumber, $f$ is the frequency of the sound source, $c(x,z)$ is the sound speed, and $p\equiv p(x,y,z)$. Since the marine environment is $y$-independent, Eq.~\eqref{eq.1} can be transformed using the following Fourier transform:
	\begin{subequations}
		\begin{align}
			\label{eq.2a}
			\tilde{p}\left(x, k_{y}, z\right) &=\int_{-\infty}^{\infty} p(x, y, z) \mathrm{e}^{-\mathrm{i} y k_{y}} \mathrm{~d} y \\
			\label{eq.2b}
			p(x, y, z) &=\frac{1}{2 \pi} \int_{-\infty}^{\infty} \tilde{p}\left(x, k_{y}, z\right) \mathrm{e}^{\mathrm{i} y k_{y}} \mathrm{d} k_{y}
		\end{align}
	\end{subequations}
	This is accomplished by applying the following operators to both sides of Eq.~\eqref{eq.1}:
	\begin{equation*}
		\int_{-\infty}^{\infty}(\cdot) \mathrm{e}^{-\mathrm{i} k_{y} y} \mathrm{d} y
	\end{equation*}
	From Eq.~\eqref{eq.2a}, Eq.~\eqref{eq.1} can be transformed into the following form:
	\begin{equation}
		\label{eq.3}
		\begin{aligned}
			&\rho(x,z) \frac{\partial}{\partial x}\left[\frac{1}{\rho(x,z)} \frac{\partial \tilde{p}}{\partial x}\right] + \\
			&\rho(x,z) \frac{\partial}{\partial z}\left[\frac{1}{\rho(x,z)} \frac{\partial \tilde{p}}{\partial z}\right]+ \left(k^{2}-k_{y}^{2}\right) \tilde{p}=-\delta\left(x-x_{s}\right) \delta\left(z-z_{s}\right)
		\end{aligned}
	\end{equation}
	For each particular $k_y$, Eq.~\eqref{eq.3} has almost the same form as the governing equation for an infinitely long line source sound field in a two-dimensional plane. The solution of Eq.~\eqref{eq.3} is described in detail in Sec.~\ref{section2.2}. After calculating the sound pressure $\tilde{p}(x,k_y,z)$ in the $k_y$--domain, the sound pressure $p(x,y,z)$ in the three-dimensional domain can be easily calculated through the inverse Fourier transform of Eq.~\eqref{eq.2b}.

	\subsection{Stepwise approximation of a two-dimensional marine environment}
	\label{section2.2}
	The previous section demonstrates that Eq.~\eqref{eq.3} is highly similar in form to the governing equation of the sound field of an infinitely long line source in a two-dimensional plane. Here, let the acoustic pressure be denoted as $P\equiv P(x,z)$. The acoustic governing equation of an infinitely long line source can be written as \cite{Jensen2011}:
	\begin{equation}
		\label{eq.4}
		\rho(x,z) \frac{\partial}{\partial x}\left[\frac{1}{\rho(x,z)} \frac{\partial P}{\partial x}\right]+\rho(x,z) \frac{\partial}{\partial z}\left[\frac{1}{\rho(x,z)} \frac{\partial P}{\partial z}\right]+k^{2} P=-\delta(x-x_s) \delta\left(z-z_{s}\right)
	\end{equation}
	The difference between Eq.~\eqref{eq.4} and Eq.~\eqref{eq.3} is only that the position of $k^2$ becomes $(k^2-k_y^2)$. The solution of Eq.~\eqref{eq.4} is described in Ref.~\cite{Tuhw2022b}. Here, we deduce the solution of Eq.~\eqref{eq.3} as the solution of the two-dimensional line source sound field.

	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale = 0.8]
			\node at (1.8,0){$0$};
			\fill[brown,opacity=0.7] (14,-8) rectangle (2,-6.7);
			\fill[cyan,opacity=0.7] (14,0)--(14,-1.7)--(12,-1.7)--(4,-5.7)--(2,-5.7)--(2,0)--cycle;
			\fill[orange,opacity=0.7] (14,-1.7)--(12,-1.7)--(4,-5.7)--(2,-5.7)--(2,-6.7)--(14,-6.7)--cycle;
			\draw[very thick, ->](2,0)--(14.5,0) node[right]{$x$};
			\draw[very thick, ->](2.02,0)--(2.02,-8.5) node[below]{$z$};
			\draw[very thick](2,-5.7)--(4,-5.7);
			\draw[dashed, very thick](4,-5.7)--(5,-5.7);
			\draw[dashed, very thick](5,-5.7)--(5,-4.7);
			\draw[dashed, very thick](5,-4.7)--(6.8,-4.7);
			\draw[dashed, very thick](6.8,-4.7)--(6.8,-3.7);
			\draw[dashed, very thick](6.8,-3.7)--(9.2,-3.7);
			\draw[dashed, very thick](8,0)--(8,-3.7);
			\draw[dashed, very thick](9.2,-3.7)--(9.2,-2.7);
			\draw[dashed, very thick](9.2,-2.7)--(11,-2.7);
			\draw[dashed, very thick](11,-2.7)--(11,-1.7);
			\draw[dashed, very thick](11,-1.7)--(12,-1.7);
			\draw[very thick](4,-5.7)--(12,-1.7);
			\draw[very thick](12,-1.7)--(14,-1.7);
			\draw[dashed, very thick](5,0)--(5,-4.7);
			\draw[dashed, very thick](6.8,0)--(6.8,-3.7);
			\draw[dashed, very thick](9.2,0)--(9.2,-2.7);
			\draw[dashed, very thick](11,0)--(11,-1.7);
			\draw[dashed, very thick](2.02,-6.7)--(14,-6.7);
			\filldraw [red] (8,-2) circle [radius=2.5pt];
			\draw[very thick](2,-2)--(2.2,-2);
			\node at (1.75,-2){$z_s$};
			\node at (7.4,-0.5){$j_s$};
			\node at (8.6,-0.5){$j_s+1$};
			\draw[very thick, ->](7,-1.5)--(7.8,-1.5);
			\node at (7.4,-1.9){$a^{j_s}$};
			\draw[very thick, ->](7.8,-2.5)--(7,-2.5);
			\node at (7.4,-2.9){$b^{j_s}$};
			\draw[very thick, ->](8.2,-1.5)--(9,-1.5);
			\node at (8.6,-1.9){$a^{j_s+1}$};
			\draw[very thick, ->](9,-2.5)--(8.2,-2.5);
			\node at (8.6,-2.9){$b^{j_s+1}$};
			\draw[very thick, ->](3.8,-1)--(3,-1);
			\node at (3.4,-1.4){$b^{1}$};
			\draw[very thick, ->](12,-1)--(12.8,-1);
			\node at (12.4,-1.4){$a^{J}$};
			\node at (8,0.2){$x_{s}$};
			\node at (7,0.2){$x_{j_s-1}$};
			\node at (9,0.2){$x_{j_s+1}$};
			\node at (4.5,-4.7){$h_{j-1}$};
			\node at (6.5,-3.7){$h_{j}$};
			\node at (11.5,-2.7){$h_{j+1}$};
			\node at (11,0.2){$x_{j+1}$};
			\node at (6,-1.0){$c_w^{j-1}$};
			\node at (6,-1.7){$\rho_w^{j-1}$};
			\node at (6,-2.4){$\alpha_w^{j-1}$};
			\node at (10.2,-0.5){$c_w^{j+1}$};
			\node at (10.2,-1.2){$\rho_w^{j+1}$};
			\node at (10.2,-1.9){$\alpha_w^{j+1}$};
			\node at (6.1,-5.5){$c_b^{j-1},\rho_b^{j-1}$};
			\node at (6,-6.3){$\alpha_b^{j-1}$};
			\node at (8,-4.5){$c_b^j,\rho_b^j,\alpha_b^j$}; \node at (10.3,-3.5){$c_b^{j+1},\rho_b^{j+1}$};
			\node at (10.2,-4.3){$\alpha_b^{j+1}$};
			\node at (8,-7.4){$c_\infty,\quad \rho_\infty,\quad \alpha_\infty$};
		\end{tikzpicture}
		\caption{Stepwise coupled modes of a line source.}
		\label{Figure2}
	\end{figure}

	Eq.~\eqref{eq.3} is physically similar to the sound field of an infinitely long line sound source in an $x$-dependent and $z$-dependent environment, as shown in Fig.~\ref{Figure2}. For an $x$-dependent marine environment, the classic solution is to divide the environment into a sufficiently large number $J$ of narrow segments along the $x$-direction \cite{Evans1983}, as shown in Fig.~\ref{Figure2}. In each flat segment, the marine environment is considered $x$-independent, so Eq.~\eqref{eq.3} can be transformed into:
	\begin{equation}
		\label{eq.5}
		\frac{\partial^2 \tilde{p}}{\partial x^2} + \rho(z) \frac{\partial}{\partial z}\left[\frac{1}{\rho(z)} \frac{\partial \tilde{p}}{\partial z}\right]+ \left[k^{2}(z)-k_{y}^{2}\right] \tilde{p}=-\delta\left(x-x_{s}\right) \delta\left(z-z_{s}\right)
	\end{equation}
	Segmenting the environment means that we solve the sound field in $J$ segments. Then, after obtaining the sound field for each of the $J$ segments, the complete sound field can be obtained by simple splicing.

	\subsection{Normal modes and coupled modes}
	Next, we derive the solution for the sound field of each segment; that is, we solve Eq.~\eqref{eq.5}. Using the technique of the separation of variables, the acoustic pressure can be decomposed into:
	\begin{equation}
		\label{eq.6}
		\tilde{p}(x,z)=X(x) \Psi(z)
	\end{equation}
	$\Psi (z)$ in Eq.~\eqref{eq.6} satisfies the following modal equation:
	\begin{subequations}
		\label{eq.7}
		\begin{gather}
			\label{eq.7a}
			\rho(z) \frac{\mathrm{d}}{\mathrm{d} z}\left[\frac{1}{\rho(z)} \frac{\mathrm{d} \Psi(z)}{\mathrm{d} z}\right]+\left[k^{2}(z)-k_y^2-\mu^{2}\right]\Psi(z)=0\\
			k(z)=(1+\mathrm{i}\eta\alpha)\omega/c(z),\quad \eta=(40\pi \log_{10}{\mathrm{e}})^{-1}
		\end{gather}
	\end{subequations}
	We redefine $k(z)$ in the above equation to account for the effect of attenuation on acoustic propagation. $\mu^2$, the constant obtained by the separation of variables, physically represents the horizontal wavenumber, and $\alpha$ is the attenuation coefficient in units of dB$/\lambda$, where $\lambda$ is the wavelength.

	Ignoring the contribution of the continuous spectrum and adding appropriate boundary conditions, Eq.~\eqref{eq.7} has a set of eigensolutions $\{(\mu_{m},\Psi_m)\}_{m=1}^\infty$, where $\Psi_m$ is also called the eigenmode. The eigenmodes of Eq.~\eqref{eq.7} should be normalized as follows:
	\begin{equation}
		\label{eq.8}
		\int_{0}^{H} \frac{\Psi_m(z) \Psi_n(z)}{\rho(z)}\mathrm {d} z=\delta_{m,n}
		\quad m,n = 1, 2, \dots
	\end{equation}
	where $H$ is the depth of the ocean and $\delta$ is the Kronecker delta function.

	$X(x)$ in Eq.~\eqref{eq.6} is related only to the horizontal range $x$ and satisfies:
	\begin{equation}
		\label{eq.9}
		\frac{\mathrm{d}^2 X(x)}{\mathrm{d} x^2}+\mu^{2} X(x)=-\frac{\delta(x-x_s) \Psi\left(z_{s}\right)}{\rho\left(z_{s}\right)}
	\end{equation}
	According to the theory of ordinary differential equations, the solution is:\begin{equation}
		\label{eq.10}
		X(x)=\frac{\mathrm{i}}{2 \rho\left(z_{s}\right)} \Psi\left(z_{s}\right) \frac{\mathrm{e}^{\mathrm{i}\mu\vert x\vert}}{\mu}
	\end{equation}
	Therefore, the fundamental solution of Eq.~\eqref{eq.5} is:
	\begin{equation}
		\label{eq.11}
		\tilde{p}(x,z)=\frac{\mathrm{i}}{2 \rho\left(z_s\right)} \sum_{m=1}^{\infty} \Psi_{m}\left(z_s\right) \Psi_{m}(z) \frac{\mathrm{e}^{\mathrm{i} \mu_m\vert x\vert}}{\mu_m}
	\end{equation}

	Eq.~\eqref{eq.11} is only the solution of the sound field of the $x$-independent segment, and for the entire two-dimensional waveguide, the normal modes of the segments are all coupled due to the dependence of the sound field on $x$. To account for the influence of coupling, referring to Eq.~\eqref{eq.11}, we define the sound field of the $j$-th segment as follows:
	\begin{subequations}
		\label{eq.12}
		\begin{gather}
			\tilde{p}^{j}(x, z)\approx \sum_{m=1}^{M}\left[a_{m}^{j} E_{m}^{j}(x)+b_{m}^{j} F_{m}^{j}(x)\right] \Psi_{m}^{j}(z) \\
			E_{m}^{j}(x)=\exp \left[\mathrm{i} \mu_ m^{j}\left(x-x_{j-1}\right)\right]\\
			F_{m}^{j}(x)=\exp \left[-\mathrm{i} \mu_ m^{j}\left(x-x_{j}\right)\right]
		\end{gather}
	\end{subequations}
	where $j=1,2,\cdots,J$; $M$ is the total number of normal modes needed to synthesize the sound field; $\{a_{m}^{j}\}_{m=1}^M$ and $\{b_{m}^{j}\}_{m=1}^M$ are the amplitudes of the forward and backward propagation modes, respectively, which are also called the coupling coefficients of the $j$-th segment; $\{\mu_m^j,\Psi_{m}^{j}(z)\}$ is the $m$-th eigensolution of the $j$-th segment; and $E_{m}^{j}(x)$ and $F_{m}^{j}(x)$ are the normalized range functions. For special cases when $j=1$, $x_{j-1}=x_1$.

	The method of coupling segments explicitly imposes two continuity conditions on the sides of the segments. The first segment condition requires that the acoustic pressure be continuous at the $j$-th side:
	\begin{equation}
		\label{eq.13}
		\tilde{p}^{j+1}\left(x_{j}, z\right)=\tilde{p}^{j}\left(x_{j}, z\right)
	\end{equation}
	Likewise, the second segment condition requires that the radial velocity of the acoustic pressure be continuous at the $j$-th side:
	\begin{equation}
		\label{eq.14}
		\frac{1}{\rho_{j+1}(z)} \frac{\partial \tilde{p}^{j+1}\left(x_{j}, z\right)}{\partial x}=\frac{1}{\rho_{j}(z)} \frac{\partial \tilde{p}^{j}\left(x_{j}, z\right)}{\partial x}
	\end{equation}
	For the first segment condition, we have:
	\begin{equation}
		\label{eq.15}
		\begin{aligned}
			&\sum_{m=1}^{M}\left[a_{m}^{j+1}E_{m}^{j+1}(x_j)+b_{m}^{j+1}F_{m}^{j+1}(x_j) \right] \Psi_{m}^{j+1}(z)=\\
			&\sum_{m=1}^{M}\left[a_{m}^{j} E_{m}^{j}\left(x_{j}\right)+b_{m}^{j}F_{m}^{j}\left(x_{j}\right)\right] \Psi_{m}^{j}(z)
		\end{aligned}
	\end{equation}
	Next, we apply the following operator to the above equation:
	\begin{equation*}
		\int_0^H (\cdot) \frac{\Psi_{\ell}^{j+1}(z)}{\rho_{j+1}(z)} \mathrm{d} z
	\end{equation*}
	Furthermore, we apply the orthogonal normalization formula in Eq.~\eqref{eq.8} to the eigenmodes of the $(j+1)$-th segment:
	\begin{subequations}
		\label{eq.16}
		\begin{gather}
			a_{\ell}^{j+1}+b_{\ell}^{j+1}F_{\ell}^{j+1}\left(x_{j}\right) =\sum_{m=1}^{M}\left[a_{m}^{j} E_{m}^{j}\left(x_{j}\right)+b_{m}^{j}\right] \tilde{c}_{\ell m}\\
			\tilde{c}_{\ell m}=\int_0^H \frac{\Psi_{\ell}^{j+1}(z) \Psi_{m}^{j}(z)}{\rho_{j+1}(z)} \mathrm{d} z, \quad \ell=1, \ldots, M
		\end{gather}
	\end{subequations}
	The above formula can be easily written in the following matrix--vector form:
	\begin{equation}
		\label{eq.17}
		\mathbf{a}^{j+1}+\mathbf{F}^{j+1}\mathbf{b}^{j+1}=\widetilde{\mathbf{C}}^{j}\left(\mathbf{E}^{j} \mathbf{a}^{j}+\mathbf{b}^{j}\right)
	\end{equation}
	Next, we find the partial derivatives on both sides of Eq.~\eqref{eq.12}:
	\begin{equation}
		\label{eq.18}
		\frac{1}{\rho_{j}} \frac{\partial \tilde{p}^{j}(x, z)}{\partial x} \simeq \frac{1}{\rho_{j}} \sum_{m=1}^{M} \mu_m^j\left[a_{m}^{j} E_{m}^{j}(x)-b_{m}^{j} F_{m}^{j}(x)\right] \Psi_{m}^{j}(z)
	\end{equation}
	For the second segment condition, we have:
	\begin{equation}
		\label{eq.19}
		\begin{aligned}
			\frac{1}{\rho_{j+1}}&\sum_{m=1}^{M}\mu_m^{j+1}\left[a_{m}^{j+1} -b_{m}^{j+1}F_{m}^{j+1}\left(x_{j}\right)\right]\Psi_{m}^{j+1}(z)=\\
			\frac{1}{\rho_{j}}&\sum_{m=1}^{M}\mu_m^j\left[a_{m}^{j} E_{m}^{j}\left(x_{j}\right)-b_{m}^{j}\right] \Psi_{m}^{j}(z)
		\end{aligned}
	\end{equation}
	Similarly, we apply the following operator to the above equation:
	\begin{equation*}
		\int_0^H(\cdot) \Psi_{\ell}^{j+1}(z) \mathrm{d} z
	\end{equation*}
	Next, we apply the orthogonal normalization formula in Eq.~\eqref{eq.8} to the eigenmodes of the $(j+1)$-th segment, yielding the following:
	\begin{subequations}
		\label{eq.20}
		\begin{gather}
			a_{\ell}^{j+1}-b_{\ell}^{j+1}F_{\ell}^{j+1} =\sum_{m=1}^{M}\left[a_{m}^{j} E_{m}^{j}\left(x_{j}\right)-b_{m}^{j}\right] \hat{c}_{\ell m}\\
			\hat{c}_{\ell m}=\frac{\mu_m^{j}}{\mu_\ell^{j+1}} \int_0^H \frac{\Psi_{\ell}^{j+1}(z) \Psi_{m}^{j}(z)}{\rho_{j}(z)} \mathrm{d} z, \quad \ell=1, \ldots, M
		\end{gather}
	\end{subequations}
	Then, the above formula can be naturally written in the following matrix--vector form:
	\begin{equation}
		\label{eq.21}
		\mathbf{a}^{j+1}-\mathbf{F}^{j+1}\mathbf{b}^{j+1}=\widehat{\mathbf{C}}^{j}\left(\mathbf{E}^{j} \mathbf{a}^{j}-\mathbf{b}^{j}\right)
	\end{equation}
	Eqs.~\eqref{eq.17} and \eqref{eq.21} can be combined into the following form:
	\begin{subequations}
		\label{eq.22}
		\begin{align}
			\left[\begin{array}{l}
				\mathbf{a}^{j+1} \\
				\mathbf{b}^{j+1}
			\end{array}\right]&=\left[\begin{array}{ll}
				\mathbf{R}_{1}^{j} & \mathbf{R}_{2}^{j} \\
				\mathbf{R}_{3}^{j} & \mathbf{R}_{4}^{j}
			\end{array}\right]\left[\begin{array}{c}
				\mathbf{a}^{j} \\
				\mathbf{b}^{j}
			\end{array}\right]\\
			\mathbf{R}_{1}^{j} &=\frac{1}{2}\left(\widetilde{\mathbf{C}}^{j}+\widehat{\mathbf{C}}^{j}\right) \mathbf{E}^{j} \\
			\mathbf{R}_{2}^{j} &=\frac{1}{2}\left(\widetilde{\mathbf{C}}^{j}-\widehat{\mathbf{C}}^{j}\right) \\
			\mathbf{R}_{3}^{j} &=\frac{1}{2}\left(\mathbf{F}^{j+1}\right)^{-1}\left(\widetilde{\mathbf{C}}^{j}-\widehat{\mathbf{C}}^{j}\right) \mathbf{E}^{j} \\
			\mathbf{R}_{4}^{j} &=\frac{1}{2}\left(\mathbf{F}^{j+1}\right)^{-1}\left(\widetilde{\mathbf{C}}^{j}+\widehat{\mathbf{C}}^{j}\right)
		\end{align}
	\end{subequations}

	Note that since the sound source is located at $(x_s, z_s)$, a virtual boundary $j_s$ needs to be introduced at the horizontal distance $x_s$. All boundaries except $j_s$ still satisfy the boundary conditions of Eqs.~\eqref{eq.13} and \eqref{eq.14}. However, due to the existence of the sound source, the boundary conditions at the virtual boundary $j_s$ need to be modified. Applying segment condition \eqref{eq.13} to this virtual boundary, we can obtain:
	\begin{subequations}
		\label{eq.23}
		\begin{gather}
			a_{m}^{j_{s}+1}+b_{m}^{j_{s}+1} F_{m}^{j_{s}+1}=a_{m}^{j_{s}} E_{m}^{j_{s}}+b_{m}^{j_{s}}\\
			\label{eq.23b}
			\mathbf{a}^{j_{s}+1}+\mathbf{F}^{j_{s}+1} \mathbf{b}^{j_{s}+1}=\mathbf{E}^{j_{s}} \mathbf{a}^{j_{s}}+\mathbf{b}^{j_{s}}
		\end{gather}
	\end{subequations}
	This equation is consistent with the form of Eq.~\eqref{eq.17}, but the acoustic properties on both sides of the boundary $j_s$ are completely consistent, resulting in $\widetilde{\mathbf{C}}=\mathbf{I}$.

	The integral of Eq.~\eqref{eq.11} yields the following formula:
	\begin{equation}
		\label{eq.24}
		\frac{\mathrm{d} X}{\mathrm{d} x}\bigg\vert_{x_{s}^{-}} ^{x_{s}^{+}}=-\frac{\Psi_{m}^{j_{s}}\left(z_{s}\right)}{\rho\left(z_{s}\right)}
	\end{equation}
	Considering Eq.~\eqref{eq.10} and $\mu_{m}^{j_s}=\mu_{m}^{j_s+1}$, we can obtain:
	\begin{subequations}
		\label{eq.25}
		\begin{gather}
			\left(\mathrm{i} \mu_{m}^{j_{s}+1} a_{m}^{j_{s}+1}-\mathrm{i} \mu_{m}^{j_{s}+1} b_{m}^{j_{s}+1} F_{m}^{j_{s}+1}\right)-\left(\mathrm{i} \mu_{m}^{j_{s}} a_{m}^{j_{s}} E_{m}^{j_{s}}-\mathrm{i} \mu_{m}^{j_{s}} b_{m}^{j_{s}}\right)=-\frac{\Psi_{m}^{j_{s}}\left(z_{s}\right)}{\rho\left(z_{s}\right)}\\
			\label{eq.25b}
			\mathbf{a}^{j_{s}+1}-\mathbf{F}^{j_{s}+1} \mathbf{b}^{j_{s}+1}=\mathbf{E}^{j_{s}} \mathbf{a}^{j_{s}}-\mathbf{b}^{j_{s}}-\mathbf{s}\\
			s_{m}=-\frac{\mathrm{i}}{\rho\left(z_{s}\right)} \frac{\Psi_{m}^{j_{s}}\left(z_{s}\right)}{\mu_{m}^{j_{s}}}, \quad m=1,2, \cdots, M
		\end{gather}
	\end{subequations}
	Eq.~\eqref{eq.25b} also corresponds to Eq.~\eqref{eq.21}. Because the acoustic properties on both sides of the virtual boundary $j_s$ are exactly the same, $\widehat{\mathbf{C}}=\mathbf{I}$. Moreover, Eqs.~\eqref{eq.23b} and \eqref{eq.25b} demonstrate that the relationship between the coupling coefficients of the $j_s$-th and $(j_s+1)$-th segments is:
	\begin{equation}
		\label{eq:26}
		\left[\begin{array}{l}
			\mathbf{a}^{j_s+1} \\
			\mathbf{b}^{j_s+1}
		\end{array}\right]=\left[\begin{array}{ll}
			\mathbf{E}^{j_s} & \mathbf{0} \\
			\mathbf{0} & (\mathbf{F}^{j_s+1})^{-1}
		\end{array}\right]\left[\begin{array}{c}
			\mathbf{a}^{j_s} \\
			\mathbf{b}^{j_s}
		\end{array}\right]+\left[\begin{array}{c}
			-\frac{1}{2}\mathbf{s} \\
			\frac{1}{2}\left(\mathbf{F}^{j_s+1}\right)^{-1}\mathbf{s}
		\end{array}\right]
	\end{equation}

	Since the sound source is located at $x_s \neq 0$, the segment conditions at $x=0$ and $x=\infty$ are the radiation condition, i.e., $\mathbf{a}^1=0$ and $\mathbf{b}^J=0$. The global matrix used to determine the coupling coefficients can be assembled as follows:
	\begin{tiny}
		\begin{equation}
			\label{eq.27}
			\left[\begin{array}{ccccccccccccc}
				\mathbf{R}_{2}^{1} & -\mathbf{I} & \mathbf{0} & & & & & \\
				\mathbf{R}_{4}^{1} & \mathbf{0} & -\mathbf{I} & & & & & \\
				& \mathbf{R}_{1}^{2} &\mathbf{R}_{2}^{2} &-\mathbf{I} &\mathbf{0} & & & & \\
				& \mathbf{R}_{3}^{2} &\mathbf{R}_{4}^{2} &\mathbf{0}  &-\mathbf{I} & & & & \\
				& & &\ddots & \ddots & \ddots & \ddots& & \\
				& & & & & \mathbf{E}^{j_s} & \mathbf{0} & -\mathbf{I} & \mathbf{0} \\
				& & & & & \mathbf{0} & (\mathbf{F}^{j_s+1})^{-1} & \mathbf{0} & -\mathbf{I}\\
				& & & & & & &  \ddots & \ddots & \ddots & \ddots&  \\
				& & & & & & & & &\mathbf{R}_{1}^{J-1} & \mathbf{R}_{2}^{J-1} & -\mathbf{I} \\
				& & & & & & & & &\mathbf{R}_{3}^{J-1} & \mathbf{R}_{4}^{J-1} & \mathbf{0}
			\end{array}\right]\left[\begin{array}{c}
				\mathbf{b}^{1} \\
				\mathbf{a}^{2} \\
				\mathbf{b}^{2} \\
				\mathbf{a}^{3} \\
				\vdots \\
				\mathbf{a}^{j_s} \\
				\mathbf{b}^{j_s} \\
				\vdots \\
				\mathbf{b}^{J-1}\\
				\mathbf{a}^{J}
			\end{array}\right]=\left[\begin{array}{c}
				\mathbf{0} \\
				\mathbf{0} \\
				\mathbf{0} \\
				\mathbf{0} \\
				\vdots \\
				\frac{1}{2}\mathbf{s} \\
				-\frac{1}{2}(\mathbf{F}^{j_s+1})^{-1}\mathbf{s} \\
				\vdots \\
				\mathbf{0}\\
				\mathbf{0} \\
			\end{array}\right]
		\end{equation}
	\end{tiny}
	Note that the global matrix is a sparse band matrix of order $(2J-2) \times M$; the upper bandwidth is $(2M-1)$, and the lower bandwidth is $(3M-1)$. The linear algebraic equations are then solved to obtain the coupling coefficients, and Eq.~\eqref{eq.12} is implemented to synthesize the sound pressure field of the entire waveguide.

	Specifically, when $x_s=0$, the segment condition at the acoustic source $x_s=0$ is:
	\begin{subequations}
		\label{eq.28}
		\begin{gather}
			\mathbf{a}^{1}=\mathbf{s}\\
			s_{m}=\frac{\mathrm{i}}{2 \rho\left(z_{s}\right)} \Psi_{m}^{1}\left(z_{s}\right) \frac{\mathrm{e}^{\mathrm{i} \mu_{m}^{1} x_{1}}}{\mu_{m}^{1}}, \quad m=1, \cdots, M
		\end{gather}
	\end{subequations}
	Therefore, the global matrix becomes:
	\begin{footnotesize}
		\begin{equation}
			\label{eq.29}
			\left[\begin{array}{ccccccccc}
				\mathbf{I} & \mathbf{0} & \mathbf{0} & & & & & \\
				\mathbf{R}_{1}^{1} & \mathbf{R}_{2}^{1} & -\mathbf{I} & \mathbf{0} & & & & \\
				\mathbf{R}_{3}^{1} & \mathbf{R}_{4}^{1} & \mathbf{0} & -\mathbf{I} & & & & \\
				& & \ddots & \ddots & \ddots & \ddots & & \\
				& & & &\mathbf{R}_{1}^{J-2} & \mathbf{R}_{2}^{J-2} & -\mathbf{I} & \mathbf{0} & \\
				& & & &\mathbf{R}_{3}^{J-2} & \mathbf{R}_{4}^{J-2} & \mathbf{0} & -\mathbf{I} & \\
				& & & & & &\mathbf{R}_{1}^{J-1} & \mathbf{R}_{2}^{J-1} & -\mathbf{I} \\
				& & & & & &\mathbf{R}_{3}^{J-1} & \mathbf{R}_{4}^{J-1} & \mathbf{0}
			\end{array}\right]\left[\begin{array}{c}
				\mathbf{a}^{1} \\
				\mathbf{b}^{1} \\
				\mathbf{a}^{2} \\
				\vdots \\
				\mathbf{b}^{J-2} \\
				\mathbf{a}^{J-1} \\
				\mathbf{b}^{J-1} \\
				\mathbf{a}^{J}
			\end{array}\right]=\left[\begin{array}{c}
				\mathbf{s} \\
				\mathbf{0} \\
				\mathbf{0} \\
				\vdots \\
				\mathbf{0} \\
				\mathbf{0} \\
				\mathbf{0} \\
				\mathbf{0}
			\end{array}\right]
		\end{equation}
	\end{footnotesize}

	\section{Numerical Discretization and Calculation}
	\subsection{Range-independent marine environment}
	In the previous section, we fully deduce the solution of a quasi-three-dimensional waveguide. Here, we solve for the eigensolutions in each segment, that is, the solutions $\{(\mu_{m}^j,\Psi_m^j)\}_{m=1}^M$ of Eq.~\eqref{eq.7}. For the sediment-covered marine environment considered in Figs.~\ref{Figure1} and \ref{Figure2}, within range-independent segments, $\rho(z)$, $c(z)$ and $\alpha(z)$ are discontinuous at the interface $z=h$; thus, the ocean is divided into a discontinuous water column layer and a bottom sediment layer. The marine environmental parameters are defined in the water column and bottom sediment as follows:
	\begin{subequations}
		\label{eq.30}
		\begin{gather}
			c(z)= \begin{cases}
				c_w(z),&0 \leq z \leq h\\
				c_b(z),&h \leq z \leq H\\
				c_\infty,&z\geq H
			\end{cases}\\
			\rho(z)= \begin{cases}
				\rho_w(z),&0 \leq z \leq h\\
				\rho_b(z),&h \leq z \leq H\\
				\rho_\infty,&z\geq H
			\end{cases}\\
			\alpha(z)= \begin{cases}
				\alpha_w(z),&0 \leq z \leq h\\
				\alpha_b(z),&h \leq z \leq H\\
				\alpha_\infty,&z\geq H
			\end{cases}
		\end{gather}
	\end{subequations}

	In this case, to solve Eq.~\eqref{eq.7}, it is necessary to impose boundary conditions at the sea surface ($z=0$) and the seabed ($z=H$) and interface conditions at the interface ($z=h$). The sea surface is usually set as a pressure-release boundary:
	\begin{equation}
		\label{eq.31}
		\Psi(z=0)=0
	\end{equation}
	In contrast, the seabed can be either a pressure-release boundary or a rigid seabed:
	\begin{subequations}
		\label{eq.32}
		\begin{gather}
			\Psi(z=H)=0\\
			\Psi'(z=H)=0
		\end{gather}
	\end{subequations}
	Furthermore, for an acoustic half-space, the lower boundary should generally satisfy the following \cite{Jensen2011}:
	\begin{equation}
		\label{eq.33}
		\Psi(H)+\frac{\rho_\infty}{\rho_b(H)\gamma_\infty} \Psi'(H)=0,\quad \gamma_{\infty}=\sqrt{\mu^{2}-\frac{\omega^{2}}{c_{\infty}^{2}}}
	\end{equation}
	Note that the eigenfunction value of the $z\in[H,+\infty]$ interval should be added to the

	orthogonal normalization
	at this time:
	\begin{equation}
		\label{eq.34}
		\int_{0}^{H} \frac{\Psi_{m}^{2}(z)}{\rho(z)} \mathrm{d} z +\frac{\Psi_{m}^{2}(H)}{2 \rho_\infty\gamma_{\infty}} =1,\quad m=1,2,\dots
	\end{equation}
	The interface conditions are defined as follows:
	\begin{subequations}
		\label{eq.35}
		\begin{gather}
			\Psi(h^{-})=\Psi(h^{+}) \\
			\frac{1}{\rho(z=h^{-})} \frac{\mathrm{d}\Psi(h^{-} )}{\mathrm {d}z}= \frac{1}{\rho(z=h^{+})}\frac{\mathrm{d}\Psi(h^{+} )}{\mathrm {d}z}
		\end{gather}
	\end{subequations}
	where the superscripts $h^{-}$ and $h^{+}$ denote above and below $h$, respectively.

	\subsection{Chebyshev--Tau spectral method}
	\label{section3.2}
	Traditionally, two numerical methods are used to solve for the local normal modes: the FDM, which is represented by the KRAKEN program \cite{Kraken2001}, and the Galerkin method, which is represented by COUPLE \cite{Couple}. The spectral method is a new type of numerical method known for its high precision that has been introduced into computational ocean acoustics in recent years. Generally, the spectral method includes both the Galerkin method (with orthogonal polynomials as the basis functions) and the collocation method (with the extreme points of orthogonal polynomials as the nodes). Numerical experiments have shown that the former method is slightly more accurate than the latter \cite{Jshen2011}. This paper employs an exceptional type of spectral method: the Chebyshev--Tau spectral method. Different from the classic Galerkin-type spectral method, the Chebyshev--Tau spectral method uses the original Chebyshev polynomial as the basis function and does not require the basis function to satisfy the boundary conditions but instead transforms the boundary conditions into the spectral space to enforce the requirement \cite{Lanczos1938}.

	Here, we employ the Chebyshev--Tau spectral method to solve for the local modes (Eq.~\eqref{eq.7a}) in the $J$ range-independent segments illustrated in Fig.~\ref{Figure3}. Since the Chebyshev polynomial $\{T_i(x)\}$, that is, the basis function, is defined in $x\in[-1,1]$, the equation to be solved, Eq.~\eqref{eq.7a}, must first be scaled to $x\in[-1,1]$ as:
	\begin{equation}
		\label{eq.36}
		\frac{4}{\vert \Delta h\vert^2}\rho(x)\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{\rho(x)}\frac{\mathrm{d}\Psi(x)}{\mathrm {d}x}\right) +\left[k^{2}(z)-k_y^2\right]\Psi(x) = \mu^2 \Psi(x)
	\end{equation}
	where $\Delta h$ denotes the length of the domain. The variable to be determined, $\Psi(x)$, is transformed into the spectral space spanned by the basis functions $\{T_i(x)\}_{i=0}^N$:
	\begin{equation}
		\label{eq.37}
		\Psi(x) \approx \sum_{i=0}^{N}\hat{\Psi}_{i}T_i(x)
	\end{equation}
	where $\{\hat{\Psi}_i\}_{i=0}^N$ represents the spectral coefficients of $\Psi(x)$. In a sense, this is a function approximation that obviously becomes increasingly accurate as $N$ increases. Due to the advantageous properties of the Chebyshev polynomial, the following is easy to prove \cite{Canuto2006}:
	\begin{subequations}
		\label{eq.38}
		\begin{gather}
			\label{eq.38a}
			\hat{\Psi}'_i \approx \frac{2}{c_i}
			\sum_{\substack{j=i+1,\\
					j+i=\text{odd}
			}}^{N} j \hat{\Psi}_j, \quad c_0=2,c_{i>1}=1
			\Longleftrightarrow \bm{\hat{\Psi}}' \approx \mathbf{D}_N \bm{\hat{\Psi}} \\
			\label{eq.38b}
			\widehat{(v\Psi)}_i \approx
			\frac{1}{2} \sum_{m+n=i}^{N} \hat{\Psi}_m\hat{v}_n +
			\frac{1}{2} \sum_{\vert m-n\vert=i}^{N} \hat{\Psi}_m\hat{v}_n  \Longleftrightarrow  \widehat{\bm{(v\Psi)}} \approx \mathbf{C}_v \bm{\hat{\Psi}}
		\end{gather}
	\end{subequations}
	Eq.~\eqref{eq.38a} describes the relationship between the spectral coefficients of a function and the spectral coefficients of its derivative. Likewise, Eq.~\eqref{eq.38b} represents the relationship between the spectral coefficients of the product of two functions and the spectral coefficients of the individual functions, where the right-hand side is the matrix--vector representation of the relationship.

	As shown in Eqs.~\eqref{eq.37} and \eqref{eq.38}, Eq.~\eqref{eq.7a} can be discretized into the following matrix--vector form:
	\begin{equation}
		\label{eq.39}
		\left(\frac{4}{\vert\Delta h\vert^2}\mathbf{C}_{\rho}\mathbf{D}_{N}\mathbf{C}_{1/\rho}\mathbf{D}_{N}+\mathbf{C}_{k^2}-k_y^2\mathbf{I}\right)\bm{\hat{\Psi}}
		= \mu^2 \bm{\hat{\Psi}}
	\end{equation}
	where $\mathbf{I}$ is the identity matrix. The above equation is obviously a matrix eigenvalue problem, but the boundary conditions are not considered at this time.

	For the waveguide in Eq.~\eqref{eq.30},

	the modal equations (Eqs.~\eqref{eq.7a} and \eqref{eq.39})
	must be established in both the water column and the bottom sediment. In the flat segments, a single set of basis functions cannot span both layers since the Chebyshev polynomial corresponding to the interface $h$ is not continuously differentiable. Thus, we apply the domain decomposition strategy \cite{Min2005} to Eq.~\eqref{eq.7a} and split the domain interval into two subintervals:
	\begin{equation}
		\label{eq.40}
		\Psi(z)= \begin{cases}
			\Psi_w(z) = \Psi_w(x) \approx\sum_{i=0}^{N_w}\hat{\Psi}_{w,i}T_i(x_w),\quad x_w=-\frac{2z}{h}+1,&
			0\leq z\leq h \\
			\Psi_b(z) = \Psi_b(x) \approx\sum_{i=0}^{N_b}\hat{\Psi}_{b,i}T_i(x_b),\quad
			x_b=-\frac{2z}{H-h}+\frac{H+h}{H-h},&
			h\leq z\leq H
		\end{cases}
	\end{equation}
	where $N_w$ and $N_b$ are the spectral truncation orders in the water column and bottom sediment, respectively, and $\{\hat{\Psi}_{w,i}\}_{i=0}^{N_w}$ and $\{\hat{\Psi}_{b,i}\}_{i=0}^{N_b}$ are the spectral coefficients in these two layers. Similar to Eq.~\eqref{eq.39}, the modal equations in the water column and bottom sediment can be discretized into matrix--vector form:
	\begin{subequations}
		\begin{align}
			\label{eq.41a}
			&\mathbf{A}
			\bm{\hat{\Psi}}_w
			= k_x^2
			\bm{\hat{\Psi}}_w,
			&\mathbf{A}
			=\frac{4}{h^2}\mathbf{C}_{\rho_w}\mathbf{D}_{N_w}\mathbf{C}_{1/\rho_w}\mathbf{D}_{N_w}+\mathbf{C}_{k_w^2}-k_y^2\mathbf{I}\\
			\label{eq.41b}
			&\mathbf{B}
			\bm{\hat{\Psi}}_b
			= k_x^2
			\bm{\hat{\Psi}}_b,&\mathbf{B}=
			\frac{4}{(H-h)^2}\mathbf{C}_{\rho_b}\mathbf{D}_{N_b}\mathbf{C}_{1/\rho_b}\mathbf{D}_{N_b}+\mathbf{C}_{k_b^2}-k_y^2\mathbf{I}
		\end{align}
	\end{subequations}
	where $\mathbf{A}$ and $\mathbf{B}$ are square matrices of order $(N_w+1)$ and $(N_b+1)$, respectively, and $\bm{\hat{\Psi}}_w$ and $\bm{\hat{\Psi}}_b$ are column vectors composed of $\{\hat{\Psi}_{w,i}\}_{i=0}^{N_w}$ and $\{\hat{\Psi}_{b,i}\}_{i=0}^{N_b}$, respectively. Since the interface conditions are related to both the water column and the bottom sediment, Eqs.~\eqref{eq.41a} and \eqref{eq.41b} should be solved simultaneously as follows:
	\begin{equation}
		\label{eq.42}
		\left[\begin{array}{cc}
			\mathbf{A}&\mathbf{0}\\
			\mathbf{0}&\mathbf{B}\\
		\end{array}\right]
		\left[\begin{array}{c}
			\bm{\hat{\Psi}}_w\\
			\bm{\hat{\Psi}}_b\\
		\end{array}
		\right]=\mu^2\left[
		\begin{array}{c}
			\bm{\hat{\Psi}}_w\\
			\bm{\hat{\Psi}}_b\\
		\end{array}\right]
	\end{equation}
	The boundary conditions and interface conditions in Eqs.~\eqref{eq.31}--\eqref{eq.35} must also be expanded to the spectral space and explicitly added to Eq.~\eqref{eq.42}. For details regarding the treatment of the boundary conditions and interface conditions, please see Eq.~(36) in Ref.~\cite{Tuhw2020b}. Rearranging and combining the modified rows by an elementary row transformation, Eq.~\eqref{eq.42} can be rewritten into the form of the following block matrix:
	\begin{equation}
		\label{eq.43}
		\left[
		\begin{array}{cc}
			\mathbf{L}_{11}&\mathbf{L}_{12}\\
			\mathbf{L}_{21}&\mathbf{L}_{22}\\
		\end{array}
		\right]\left[
		\begin{array}{c}
			\bm{\hat{\Psi}}_1\\
			\bm{\hat{\Psi}}_2
		\end{array}
		\right]=\mu^2\left[
		\begin{array}{c}
			\bm{\hat{\Psi}}_1\\
			\mathbf{0}
		\end{array}\right]
	\end{equation}
	where $\mathbf{L}_{11}$ is a square matrix of order $(N_w+N_b-2)$ and $\mathbf{L}_{22}$ is a square matrix of order 4. Solving this mixed linear algebraic system can yield the wavenumbers and spectral coefficients of the eigenmodes $(\mu,\bm{\hat{\Psi}}_w^j,\bm{\hat{\Psi}}_b^j)$.

	For the acoustic half-space boundary condition in Eq.~\eqref{eq.33}, since $\gamma_\infty$ contains the eigenvalue to be solved, $\mu$, the elements of $\mathbf{L}$ on the left side of Eq.~\eqref{eq.43} contain $\mu$; thus, Eq.~\eqref{eq.43} is no longer a linear algebraic eigenvalue problem and can be solved iteratively only by a root-finding algorithm \cite{Sabatini2019}. Since a prior estimate of $\mu$ is usually not available, many of the existing numerical programs following similar principles fail to converge to a specific root in some cases. To avoid the same dilemma, we use $k_{z,\infty}=\sqrt{k_\infty^{2}-\mu^{2}}$ to transform the modal equation and Eq.~\eqref{eq.33} as follows:
	\begin{subequations}
		\label{eq.44}
		\begin{gather}
			\label{eq.44a}
			\rho(z) \frac{\mathrm{d}}{\mathrm{d} z}\left(\frac{1}{\rho(z)} \frac{\mathrm{d} \Psi}{\mathrm{d}z}\right)+\left[k^{2}(z)-k_y^2-k_{\infty}^{2}+k_{z,\infty}^{2}\right] \Psi=0 \\
			\label{eq.44b}
			\frac{\mathrm{i}\rho_\infty}{\rho_b(H)}\frac{\mathrm{d} \Psi(z)}{\mathrm{d} z}\bigg\vert_{z=H}+k_{z,\infty} \Psi(H)=0
		\end{gather}
	\end{subequations}
	where $k_\infty=(1+\mathrm{i}\eta\alpha_\infty)\omega/c_\infty$ is a complex constant. In the Chebyshev--Tau spectral method, Eq.~\eqref{eq.44a} can be discretized into:
	\begin{equation}
		\label{eq.45}
		\left[\mathbf{U}+k_{z,\infty}^{2} \mathbf{I}\right] \bm{\hat{\Psi}}=\mathbf{0},\quad \mathbf{U}=\mathbf{L}-k_{\infty}^{2}\mathbf{I}
	\end{equation}
	After adding the boundary condition in Eq.~\eqref{eq.44b} related to $k_{z,\infty}$, Eq.~\eqref{eq.45} finally becomes the following polynomial eigenvalue problem:
	\begin{equation}
		\label{eq.46}
		\left[\mathbf{U}+k_{z,\infty} \mathbf{V}+k_{z,\infty}^{2} \mathbf{W}\right] \bm{\hat{\Psi}}=\mathbf{0}
	\end{equation}
	where $\mathbf{U}$ in Eq.~\eqref{eq.46} is not exactly identical to that in Eq.~\eqref{eq.45}, as it has been modified by boundary conditions and interface conditions; nevertheless, we denote it as $\mathbf{U}$. This polynomial eigenvalue problem can be efficiently solved by transforming into a general matrix eigenvalue problem, although the matrices are twice as large:
	\begin{subequations}
		\label{eq:47}
		\begin{gather}
			\tilde{\mathbf{U}} \tilde{\bm{\Psi}}=k_{z,\infty}^{2} \tilde{\mathbf{V}} \tilde{\bm{\Psi}}, \\
			\tilde{\mathbf{U}}=\left[\begin{array}{cc}
				-\mathbf{V} & -\mathbf{U} \\
				\mathbf{I} & 0
			\end{array}\right], \quad \tilde{\mathbf{V}}=\left[\begin{array}{cc}
				\mathbf{W} & 0 \\
				0 & \mathbf{I}
			\end{array}\right], \quad \tilde{\bm{\Psi}}=\left[\begin{array}{c}
				k_{z,\infty} \bm{\hat{\Psi}} \\
				\bm{\hat{\Psi}}
			\end{array}\right]
		\end{gather}
	\end{subequations}
	Obviously, $\mathbf{U}$, $\mathbf{V}$, and $\mathbf{W}$ are all very sparse.

	Using Eq.~\eqref{eq.37}, the eigenvectors $\bm{\hat{\Psi}}_w$ and $\bm{\hat{\Psi}}_b$ can be easily transformed into $\bm{\Psi}_w(z),z\in[0,h]$ and $\bm{\Psi}_b(z),z\in[h,H]$. The vectors $\bm{\Psi}_w$ and $\bm{\Psi}_b$ are stacked into a single column vector to form $\bm{\Psi}$. Then, either Eq.~\eqref{eq.8} or Eq.~\eqref{eq.34} is used to normalize $\bm{\Psi}$, and finally, a set of eigenpairs $\{\mu,\Psi(x)\}$ is obtained.

	\subsection{Numerical Sound Field Synthesis}

	According to the above analysis, after calculating $\tilde{p}(x,k_y,z)$, the spatial sound field can be obtained by the inverse Fourier transform in Eq.~\eqref{eq.2b}. However, $\tilde{p}(x,k_y,z)$ has singularities along the real $k_y$ axis; therefore, we refer to the approach in Ref.~\cite{Luowy2016} and convert the integral of Eq.~\eqref{eq.2b} into a contour integral on the complex plane. The contour, as shown in Fig.~\ref{Figure3}, takes the following form:
	\begin{equation}
		\label{eq.48}
		k_{y}(q)=q-\mathrm{i} \varepsilon \tanh (\delta q), \quad-\infty<q<\infty
	\end{equation}

	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale = 0.8,samples=600,domain=-6:6]
			\draw[very thick, ->](2,-5)--(14,-5) node[right]{Re($k_y$)};
			\draw[very thick, ->](8,-9)--(8,0) node[above]{{Im($k_y$)}};
			\filldraw [red] (8.2,-2) circle [radius=2.5pt];
			\filldraw [red] (8.7,-2.8) circle [radius=2.5pt];
			\filldraw [red] (9.4,-3.5) circle [radius=2.5pt];
			\filldraw [red] (10.2,-4.1) circle [radius=2.5pt];
			\filldraw [red] (11,-4.5) circle [radius=2.5pt];
			\filldraw [red] (12,-4.8) circle [radius=2.5pt];
			\filldraw [red] (12.9,-4.9) circle [radius=2.5pt];

			\draw[color=blue,very thick,smooth] plot(\x+8,{-3*((exp(\x))-(exp(-\x)))/((exp(\x))+(exp(-\x)))-5})
			node[below left] {$k_{y}(q)=q-\mathrm{i} \varepsilon \tanh (\delta q)$};

		\end{tikzpicture}
		\caption{Complex integration contour strategy.}
		\label{Figure3}
	\end{figure}

	The integral in Eq.~\eqref{eq.2b} is transformed into the following form:
	\begin{equation}
		\label{eq.49}
		p(x, y, z)=\frac{1}{\pi} \int_{0}^{\infty} \tilde{p}\left(x, k_{y}(q), z\right) \cos \left(y k_{y}(q)\right)\left[1-\mathrm{i} \varepsilon \delta \operatorname{sech}^{2}(\delta q)\right] \mathrm{d} q
	\end{equation}
	Since the interval is infinite, integrating over this interval is obviously impossible in an actual numerical calculation. Luo et al demonstrated that sufficiently accurate results can be obtained by setting the integration interval of $q$ in the above formula as $[0, 1.5k]$. $q$ takes equidistant, discrete values in the interval; thus, the above formula can be transformed into a general summation problem \cite{Luowy2015}.

	\section{Algorithm}
	Based on the above derivation, the proposed algorithm can be summarized as follows:

	\textbf{Input:} Data from the marine environment and the program parameters.

	\textbf{Output:} The three-dimensional sound pressure field.
	\begin{enumerate}
		\item
		Set the parameters.

		The parameters include the frequency $f$, the truncation interval $q$ of the Fourier integral (Eq.~\eqref{eq.48}), the number of sampling points $N_q$ in the $k_y$ wavenumber domain, the location of the source $(x_s,z_s)$, the total depth of the ocean $H$, the topography of the seabed, the number of acoustic profiles, and the specific information pertaining to each group of acoustic profiles. In addition, the parameters should include the spectral truncation orders ($N_w$ and $N_b$), the horizontal and vertical resolutions, and the number of modes $M$ to be coupled. If the lower boundary is the upper interface of an acoustic half-space, the speed $c_\infty$, density $\rho_\infty$ and attenuation $\alpha_\infty$ of the half-space must also be specified.

		\item
		Segment the marine environment based on the seabed topography, frequency of the source and acoustic profiles.

		Jensen confirmed that the strictest segmentation criterion is $\Delta r \le \lambda/4$, where $\lambda=\min\{c_w,c_b\}/f$ \cite{Jensen1998}. We assume that the $x$-dependence is divided into $J$ segments.

		\item
		Discretize the $q$ wavenumber domain into $N_q$ equidistant points, and perform the calculation from the fourth to the seventh step for each $k_y(q)$.

		\item
		Apply the Chebyshev--Tau spectral method to solve for the eigenpairs $\{(\mu_m^j,\Psi_m^j)\}_{m=1}^{M}$ of the $J$ $x$-independent segments (see Sec.~\ref{section3.2}).

		\item
		Calculate the coupling submatrices $\{\mathbf{R}_{1}^{j}\}_{j=1}^{J-1}$, $\{\mathbf{R}_{2}^{j}\}_{j=1}^{J-1}$, $\{\mathbf{R}_{3}^{j}\}_{j=1}^{J-1}$, and $\{\mathbf{R}_{4}^{j}\}_{j=1}^{J-1}$.

		\item
		Calculate $\mathbf{s}$, construct the global matrix according to either Eq.~\eqref{eq.27} or Eq.~\eqref{eq.29}, and solve the system of linear equations to obtain the coupling coefficients $(\{\mathbf{a}^j\}_{j=1}^J,\{\mathbf{b}^j\}_{j=1}^J)$ of the $J$ segments.

		\item
		Calculate the $k_y$--domain sound field $\tilde{p}(x,k_y(q),z)$.

		The sound pressure subfields of the segments are calculated according to Eq.~\eqref{eq.12}. The subfields of the $J$ segments are then spliced to form the sound pressure field of the entire $xoz$--plane. The fourth to the seventh steps are naturally parallel because the solution of $\tilde{p}(x,k_y(q),z)$ is independent for each $k_y(q)$ (see Sec.~\ref{section2.1}).

		\item
		Calculate the sound field $\tilde{p}(x,y,z)$.

		According to Eq.~\eqref{eq.50}, the $N_q$ discrete $\tilde{p}(x,k_y(q),z)$ fields are accumulated in the following manner to obtain the spatial sound field $\tilde{p}(x,y,z)$:
		\begin{equation}
			\label{eq.50}
			\begin{split}
				p(x, y, z)=\frac{3k}{2\pi(N_q-1)} \sum_{i=1}^{N_q} \tilde{p}\left(x, k_{y,i}, z\right) \cos \left(y k_{y,i}\right)\\
				\times\left\{1- \frac{\mathrm{i}}{4\pi\lg \mathrm{e}} \operatorname{sech}^{2}\left[\frac{(N_q-1)k_{y,i}}{9}\right]\right\}\\
			\end{split}
		\end{equation}
		The above formula is equivalent to taking $\frac{3k}{2(N_q-1)}$ for $\mathrm{d}q$, $\frac{3\mathrm{d}q}{2\pi\lg \mathrm{e}}$ for $\varepsilon$ and $\frac{1}{6\mathrm{d}q}$ for $\delta$ in Eq.~\eqref{eq.50} \cite{Luowy2016}.
	\end{enumerate}

	\section{Numerical Simulation}
	Here, we present a program named SPEC3D that we developed based on the above algorithm, and we verify the correctness of the algorithm through several numerical experiments. The transmission loss (TL) field is often used in actual displays to compare and analyze sound fields \cite{Jensen2011}. Therefore, to present the sound field results, the TL of the acoustic pressure is defined as $\text{TL}=-20\log_{10}(\vert p\vert/\vert p_0\vert)$ in dB, where $p_0=\exp(\mathrm{i}k_s)/4\pi$ is the acoustic pressure 1 m from the source and $k_s$ is the wavenumber of the medium at the location of the source.

	For convenience, in SPEC3D, the spectral truncation orders of all the examples are set to $N_w=N_b=20$; in an actual simulation, however, the user can arbitrarily specify the spectral truncation orders in the input file of SPEC3D. The number of discrete points in the $k_y$--domain is $N_q=512$. In the following seven numerical experiments, the speed of sound in seawater is set to 1500.0 m/s, the density of seawater is set to 1.0 g/cm$^3$, the speed of sound in the sediment is set to 1700.0 m/s, the sediment density is set to 1.5 g/cm$^3$, the attenuation coefficient in the sediment is set to 0.5 dB$/\lambda$, the sound speed in the acoustic half-space is set to 2000.0 m/s, the density is set to 2.0 g/cm$^3$, and the attenuation coefficient is set to 1.0 dB$/\lambda$.

	\subsection{Analytical example: an ideal three-dimensional fluid waveguide}
	Here, we use an ideal fluid waveguide, which is a highly simplified case in underwater acoustics, to validate SPEC3D. The structure of the three-dimensional ideal fluid waveguide is shown in Fig.~\ref{Ideal}. The frequency of the sound source is $f=20$ Hz, and the source is located at $z_s=36$ m. The number of coupled modes is $M=2$. Due to the horizontal independence of the ocean environment, the sound field must be cylindrically symmetrical about the axis where the sound source is located. First, consider the case when the lower boundary is also the pressure-release boundary. Fig.~\ref{Figure5} presents the sound field calculated by SPEC3D, featuring obvious cylindrical symmetry.
	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale = 0.8,domain=-5:5]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(4,-2.5)--(4,-5.5)--(2,-3)--cycle;
			\fill[cyan,opacity=0.7] (4,-5.5)--(14,-5.5)--(14,-2.5)--(4,-2.5)-- cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-5.5) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(5,-3.75) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-3.02);
			\draw[very thick](4,-2.5)--(4,-5.5);
			\draw[very thick](2,-3)--(4,-5.5);
			\draw[dashed, very thick](12,0)--(12,-3.0);
			\draw[very thick](14,-2.5)--(14,-5.5);
			\draw[dashed, very thick](12,-3.0)--(14,-5.5);
			\draw[very thick] (4,-5.5)--(14,-5.5);
			\draw[dashed, very thick](2,-3)--(12,-3);
			\filldraw [red] (3,-2.5) circle [radius=2.5pt];
			\node at (2.7,-1.3){$o$};
		\end{tikzpicture}
		\caption{Schematic diagram of the three-dimensional ideal fluid waveguide.}
		\label{Ideal}
	\end{figure}

	\begin{figure}[htbp]
		\centering\includegraphics[width=10cm]{Figure5}
		\caption{Sound field of the three-dimensional ideal fluid waveguide calculated by SPEC3D.}
		\label{Figure5}
	\end{figure}

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure6a}\label{Figure6a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure6b}\label{Figure6b}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure6c}\label{Figure6c}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure6d}\label{Figure6d}}
		\subfigure[]{\includegraphics[width=13cm]{Figure6e}\label{Figure6e}}
		\caption{Sound fields of the three-dimensional ideal fluid waveguide with free bottom calculated by SPEC3D at depths of 36 m (a) and 80 m (b); analytical solution (c) and sound field calculated by SPEC3D on the $y$--plane at $y=0$ m (d); TL along the $x$-direction at the depth of the sound source.}
		\label{Figure6}
	\end{figure}

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure7a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure7b}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure7c}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure7d}}
		\subfigure[]{\includegraphics[width=13cm]{Figure7e}}
		\caption{Sound fields of the three-dimensional ideal fluid waveguide with rigid bottom calculated by SPEC3D at depths of 25 m (a) and 80 m (b); analytical solution (c) and sound field calculated by SPEC3D on the $y$--plane at $y=0$ m (d); TL along the $x$-direction at the depth of the sound source.}
		\label{Figure7}
	\end{figure}

	Figs.~\ref{Figure6a} and \ref{Figure6b} depict slices of the sound field calculated by SPEC3D at two depths. According to the above analysis, the rings formed by the sound field on each slice should be strictly concentric, but the sound fields are not standard rings, so there are some errors. There is an analytical solution for the ideal fluid waveguide with a pressure-release boundary. The analytical solution for each plane of symmetry is:
	\begin{subequations}
		\begin{gather}
			\label{eq.51a}
			p(x,z)=\frac{2\pi \mathrm{i}}{H}\sum_{n=1}^{\infty}\sin{(k_{z,n}z_s)}\sin{(k_{z,n}z)}\mathcal{H}_0^{(1)}(k_{x,n}x)\\
			k_{z,n}=\frac{n\pi}{H},\quad k_{x,n}=\sqrt{k_0^2-k_{z,n}^2}, \quad n=1,2,3,\cdots
		\end{gather}
	\end{subequations}
	Figs.~\ref{Figure6c} and \ref{Figure6d} show the analytical sound field and the sound field calculated by SPEC3D, respectively, indicating that the two sound fields are reasonably consistent. Fig.~\ref{Figure6e} plots the TL curve along the $x$-direction at the depth of the sound source. The consistency between the SPEC3D solution and the analytical solution verifies the correctness and reliability of SPEC3D.

	When the lower boundary is perfectly rigid, the analytical solution of the ideal fluid waveguide is the same as Eq.~\eqref{eq.51a}, except that the vertical wavenumber becomes:
	\begin{equation}
		k_{z,n}=\left(n-\frac{1}{2}\right)\frac{\pi}{H},\quad k_{x,n}=\sqrt{k_0^2-k_{z,n}^2},\quad n=1,2,3\dots
	\end{equation}
	The number of coupled modes is $M=3$. Similarly, the results of Fig.~\ref{Figure7} can lead to the same conclusion as Fig.~\ref{Figure6}.

	\subsection{Analytical example: an ideal wedge-shaped waveguide}

	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale = 0.8]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(4,-8.5)--(2,-6)--cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-8) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(4.5,-3.125) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-6.02);
			\draw[very thick](4,-2.5)--(4,-8.5);

			\draw[dashed, very thick](2,-6)--(12,0);
			\draw[very thick](4,-8.5)--(14,-2.5);

			\draw[very thick](2,-6)--(4,-8.5);
			\draw[dashed, very thick](4,-8.5)--(8,-8.5);
			\filldraw [red] (8,-2.7) circle [radius=2.5pt];
			\node at (2.7,-1.3){$0$};
			\draw[very thick](5,-8.5) arc (0:30:1);
			\node at (6.2,-8.2){$\theta_0$=2.86°};
		\end{tikzpicture}
		\caption{Schematic of the ideal wedge-shaped waveguide.}
		\label{Wedge}
	\end{figure}

	Consider the ideal wedge-shaped waveguide shown in Fig.~\ref{Wedge}, which is a primary benchmark problem for range-dependent waveguides. Both the sea surface and the seabed are pressure-release boundaries, and the source is located at $(x_s, 0, z_s)$. Buckingham \cite{Buckingham1984,Buckingham1987,Doolittle1988} developed an analytical solution to this benchmark problem, which was originally proposed and discussed at two consecutive conferences of the Acoustic Society of America (ASA). Here, we briefly introduce this problem. The velocity potential in the water body is:
	\begin{subequations}
		\label{eq.53}
		\begin{gather}
			\Phi=\frac{1}{\theta_{0}} \sum_{n=1}^{\infty} I_{\nu_n}\left(r, r_s, z\right) \sin (\nu_n \theta) \sin \left(\nu_n \theta_s\right) \\
			\label{eq.53b}
			I_{\nu_n}\left(r, r_s, z\right)=\mathrm{i} \int_{0}^{\infty} k_r \frac{\exp (\mathrm{i} \eta|z|)}{\eta} J_{\nu_n}(k_r r) J_{\nu_n}\left(k_r r_s\right) \mathrm{d} k_r\\
			\eta=\sqrt{k^2-k_r^2},\quad \nu_n = \frac{n\pi}{\theta_0}, \quad n=1,2,\cdots
		\end{gather}
	\end{subequations}
	where $r$ and $r_s$ are the distances from the receiver and source to the apex of the wedge, respectively; $\theta$ and $\theta_s$ are the angles measured from the apex to the depths of the receiver and source, respectively; $\theta_0$ is the wedge angle; $k$ is the wavenumber of the seawater; and $J_{\nu_n}(\cdot)$ is the Bessel function of order $\nu_n$. The TL field is calculated by the following formula:
	\begin{subequations}
		\label{eq.54}
		\begin{gather}
			\mathrm{TL}=-20 \log _{10}\bigg\vert\frac{\Phi\left(r, r_s, \theta, \theta_s\right)}{\Phi_{0}(1)}\bigg\vert \\
			\Phi_{0}(r)=\frac{\exp(\mathrm{i}k r)}{4\pi r}
		\end{gather}
	\end{subequations}
	Note that directly calculating the integral in Eq.~\eqref{eq.53b} may trigger numerical overflow, which can be approximated by taking the Debye asymptotic expansion \cite{Milton1972}. In this paper, we use the Gauss--Kronrod integration.

	\begin{figure}[htbp]
		\centering\includegraphics[width=13cm]{Figure9}
		\caption{Sound field of the three-dimensional wedge-shaped waveguide with a free bottom calculated by SPEC3D.}
		\label{Figure9}
	\end{figure}

	In this case, the frequency of the sound source is $f=25$ Hz, and the source is located at $x_s=2000$ m and $z_s=50$ m. The wedge is $x=4000$ m long in the horizont direction and $z=200$ m deep in the vertical direction. The number of coupled modes is $M=6$. Fig.~\ref{Figure9} shows the sound field of the three-dimensional wedge-shaped waveguide calculated by SPEC3D; the sound field exhibits obvious three-dimensional characteristics.

	Next, we compare slices through the sound field with slices from the analytical solution in detail. As observed from the slices through the sound field in Fig.~\ref{Figure10}, the sound field calculated by SPEC3D effectively matches the analytical sound field. Thus, despite the small (typically less than 1 dB) error at long distances, SPEC3D can can successfully reproduce the sound field of the waveguide.

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10b}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10c}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10d}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10e}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure10f}}
		\caption{Sound fields of the three-dimensional wedge-shaped waveguide with a free bottom calculated by the analytical solution (a) and SPEC3D (b) at a depth of $z=25$ m; sound fields calculated by the analytical solution (c) and SPEC3D (d) on the $y$--plane at $y=1000$ m; TLs along the $y$-direction (e) and $x$-direction at a depth of $z=25$ m.}
		\label{Figure10}
	\end{figure}

	The analytical solution for the three-dimensional wedge-shaped waveguide with a rigid bottom has the same form as Eq.~\eqref{eq.53}, but:
	\begin{equation}
		\nu_n = \left(n-\frac{1}{2}\right)\frac{\pi}{\theta_0}, \quad n=1,2,\cdots
	\end{equation}
	Fig.~\ref{Figure11} displays a comparison between the results of SPEC3D and the analytical solution, the number of coupled modes is $M=7$, revealing obvious agreement, which further confirms the correctness and reliability of SPEC3D.
	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11b}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11c}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11d}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11e}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure11f}}
		\caption{Sound fields of the three-dimensional wedge-shaped waveguide with a rigid bottom calculated by the analytical solution (a) and SPEC3D (b) at a depth of $z=10$ m; sound fields calculated by the analytical solution (c) and SPEC3D (d) on the $y$--plane at $y=1000$ m; TLs along the $y$-direction (e) and $x$-direction at a depth of $z=10$ m.}
		\label{Figure11}
	\end{figure}

	\subsection{Acoustic half-space}
	Here, we consider a example involving an acoustic half-space, as shown in Fig.~\ref{Halfspace}, where the seafloor topography is also $x$-independent. In this case, the frequency of the sound source is $f=50$ Hz, and the source is located at $x_s=0$ m and $z_s=36$ m, $h=50$ m and $H=100$ m. The number of coupled modes is $M=6$. Fig.~\ref{Figure13} roughly shows the three-dimensional sound field structure of the waveguide.
	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale=0.8,domain=-5:5]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(4,-2.5)--(4,-5.5)--(2,-3)--cycle;
			\fill[cyan,opacity=0.7] (4,-5.5)--(14,-5.5)--(14,-2.5)--(4,-2.5)-- cycle;
			\fill[orange,opacity=0.7] (2,-3)--(4,-5.5)--(4,-8.5)--(2,-6)--cycle;
			\fill[orange,opacity=0.7]  (4,-5.5)--(14,-5.5)--(14,-8.5)--(4,-8.5)--cycle;
			\fill[brown,opacity=0.7] (2,-6)--(4,-8.5)--(14,-8.5)--(14,-9.5)--(4,-9.5)--(2,-7)--cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-9.5) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(5,-3.75) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-7);
			\draw[very thick](4,-2.5)--(4,-9.5);
			\draw[very thick](2,-3)--(4,-5.5);

			\draw[dashed, very thick](12,0)--(12,-6.0);
			\draw[very thick](14,-2.5)--(14,-9.5);
			\draw[dashed, very thick](12,-3.0)--(14,-5.5);
			\draw[very thick] (4,-5.5)--(14,-5.5);
			\draw[dashed, very thick](2,-3)--(12,-3);
			\draw[dashed, very thick](12,-6.0)--(14,-8.5);
			\draw[dashed, very thick](2,-6.0)--(12,-6.0);

			\draw[very thick](2,-6)--(4,-8.5);
			\draw[very thick](4,-8.5)--(14,-8.5);
			\filldraw [red] (3,-3.5) circle [radius=2.5pt];
			\node at (2.7,-1.3){$o$};
		\end{tikzpicture}
		\caption{Schematic diagram of the three-dimensional horizontally independent waveguide with an acoustic half-space.}
		\label{Halfspace}
	\end{figure}

	\begin{figure}[htbp]
		\centering\includegraphics[width=10cm]{Figure13}
		\caption{Sound field of the three-dimensional horizontally independent waveguide with an acoustic half-space calculated by SPEC3D.}
		\label{Figure13}
	\end{figure}

	Figs.~\ref{Figure14a} and \ref{Figure14b} present cross-sections of the sound field calculated by SPEC3D at depths of 36 m and 50 m, respectively. The sound field still clearly exhibits the symmetry of horizontally independent terrain, but a certain degree of error remains at different azimuth angles. However, considering the additional energy loss caused by the leakage of the seabed in the half-space and the error due to the numerical calculation, this error is inevitable and acceptable. Figs.~\ref{Figure14c} and \ref{Figure14d} display the two-dimensional sound fields calculated by KRAKEN \cite{Kraken2001} and NM-CT \cite{Tuhw2020b,NM-CT}, respectively, and Fig.~\ref{Figure14e} plots a slice of the sound field computed by SPEC3D at $y=0$ m; the simulation results of all three programs exhibit good agreement. Additionally, Fig.~\ref{Figure14f} shows various slices through the three-dimensional sound field of the waveguide, and the cross-sections further illustrate the significant three-dimensional characteristics of the sound field.

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14a}\label{Figure14a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14b}\label{Figure14b}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14c}\label{Figure14c}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14d}\label{Figure14d}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14e}\label{Figure14e}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure14f}\label{Figure14f}}
		\caption{Sound fields of the three-dimensional horizontally independent waveguide with an acoustic half-space calculated by SPEC3D at depths of $z=36$ m (a) and $z=50$ m (b); sound fields calculated by KRAKEN (c), NM-CT (d) and SPEC3D (e) on the $y$--plane at $y=0$ m; three-dimensional sound field slices calculated by SPEC3D in the planes of $x=1000$ m, $y=0$ m, and $z=36$ m (f).}
	\end{figure}

	\subsection{Ridge-shaped waveguide}
	Consider a ridge-shaped waveguide, as shown in Fig.~\ref{Ridge}. In the $xoz$--plane, the expression for the seabed topography is given by the following equation:
	\begin{equation}
		h(x)= \begin{cases}65+25 \cos \left[\frac{\pi(x-500)}{200}\right], & 300<x<700, \\ 40, & \text { elsewhere, }\end{cases}
	\end{equation}
	The frequency of the sound source is $f=25$ Hz, and the source is located at $x_s=0$ m and $z_s=25$ m. The ocean depth $H$ is 100 m, and the lower boundary is the rigid seabed. Fig.~\ref{Figure16} is the three-dimensional sound field structure of the waveguide calculated by SPEC3D. Due to the flat terrain in both the near field and the far field, the sound field still exhibits a certain symmetry. The uplift of the ridge at 300--700 m induces obvious changes in the sound field within this region.

	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale=0.8,samples=500,domain=7:10.5]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(4,-2.5)--(4,-5.5)--(2,-3)--cycle;
			\fill[cyan,opacity=0.7]  (4,-2.5)--(4,-5.5)--(7,-5.5)--plot(\x,{-4.4-1.1*cos(103*(\x))})--(10.5,-5.5)--(14,-5.5)--(14,-2.5)-- cycle;
			\fill[orange,opacity=0.7] (2,-3)--(4,-5.5)--(4,-8.5)--(2,-6)--cycle;
			\fill[orange,opacity=0.7]  (4,-5.5)--(7,-5.5)--plot(\x,{-4.4-1.1*cos(103*(\x))})--(10.5,-5.5)--(14,-5.5)--(14,-8.5)--(4,-8.5)--cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-8) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(5,-3.75) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-6);
			\draw[very thick](4,-2.5)--(4,-8.5);
			\draw[dashed, very thick](12,0)--(12,-6);
			\draw[very thick](14,-2.5)--(14,-8.5);
			\draw[dashed, very thick](12,-3)--(14,-5.5);
			\draw[dashed, very thick](2,-6)--(12,-6);
			\draw[dashed, very thick](12,-6)--(14,-8.5);

			\draw[very thick](2,-6)--(4,-8.5);
			\draw[very thick](4,-8.5)--(14,-8.5);
			\filldraw [red] (3,-2.7) circle [radius=2.5pt];
			\node at (2.7,-1.3){$o$};

			\draw[very thick](2,-3)--(4,-5.5);
			\draw[very thick](4,-5.5)--(7,-5.5);
			\draw[very thick, smooth] plot(\x,{-4.4-1.1*cos(103*(\x))});
			\draw[very thick](10.48,-5.5)--(14,-5.5);
			\draw[dashed, very thick](2,-3)--(5,-3);
			\draw[dashed, very thick, smooth] plot(\x-2,{-1.9-1.1*cos(103*(\x))});
			\draw[dashed, very thick](8.48,-3)--(12,-3);

		\end{tikzpicture}
		\caption{Schematic diagram of the three-dimensional ridge-shaped waveguide.}
		\label{Ridge}
	\end{figure}

	\begin{figure}[htbp]
		\centering\includegraphics[width=10cm]{Figure16}
		\caption{Sound field of the three-dimensional ridge-shaped waveguide calculated by SPEC3D.}
		\label{Figure16}
	\end{figure}

	For the cases where there is no analytical solution, we use the results of the benchmark program `COACH' developed by Liu et al as a reference \cite{Liuw2021}. COACH is a three-dimensional finite difference program with fourth-order accuracy developed for the ocean acoustic Helmholtz equation, which can be used to address arbitrary bathymetry and provide more accurate benchmark solutions for other three-dimensional underwater acoustic approximation models. The derivatives in the acoustic Helmholtz equation are numerically discretized based on regular grids, and the perfectly matched layer is introduced to absorb unphysical reflections from the boundaries where Sommerfeld radiation conditions are deployed. Here, running COACH is equivalent to directly numerically solving Eq.~\eqref{eq.1}.

	Figs.~\ref{Figure17a} and \ref{Figure17b} show the sound field slices calculated by COACH and SPEC3D at $z=25$ m, respectively. The results of SPEC3D and COACH are in good agreement, and even in a finer comparison of the TL curves in Fig.~\ref{Figure17c}, we can see that the errors of both are satisfactory. In addition, the backscattering and forward scattering generated by the uplifted topography of the ridge show how the ridge perturbs the near-field and far-field sound fields.

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure17a}\label{Figure17a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure17b}\label{Figure17b}}
		\subfigure[]{\includegraphics[width=13cm]{Figure17c}\label{Figure17c}}
		\caption{Sound fields of the three-dimensional ridge-shaped waveguide calculated by COACH (a) and SPEC3D (b) at $z=25$ m; TLs along the $x$-direction at a depth of $z=25$ m (c).}
	\end{figure}

	\subsection{Trench-shaped waveguide}
	In contrast to the previous example, Fig.~\ref{Trench} shows a marine environment with a trench-shaped waveguide. The seafloor topography is given by:
	\begin{equation}
		h(x)= \begin{cases}50-25 \cos \left[ \frac{\pi(x-500)}{200}\right], & 300<x<700, \\ 75, & \text { elsewhere, }\end{cases}
	\end{equation}
	The remainder of the configuration is exactly the same as the above example. Fig.~\ref{Figure19} is the three-dimensional sound field structure of the waveguide calculated by SPEC3D. The horizontal independence of both the near field and the far field makes the sound field retain a certain symmetry. However, the occurrence of refraction in the trench causes the sound field to exhibit an obvious three-dimensional effect, which further affects the near and far fields through forward scattering and backscattering.

	Likewise, Figs.~\ref{Figure20a} and \ref{Figure20b} plot sound field slices calculated by COACH and SPEC3D at $z=25$ m, respectively. Fig.~\ref{Figure20c} displays the TL curves calculated by the two programs along the $x$ direction at a depth of $z=25$ m. In both slice and line plots, the results of SPEC3D are very similar to those of COACH, which directly solves the Helmholtz equation. Therefore, we believe that the algorithm proposed in this paper is accurate and efficient.

	\begin{figure}[htbp]
		\centering
		\begin{tikzpicture}[node distance=2cm,scale=0.8,samples=500,domain=7:10.5]
			\fill[cyan,opacity=0.7] (2,0)--(12,0)--(14,-2.5)--(4,-2.5)--(4,-5.5)--(2,-3)--cycle;
			\fill[cyan,opacity=0.7]  (4,-2.5)--(4,-5.5)--(7,-5.5)--plot(\x,{-6.6+1.1*cos(103*(\x))})--(10.5,-5.5)--(14,-5.5)--(14,-2.5)-- cycle;
			\fill[orange,opacity=0.7] (2,-3)--(4,-5.5)--(4,-8.5)--(2,-6)--cycle;
			\fill[orange,opacity=0.7]  (4,-5.5)--(7,-5.5)--plot(\x,{-6.6+1.1*cos(103*(\x))})--(10.5,-5.5)--(14,-5.5)--(14,-8.5)--(4,-8.5)--cycle;
			\draw[thick, ->](3,-1.25)--(14,-1.25) node[right]{$x$};
			\draw[thick, ->](3,-1.25)--(3,-8) node[below]{$z$};
			\draw[very thick](1.98,0)--(12.02,0);
			\draw[very thick](4,-2.5)--(14.02,-2.5);
			\draw[very thick](2,0)--(4,-2.5);
			\draw[thick, ->](2,0)--(5,-3.75) node[below]{$y$};
			\draw[very thick](12,0)--(14,-2.5);
			\draw[very thick](2,0.02)--(2,-6);
			\draw[very thick](4,-2.5)--(4,-8.5);
			\draw[dashed, very thick](12,0)--(12,-6);
			\draw[very thick](14,-2.5)--(14,-8.5);
			\draw[dashed, very thick](12,-3)--(14,-5.5);
			\draw[dashed, very thick](2,-6)--(12,-6);
			\draw[dashed, very thick](12,-6)--(14,-8.5);

			\draw[very thick](2,-6)--(4,-8.5);
			\draw[very thick](4,-8.5)--(14,-8.5);
			\filldraw [red] (3,-2.7) circle [radius=2.5pt];
			\node at (2.7,-1.3){$o$};

			\draw[very thick](2,-3)--(4,-5.5);
			\draw[very thick](4,-5.5)--(7,-5.5);
			\draw[very thick, smooth] plot(\x,{-6.6+1.1*cos(103*(\x))});
			\draw[very thick](10.48,-5.5)--(14,-5.5);
			\draw[dashed, very thick](2,-3)--(5,-3);
			\draw[dashed, very thick, smooth] plot(\x-2,{-4.1+1.1*cos(103*(\x))});
			\draw[dashed, very thick](8.48,-3)--(12,-3);

		\end{tikzpicture}
		\caption{Schematic diagram of the three-dimensional trench-shaped waveguide.}
		\label{Trench}
	\end{figure}

	\begin{figure}[htbp]
		\centering\includegraphics[width=10cm]{Figure19}
		\caption{Sound field of the three-dimensional trench-shaped waveguide calculated by SPEC3D.}
		\label{Figure19}
	\end{figure}

	\begin{figure}[htbp]
		\centering
		\subfigure[]{\includegraphics[width=6.5cm]{Figure20a}\label{Figure20a}}
		\subfigure[]{\includegraphics[width=6.5cm]{Figure20b}\label{Figure20b}}
		\subfigure[]{\includegraphics[width=13cm]{Figure20c}\label{Figure20c}}
		\caption{Sound fields of the three-dimensional trench-shaped waveguide calculated by COACH (a) and SPEC3D (b) at $z=25$ m; TLs along the $x$-direction at a depth of $z=25$ m (c).}
	\end{figure}

	\section{Discussion, Remarks and Conclusion}
	\subsection{Discussion}
	According to the above derivation and analysis, we can obtain a quasi-three-dimensional sound field through an inverse Fourier transform as long as we solve the sound field in Fig.~\ref{Figure2} corresponding to different $k_y$. It is worth noting that $k_y$ is a set of complex constants and that the marine environment parameters in Fig.~\ref{Figure2} do not change for different $k_y$, which means that the stepwise segmentation of the $x$-dependences in Fig.~\ref{Figure2} also does not change. Let us re-examine Eq.~\eqref{eq.7}, which becomes the matrix eigenvalue problem in Eq.~\eqref{eq.39} after spectral discretization. According to linear algebra theory, Eq.~\eqref{eq.39} is equivalent to the following matrix eigenvalue problem:
	\begin{equation}
		\label{eq.58}
		\left(\frac{4}{\vert\Delta h\vert^2}\mathbf{C}_{\rho}\mathbf{D}_{N}\mathbf{C}_{1/\rho}\mathbf{D}_{N}+\mathbf{C}_{k^2}\right)\bm{\hat{\Psi}}
		= (k_y^2+\mu^2) \bm{\hat{\Psi}}
	\end{equation}
	Intriguingly, the above indicates that the sum of squares of $k_y$ and $\mu$ is always equal to the eigenvalue in the above equation and that the eigenvectors do not change with $k_y$. Therefore, it seems that Eq.~\eqref{eq.58} needs to be solved only once after the derivation of Sec.~\ref{section3.2} to obtain the eigenvalues of $\mu$ values corresponding to different $k_y$. Since the eigenvectors are invariant, the coupling submatrix in Eq.~\eqref{eq.22} is also shared for different $k_y$.

	However, in fact, the above conjecture includes an obvious misunderstanding, that is, for different $k_y$, it is not Eq.~\eqref{eq.39} that needs to be solved, but Eqs.~\eqref{eq.43} or \eqref{eq.46} after boundary conditions are applied. Obviously, Eqs.~\eqref{eq.43} and \eqref{eq.46} cannot be changed, similar to Eq.~\eqref{eq.48}, so in the $k_y$--domain, $N_q$ sound fields still need to be calculated strictly according to the derivation in Secs.~\ref{section2.2} and \ref{section3.2}.

	\subsection{Remarks}
	From the above analysis, we can directly summarize the following advantages of the algorithm and program developed in this article:
	\begin{enumerate}

		\item
		The two-dimensional problem formed after implementing the Fourier transform can be solved parallel because the calculation of $\tilde{p}(x,k_y(q),z)$ is independent. This computational efficiency makes this algorithm more useful than similar algorithms.

		\item
		Furthermore, the two-dimensional problem formed after implementing the Fourier transform can be solved efficiently by using the global matrix method. This model is considered to involve two-way coupled normal modes, and its accuracy is higher than that of a model with only one-way coupled modes.

		\item
		The Chebyshev--Tau spectral method used for the single-segment solution can accurately simulate free, rigid and half-space seafloor boundaries, especially for a half-space, and this approach avoids the leakage problem of traditional iterative root-finding algorithms.

	\end{enumerate}

	\subsection{Conclusion}
	In this paper, we develop a novel algorithm and model that can solve for quasi-three-dimensional sound fields. The results of numerical simulations verify the correctness and reliability of the model and code. In the algorithm, the Fourier transform is used to transform the equation governing quasi-three-dimensional acoustic propagation from the $y$--domain to the $k_y$--domain. After the sound pressure at each discrete $k_y$ is obtained, the sound pressure in the $y$--domain can be obtained by applying an inverse Fourier transformation; this process can also be regarded as $k_y$ wavenumber integration to some extent. The governing equation in the $k_y$--domain adopts a similar form to that of a two-dimensional line source sound field. The stair-step approximation strategy is used to incorporate the $x$-dependence of the two-dimensional sound field, and the coupling coefficients are assembled by means of a global matrix, a good sparse band matrix that can be solved efficiently. In each $x$-independent segment, the Chebyshev--Tau spectral method is used to solve for the eigenpairs, and the proposed method achieves good accuracy, speed and robustness.

	In terms of its application scope, this model requires that the marine environment be independent in the $y$ direction, which limits the practicality of SPEC3D to a certain extent. Therefore, there is considerable promise in the development of a highly precise algorithm for determining the real three-dimensional sound field of an ocean environment with random and arbitrary changes based on spectral methods.

	\section*{Acknowledgments}
	The authors would like to thank Prof. Wenyu Luo from the Institute of Acoustics, Chinese Academy of Sciences, for his valuable guidance on the analytical solution of the sound field for the wedge-shaped waveguide.

	This work was supported by the National Key Research and Development Program of China [grant number 2016YFC1401800] and the National Natural Science Foundation of China [grant numbers 61972406].

\input{main.bbl}
\end{document}

redpajama_set_name: RedPajamaArXiv


text:
AlarmForce Announces Record Fourth Quarter Results for 2010
January 27, 2011 16:30 ET | Source: AlarmForce Industries Inc.
TORONTO, ONTARIO--(Marketwire - Jan. 27, 2011) - AlarmForce Industries Inc. (TSX:AF), Canada's largest manufacturer and installer of live two-way voice home alarms systems, is pleased to announce annual results for the year ended October 31, 2010:
October 31, 2010 October 31, 2009 Change
(year ended) (year ended) (year over year)
Revenue $ 37,187,000 $ 34,133,000 +9 %
Gross profit $ 28,954,000 $ 26,101,000 +11 %
Income before income taxes $ 7,052,000 $ 5,673 ,000 +24 %
Net income $ 4,782,000 $ 3,780 ,000 +27 %
Cash flow from operations $ 7,209,000 $ 6,467,000 +11 %
Cash flow from operations excluding working capital $ 8,294 ,000 $ 6,718,000 +23 %
EBITDA $ 10,813,000 $ 9,057,000 +19 %
EBITDA/share $ 0.88 $ 0.74 +19 %
Basic Net income / share $ 0.39 $ 0.31 +25 %
Diluted Net income / share $ 0.39 $ 0.31 +25 %
* EBITDA (Earnings Before Interest, Income Taxes, Depreciation, and Amortization) is a key measure in the security industry and should not be interpreted as GAAP)
EBITDA is defined as earnings before interest expenses, income taxes, depreciation and amortization. EBITDA is a standard measure used in the security industry to assist in understanding and comparing operating results and is often referred to by our competitors. Management views EBITDA as an important measure of operating performance of the Company. Yet, since it does not have any standardized meaning defined by Canadian GAAP, it may not be considered in isolation of GAAP measures such as net income/loss or cash flows, as a measure of liquidity. The Company, however, utilizes these measures in making operating decisions and assessing its performance. Management believes this to be an important measure, as it allows the Company to assess its ongoing business without the impact of depreciation or amortization expenses. Since EBITDA is not a defined term under Canadian GAAP, it is unlikely to be comparable to similar measures presented by other issuers.
Joel Matlin, President and CEO, announced that AlarmForce Industries Inc. added over 23,800 new accounts in the fiscal year ended October 31, 2010, closing with a subscriber base of 113,500 accounts.
He stated that "The net growth rate in the total subscriber base was 11%, which generated additional monthly recurring revenues for the Company and resulted in a record $37.2 million of total revenues. The increase in revenues from the previous year was $3.0 million or 9%. Gross profit increased by 11% for the year, and as a percentage of sales it increased from 76% to 78%. The growth in subscribers, accompanied by the increase in operating efficiency reflected in the gross profit margin, helped us to achieve a record for income and earnings per share in 2010."
Earnings per share increased from 31 cents per share to 39 cents per share, the highest recorded in our 22 years of business. Net income was up by 27%, increasing by approximately $1 million to $4.8 million in 2010. We are continuing to grow new subscribers profitably, fueling growth in earnings at a higher rate than subscriber growth and revenues.
Corresponding to the growth in subscribers and recurring revenues, the income before taxes increased from $5.7 million to $7.1 million, a record figure, up by 24%. EBITDA increased by 19% from $9.1 million to $10.8 million. We achieved these increases in spite of the impact of higher research and development expenses in 2010. The reported figures are net of the cost of research and development expenses on the video product which amounted to $0.8 million in 2010.
Due to the fact that we create organic accounts, rather than acquiring existing accounts, the costs of all our advertising programs are expensed, reducing operating income and EBITDA in that given year. Therefore, as an alternative indicator, the results may also be compared by excluding marketing expenses. Excluding these expenses, EBITDA increased from $18.1 million to $20.1 million in 2010, an 11% increase.
Our cash flow from operations increased by 23% to $8.3 million before the impact of increases in working capital, and increased by 11% to $7.2 million after the impact of working capital requirements, which reflects the timing of expenditures related to growth. Cash used for investing activities totaled $8.9 million for the year and cash used for debt repayment totaled $1.1 million.
We have continued to finance all expenditures for growth in subscriber installations and brand expansion in the US out of internal cash flows, resulting in no outstanding debt and a strong balance sheet, which will allow us to pursue our goals in the next few years. Our cash requirements include expenditures for the development and marketing of our application- based remote video system,. In addition, we are continually enhancing our sales and customer support staff to facilitate the growth.
We estimate the development of the application-based remote video system to be completed in the first half of 2011 and we intend to advertise and launch our enhanced services in the second half of the year. The enhanced services will be offered to both existing AlarmForce subscribers and to new subscribers in the general consumer market, in all our service territories. This expansion of our market reach is a significant, natural extension of our brand.
AlarmForce provides security and personal emergency response monitoring and related services to residential and commercial subscribers throughout Canada and selected centers across the United States. The Company is a leading provider of two-way voice alarm systems in Canada and sole provider of AlarmPlus, a wireless line-cut protection system. More information about the Company's products and services can be found at www.alarmforce.com.
Contact Information: AlarmForce Industries Inc.
Investor Relations Department
investorrelations@alarmforce.com
www.alarmforce.com
AlarmForce Industries Inc.

redpajama_set_name: RedPajamaCommonCrawl


Are you contemplating the kind of drivetrain in your next automobile? Toyota – The second largest automaker in the world is Toyota, maker of the Camry, the Corolla, and a number of SUVs, vans, passenger automobiles, and a van. It was once mentioned that a retailer could earn extra revenue from a customer over the three years after buying their car than was realised from the original sale.
Lubrication Techniques: Automotive lubrication programs include several different types of oil meters, ATF meters, gear lube meters, dispense valves and grease dispenses valves. Income from the sale of spare elements; oils and lubricants; paint and materials; and sublet and varied are all subsidiary to the buying and promoting of productives' time.
Though car maintenance is intently associated to your vehicle make and mannequin, there are some upkeep points which can be pretty much frequent whether you might have a Cadillac, VW or a Mercedes, for instance changing oil, changing air filters, engine wash, fluid maintenance, tier maintenance and so forth.
New autos repairs are sometimes covered by warranties. However their buyer base is thousands and thousands of motorists fairly than a dozen insurance firms, so service managers can set whatever times they need – within reason, and naturally, subject to competitors.
Heavy automobiles might be raised by way of this equipment. Due to this fact the automobile corporations begin slicing the prices, especially fixed costs, like closing the unprofitable facilities, eliminating jobs. Toyota's Lexus division produces luxurious cars, while its Scion division manufactures youth oriented autos.


Hello Connie and thanks for asking the Disney Parks Moms Panel such a great question.
Perhaps you could try making your FastPass+ selections online with your My Disney Experience account; Walt Disney World Resort Guests are able to make FastPass+ selections up to 60 days in advance of their trip via the easy to use My Disney Experience online service.
As you rightly suggest if you wish, you can make your FastPass+ choices on the day of your visit at one of the many conveniently placed FastPass+ kiosks available at all the Disney Theme Parks. For your added convenience be sure to pick up a Park Guide Map where kiosks are clearly located and should you require any further assistance you will find plenty of knowledgeable Cast Members on hand to help.
If you think of any more Walt Disney World Resort vacation planning questions that I can assist you with please don't hesitate to ask, I'm always happy to help.


text:
Become a puzzling Mastermind with Peaky Blinders on Xbox One, PS4, Switch and PC
Neil Watton (neil363)
It's pretty safe to say that Peaky Blinders has taken the TV watching audience by storm. Bursting onto the BBC and then Netflix, the exploits of the Shelby crew have captivated many, with no less than five series of the show already under its belt and a couple more still promised. But if you can't wait for the show to return for more viewing pleasures, there's another route you can take – get involved with the Shelby's in Peaky Blinders: Mastermind on Xbox One, PS4, Nintendo Switch and PC.
Available to purchase and download today, Peaky Blinders: Mastermind comes to Xbox One, PS4, Nintendo Switch and PC via the teams at FuturLab and Curve Digital. Having been created in close collaboration with the production teams behind the show, Mastermind will see you partaking in a puzzle-adventure that follows the goings on of the Shelby's as they hear of a plot which looks set to put Tommy out of business.
Working multiple strategies is going to be par for the course here with Peaky Blinders: Mastermind as you take direct control of key members of the Shelby family simultaneously over the same timeline, moving characters into position at exactly the right moment in order to navigate challenging scenarios. With the opportunity to rewind each character's path in order to fine tune your tactics, you'll be left to slowly solve all manner of mysteries. And if you don't? Well, that could well be the end of the Shelbys.
TRUE TO THE TV SHOW Take control of Tommy, Arthur, Polly, John, Ada and Finn as they navigate through familiar locations including The Garrison Tavern, Charlie Strong's Yard, Shelby Parlour and more.
THE RISE OF TOMMY When a Peaky Blinder is implicated in the murder of a Chinese opium dealer, Tommy embarks on a mission to find out the truth, slowly growing into the true Peaky Blinders' Mastermind.
BECOME THE MASTERMIND Exercise Tommy's power to plan complex scenarios in his head, by taking control of each key member of the Shelby family and perform their actions in parallel.
MANIPULATE TIME Each character's actions are automatically logged in a timeline which you can freely reset, fast-forward and rewind, until the plan is perfectly adjusted.
STEALTH AND STRATEGY Each mission is an increasingly complex puzzle to be solved undetected and on time. Choose the right characters and abilities to coordinate and execute the optimum plan.
Priced at £19.99, you'll find Peaky Blinders: Mastermind available today on Xbox One, PS4, Nintendo Switch and PC. The usual digital stores will sort you out with a download. Our full review is right here if you wish to know more about it.
Peaky Blinders: Mastermind is a puzzle-adventure game, based on the multi-award-winning BBC and Netflix TV show. Welcome to Birmingham, during the aftermath of the Great War. Set right before the events of Season One, Peaky Blinders: Mastermind lets you join the Shelby family's criminal gang. Follow the rise of Tommy as he successfully uncovers a sinister plot to put the family out of business, proving himself worthy of being the true mastermind of the Peaky Blinders. Mastermind refers to Tommy's ability to plan complex scenarios in his head. As the player, you get to exercise this power by taking control of key members of the Shelby family including Tommy, Arthur, Polly and more. Become the Mastermind as you freely reset and rewind each character's path in order to tune all of their actions for perfect coordination. Achieve the highest mission rating by strategically utilising characters and stealth to solve increasingly complex puzzles with optimum timing.
Previous articleMultiple new DLC pieces arrive in MY HERO ONE'S JUSTICE 2
Next articleStealth SX-C100 Twin Charging Dock for Xbox One Review
https://www.thexboxhub.com
An Xbox gamer since 2002, I bought the big black box just to play Oddworld: Munch's Oddysee. I have since loved every second of the 360's life and am now just as obsessed with the Xbox One - mostly with the brilliant indie scene that has come to the fore. Gamertag is neil363, feel free to add me to your list.
eight + 2 =
Q.U.B.E. 2 Review
Animal Super Squad Review


Vacant for over 20 years, the Clyde Theatre, in Quimby Village, is on its way to becoming a premier arts and entertainment destination.
Portion of interior reflecting theatre's original grandeur.
With the recent flurry of downtown redevelopment, renovation of the long-vacant Clyde Theatre, located on Fort Wayne's south side is going full speed ahead.
Rick Kinney is honored to have new partners, Chuck and Lisa Surack on board to see his dreams for the Clyde Theatre fast becoming a reality.
Rick Kinney, founder of Even Keel Event Productions, and his new equity partners, Chuck Surack, CEO of Sweetwater Sound and his wife Lisa, are eager to make Rick's dream of establishing one of Indiana's premier arts and entertainment venues a reality.
Once converted into a performance hall, the space will accommodate a standing-room-only crowd of 2,220 or 612 in seated format. Longer term plans will address refurbishing of the surrounding Quimby Village retail spaces to attract local entrepreneurs. View the entire Quimby Village complex here.
Earlier this week, WBOI's Julia Meek met with Kinney to discuss the scope of the project, and the impact it will make on Fort Wayne's cultural landscape.
Listen to an extended version of Julia Meek's interview with Rick Kinney on the Clyde Theatre renovation project.
For more information on the renovation and project updates, visit the Clyde Theatre website.
The African Children's Choir, with its parent organization, Music for Life, aims to help their country's most vulnerable children today, so they can help Africa grow in the future.
This Saturday and Wednesday, these 18 talented young singers will be performing in Harlan and Berne, Indiana.
Since 1904, the 1000 block of Rudisill Boulevard, on the city's south side, has been home to a variety of colleges and universities, and thousands of local graduates.
In 2011, Ambassador Enterprises purchased the campus and created The Summit, with the purpose of creating a space for like-minded partner organizations to honor that past, while embracing the whole community, and its future.
Fort Wayne landscaper, Jim Saxton, has a passion for rescuing historic lighting fixtures, which he traces back to his early teen years. That's when his job at the Paramount Theatre sparked an interest in lighting, as well as the city's rich past.
Enter Fort Wayne's Monthly Storytelling Showcase By Way Of "The Trap Door"
"The Trap Door," a monthly storytelling showcase here in Fort Wayne, has been going strong since it recently opened in September. So strong, in fact that as of Thursday evening it will find a permanent home in Artlink, to accommodate the ever-growing audiences in attendance.


text:
Congregation Sharey Tzedek Synagogue (שערי צדק) is a historic Jewish synagogue at 833 South 200 East in Salt Lake City, Utah, United States.

History
A group of mostly of Yiddish-speaking immigrants from Russia and elsewhere in Eastern Europe split off from Congregation Montefiore in 1916 to form Congregation Sharey Tzedek, Utah's third Jewish congregation. The group raised $33,200 to buy the land and build a synagogue. The building was built by general contractor John E. Anderson, and it was completed in 1920. Utah Governor Simon Bamberger spoke at the dedication on March 28, 1920. When the synagogue was built, its front had circular windows that contained colored glass with a star of David. Rabbi Joseph Strinkomsky served as the congregation's first rabbi. Services were in line with Orthodox Judaism, with different seating areas for men and women. Most of the members lived within a short walk of the synagogue, and they shopped at a kosher butcher nearby.

The congregation disbanded in the 1930s, and the building was sold to the Veterans of Foreign Wars in 1948. At the time of the building's sale, many of the original members of the congregation had died, and others were attending religious services at B'nai Israel, which followed Reform Judaism, or at Montefiore Synagogue, which followed Conservative Judaism.

A new Veterans of Foreign Wars post, led by Guy Snyder and Post Commander Lloyd S. Grover, began using the building.

The building was added to the National Register of Historic Places in 1985.

According to the Utah Historical Society, the building is significant "for its historical association with Utah's pluralistic community."

See also

B'nai Israel Temple - Synagogue belonging to the first Jewish congregation in the Salt Lake City area.
Congregation Montefiore Synagogue - Synagogue belonging to the second Jewish congregation in the Salt Lake City area.

References

External links

Ashkenazi Jewish culture in the United States
European American culture in Utah
National Register of Historic Places in Salt Lake City
Synagogues completed in 1920
Romanesque Revival architecture in Utah
Russian-Jewish culture in the United States
Synagogues in Salt Lake City
Synagogues on the National Register of Historic Places in Utah

redpajama_set_name: RedPajamaWikipedia