# Scattering of elastic waves by a spheroidal inclusion

Scattering of elastic waves by a spheroidal inclusion Summary An analytical solution is presented for scattering of elastic waves by prolate and oblate spheroidal inclusions. The problem is solved in the frequency domain where separation of variables leads to a solution involving spheroidal wave functions of the angular and radial kind. Unlike the spherical problem, the boundary equations remain coupled with respect to one of the separation indices. Expanding the angular spheroidal wave functions in terms of associated Legendre functions and using their orthogonality properties leads to a set of linear equations that can be solved to simultaneously obtain solutions for all coupled modes of both scattered and interior fields. To illustrate some of the properties of the spheroidal solution, total scattering cross-sections for P, SV and SH plane waves incident at an oblique angle on a prolate spheroid, an oblate spheroid and a sphere are compared. The waveforms of the scattered field exterior to the inclusion are calculated for these same incident waves. The waveforms scattered by a spheroid are strongly dependent upon the angle of incidence, are different for incident SV and SH waves and are asymmetrical about the centre of the spheroid with the asymmetry different for prolate and oblate spheroids. Wave propagation, Wave scattering and diffraction 1 INTRODUCTION The situation of an inclusion having material properties different from the surrounding medium is commonly encountered in studies of elastic wave propagation. When the shape of the inclusion has all dimensions approximately equal, the known solutions for scattering of elastic waves by a sphere are useful (Korneev & Johnson 1993a, 1996). However, in many cases, such as markedly elongated or flattened inclusions, the approximation by a sphere is questionable and more appropriate solutions are needed. The next step up from a sphere in terms of complexity is a spheroidal inclusion, which has the advantage of being a good approximation to many elongated or flattened shapes while still leading to differential equations that are separable. It appears that a complete solution for the scattering of elastic waves by a spheroidal inclusion has not yet appeared in the literature, so that is the purpose of the present study. Unlike the case of elastic waves, scattering of electromagnetic waves by a spheroidal inclusion has a well-developed literature (see Mishchenko et al. 2000, for a review). The elastic wave problem is more complicated than the electromagnetic wave problem in that it involves two velocities instead of one, requires two potential functions instead of one, allows three types of incident waves instead of two, and the relationship between displacement and stress fields lacks the symmetry of the relationship between electric and magnetic fields. However, essentially all of the methods and mathematics that have been developed for the electromagnetic problem can be applied to the elastic problem. Of the various methods that have been developed for electromagnetic problems, the present study follows most closely that of Asano & Yamamoto (1975). The approach followed in this study is to write the differential equations for elastic displacement in spheroidal coordinates and construct separate solutions to these equations for the incident field, the scattered field outside the inclusion and the interior field within the inclusion. The arbitrary constants in these fields are determined by requiring that boundary conditions on the surface of the inclusion are satisfied. This leads to a set of matrix equations than can be solved to yield an exact analytical solution to the scattering problem. There are a number of alternative approaches to the elastic scattering problem for a spheroidal inclusion that yield approximate solutions. One such approach is the transition matrix method that has been used by Hackman (1984) and Hackman & Lim (1994) to obtain the scattering solutions for inclusions of various shapes including a prolate spheroid. Another possible approach is to formulate the scattering problem in terms of an equivalent integral equation. While the integral formulation does not generally lead to an exact analytical solution, it is useful as the starting point for various approximations such as the Born series, low frequency, low contrast in material properties and far field. Examples of such approximate solutions can be found in Miles (1960), Gubernatis et al. (1977a,b), Datta (1977), Wu & Aki (1985) and Dassios & Kleinman (2000). 2 BASIC EQUATIONS IN SPHEROIDAL COORDINATES Consider a spheroidal inclusion with a polar axis rp having rotational symmetry and a perpendicular equatorial axis re. The degree to which the shape of the spheroid departs from that of a sphere is characterized by the interfocal distance h defined as   $$h = 2 \sqrt{r_p^2 - r_e^2} .$$ (1)It is also convenient to use the aspect ratio defined as   $$R_a = \frac{\text{max} \lbrace r_p , r_e\rbrace }{ \text{min} \lbrace r_p , r_e \rbrace } .$$ (2)When rp > re the spheroid is prolate and when rp < re it is oblate. Prolate and oblate spheroids are sketched in Figs 1 and 2, respectively, along with the geometry that is used later in numerical calculations. Figure 1. View largeDownload slide 2-Dview of a 3-D prolate spheroid having a polar axis rp and an equatorial axis re. An incident plane wave propagates from below in a positive z-direction that makes an angle ζ with the rp axis. Observations of the scattered waves are registered in the x, y and z directions on a plane that is displaced a distance zo from the centre of the spheroid. Figure 1. View largeDownload slide 2-Dview of a 3-D prolate spheroid having a polar axis rp and an equatorial axis re. An incident plane wave propagates from below in a positive z-direction that makes an angle ζ with the rp axis. Observations of the scattered waves are registered in the x, y and z directions on a plane that is displaced a distance zo from the centre of the spheroid. Figure 2. View largeDownload slide A sketch similar to Fig. 1 for an oblate spheroid. Figure 2. View largeDownload slide A sketch similar to Fig. 1 for an oblate spheroid. The spheroidal inclusion is imbedded in a homogeneous background having material properties density ρ, bulk modulus κ and shear modulus μ. The elastic compressional velocity is $$\alpha = [ (\kappa +\frac{4}{3}\mu )/\rho ]^{1/2}$$ and the shear velocity is β = [μ/ρ]1/2. Material properties within the inclusion will in general be different from the background and are denoted with a prime, so they are given by ρ΄, κ΄, μ΄, α΄, β΄. Considering first the prolate case, the spheroidal coordinates (η, ξ, ϕ) can be defined as an orthogonal right-handed system with their relationship to rectangular coordinates (x1, x2, x3) given by (Morse & Feshbach 1953; Flammer 1957)   $$x_1 = \frac{h}{2} \left[(\xi ^2-1)(1-\eta ^2)\right]^{1/2} \cos \phi ,$$ (3a)  $$x_2 = \frac{h}{2} \left[(\xi ^2-1)(1-\eta ^2)\right]^{1/2} \sin \phi ,$$ (3b)  $$x_3 = \frac{h}{2} \eta \xi .$$ (3c) with x3 parallel to the axis of rotational symmetry. These spheroidal coordinates have the ranges −1 ≤ η ≤ 1, 1 ≤ ξ ≤ ∞ and 0 ≤ ϕ ≤ 2π. The particle displacement is the vector u(x, t) where x is the position vector and t is time. In terms of angular frequency ω, this can be represented as   $${\bf u}({\bf x},t) = {\bf u}({\bf x},\omega ) \, e^{i \omega t} .$$ (4)The problem will be solved in the frequency domain, so the objective is to determine u(x, ω). In a source-free region conservation of linear momentum requires that   $$\alpha ^2 \nabla \left[ \nabla \cdot {\bf u}({\bf x},\omega ) \right] - \beta ^2 \nabla \times \left[ \nabla \times {\bf u}({\bf u},\omega ) \right] + \omega ^2 {\bf u}({\bf x},\omega ) = 0 .$$ (5)A general solution of this equation can be written in terms of three vectors (L, M, N)   $${\bf u}({\bf x},\omega ) = a^{(L)} {\bf L}({\bf x},\omega ) + a^{(M)} {\bf M}({\bf x},\omega ) + a^{(N)} {\bf N}({\bf x},\omega ) ,$$ (6)where a(L), a(M) and a(N) are constants. These vectors are related to two scalar potential functions ψ and χ that satisfy Helmholtz equations by   $${\bf L}({\bf x},\omega ) = \frac{1}{k_{\alpha }} \nabla \psi ({\bf x},\omega ), \ \ \left[ \nabla ^2 + k_{\alpha }^2 \right] \psi ({\bf x},\omega ) = 0 ,$$ (7a)  $${\bf M}({\bf x},\omega ) = \nabla \times [\chi ({\bf x},\omega ) {\bf r}], \ \ \left[ \nabla ^2 + k_{\beta }^2 \right] \chi ({\bf x},\omega ) = 0 ,$$ (7b)  $${\bf N}({\bf x},\omega ) = \frac{1}{k_{\beta }} \nabla \times {\bf M}({\bf x},\omega ) = \frac{1}{k_{\beta }} \nabla \times ( \nabla \times [ \chi ({\bf x},\omega ) {\bf r} ] ) ,$$ (7c) with wavenumbers kα = ω/α, kβ = ω/β and r the radial vector. Using this representation of u allows separation of variables to be used on eq. (5) in spheroidal coordinates (Morse & Feshbach 1953), so solutions can be written (Flammer 1957) in terms of angular wave functions eimϕ, angular wave functions Smn(c, η), and radial wave functions $$R_{mn}^{(j)}(c,\xi )$$, where m and n are separation indices. Solutions for the compressional and shear potential functions are   $$\psi _{mn}^{(j)}({\bf x},\omega ) = S_{mn}(c_\alpha ,\eta ) R_{mn}^{(j)}(c_\alpha ,\xi ) \, e^{im\phi } ,$$ (8a)and   $$\chi _{mn}^{(j)}({\bf x},\omega ) = S_{mn}(c_\beta ,\eta ) R_{mn}^{(j)}(c_\beta ,\xi ) \, e^{im\phi } ,$$ (8b) where   $$c_\alpha = \frac{h}{2} k_\alpha , \ \ c_\beta = \frac{h}{2} k_\beta .$$ (9)While the ϕ dependence is identical to the case of spherical coordinates, the other two wave functions satisfy the ordinary differential equations   $$\left[ \frac{\text{d}}{\text{d}\eta } (1-\eta ^2)\frac{\text{d}}{\text{d}\eta } - \frac{m^2}{1-\eta ^2} - c^2 \eta ^2 \right] S_{mn}(c,\eta ) = - \Lambda _{mn}(c) S_{mn}(c,\eta ) ,$$ (10a)  $$\left[ \frac{\text{d}}{\text{d}\xi } (\xi ^2-1)\frac{\text{d}}{\text{d}\xi } - \frac{m^2}{\xi ^2-1} + c^2 \xi ^2 \right] R_{mn}^{(j)}(c,\xi ) = \Lambda _{mn}(c) R_{mn}^{(j)}(c,\xi ) ,$$ (10b) where Λmn(c) is an eigenvalue. The two independent radial wave functions are denoted $$R_{mn}^{(1)}$$ and $$R_{mn}^{(2)}$$, where, similar to spherical Bessel and Neumann functions, the former is convergent at the origin and the latter is divergent. Similar to the spherical Hankel function of the second kind, a fourth radial wave function is defined as   $$R_{mn}^{(4)}(c,\xi ) = R_{mn}^{(1)}(c,\xi ) - i R_{mn}^{(2)}(c,\xi ) ,$$ (11)so that the combination of this with eq. (4) corresponds to an outward propagating wave at large ξ. The spheroidal angular function can be related to associated Legendre functions $$P_n^m(.)$$ by   $$S_{mn}(c,\eta ) = \sum _{r=0,1}^{\infty \, ^{\prime \prime }} d_r^{mn} (c) P_{m+r}^m(\eta ) ,$$ (12)where the $$d_r^{mn}(c)$$ are connection coefficients described by Flammer (1957). The double quotes over the summation sign indicates that when n − m is even the sum starts with r = 0 and increments by 2 so that only even values are included, whereas when n − m is odd it starts with r = 1 and increments by 2 so that only odd values are included. This convention will be used throughout this study. The normalization of these angular functions comes from the equation   $$\int _{-1}^{1} S_{mn}(c,\eta ) S_{ml}(c,\eta ) \, \text{d}\eta = N_{mn}^S(c) \delta _{nl} ,$$ (13)where   $$N_{mn}^S(c) = 2 \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{(2m+r)!}{(2m+2r+1)r!} \left[d_r^{mn}(c)\right]^2 .$$ (14) The stress tensor $$\tau$$ generated by a displacement u is (Morse & Feshbach 1953)   $$\tau \left\lbrace {\bf u} \right\rbrace = \left(\kappa - \frac{2}{3} \mu \right) (\nabla \cdot {\bf u}) {\bf I} + \mu \left[ \nabla {\bf u}+{\bf u}\nabla \right].$$ (15)What is actually needed for the boundary conditions is the traction on a surface of constant ξ having unit normal $$\hat{\bf a}_{\xi }$$ and this is the vector   $${\bf T} \left\lbrace {\bf u} \right\rbrace = \tau \left\lbrace {\bf u} \right\rbrace \cdot \hat{\bf a}_{\xi } .$$ (16)Specific expressions in spheroidal coordinates for the displacement vectors L, M, N and the associated tractions are given in Appendix A. Consider the boundary value problem for a general background material having material properties (ρ, α, β) surrounding an inclusion having material properties (ρ΄, α΄, β΄). The displacement of the scattered field outside the inclusion is   $${\bf u}^{\text{sca}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ a_{mn}^{(L)} {\bf L}_{mn}^{(4)}(c_{\alpha }) + a_{mn}^{(N)} {\bf N}_{mn}^{(4)}(c_{\beta }) + a_{mn}^{(M)} {\bf M}_{mn}^{(4)}(c_{\beta }) \right] ,$$ (17a)that of the field interior to the inclusion is   $${\bf u}^{\text{int}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ b_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha ^{\prime }}) + b_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta ^{\prime }}) + b_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta ^{\prime }}) \right] ,$$ (17b)and that of the incident field is   $${\bf u}^{\text{inc}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ s_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta }) \right] ,$$ (17c) The coefficients $$a_{mn}^{(.)}$$ and $$b_{mn}^{(.)}$$ are determined by the boundary conditions. The incident field is taken to be a plane wave parallel to the xy plane that propagates in a positive z-direction that makes an angle ζ with the rp axis of the spheroid (see Figs 1 and 2). First define three functions of the incident angle ζ as   $$f_{mn}^{(0)}(\zeta ) = \frac{2}{N_{mn}^S(c_\alpha )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c_\alpha ) P_{m+r}^m(\cos \zeta ) ,$$ (18a)  $$f_{mn}^{(1)}(\zeta ) = \frac{2m}{N_{mn}^S(c_\beta )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{d_r^{mn}(c_\beta )}{(r+m)(r+m+1)} \frac{P_{m+r}^m(\cos \zeta )}{\sin \zeta } ,$$ (18b)  $$f_{mn}^{(2)}(\zeta ) = \frac{2}{N_{mn}^S(c_\beta )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{d_r^{mn}(c_\beta )}{(r+m)(r+m+1)} \frac{\text{d}}{\text{d}\zeta } P_{m+r}^m(\cos \zeta ) .$$ (18c) Three different types of incident waves are possible with their properties listed in Table 1. For an incident P wave the non-zero expansion coefficients are   $$s_{mn}^{(L)}(\zeta ) = - i f_{mn}^{(0)}(\zeta ) .$$ (19a)For an incident SV wave the non-zero expansion coefficients are   $$s_{mn}^{(N)}(\zeta ) = - i f_{mn}^{(2)}(\zeta ) , \ \ s_{mn}^{(M)}(\zeta ) = - i f_{mn}^{(1)}(\zeta ) .$$ (19b)For an incident SH wave the non-zero expansion coefficients are   $$s_{mn}^{(N)}(\zeta ) = - f_{mn}^{(1)}(\zeta ) , \ \ s_{mn}^{(M)}(\zeta ) = - f_{mn}^{(2)}(\zeta ) .$$ (19c) Table 1. Types of incident waves.   Velocity  Displacement direction  Displacement plane  P  α  ∥ to direction of propagation  in plane containing ζ  SV  β  ⊥ to direction of propagation  in plane containing ζ  SH  β  ⊥ to direction of propagation  ⊥ to plane containing ζ    Velocity  Displacement direction  Displacement plane  P  α  ∥ to direction of propagation  in plane containing ζ  SV  β  ⊥ to direction of propagation  in plane containing ζ  SH  β  ⊥ to direction of propagation  ⊥ to plane containing ζ  View Large 3 BOUNDARY EQUATIONS The boundary conditions are that both displacement and traction should be continuous on the boundary of the inclusion. Letting the boundary of the inclusion be ξ = ξB the displacement boundary equation is   $${\bf u}^{\text{inc}}(\eta ,\xi _B,\phi ,\omega ) + {\bf u}^{\text{sca}}(\eta ,\xi _B,\phi ,\omega ) = {\bf u}^{\text{int}}(\eta ,\xi _B,\phi ,\omega ) ,$$ (20)and the traction boundary equation is   $${\bf T} \left\lbrace {\bf u}^{\text{inc}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace + {\bf T} \left\lbrace {\bf u}^{\text{sca}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace = {\bf T} \left\lbrace {\bf u}^{\text{int}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace .$$ (21) Written in component form, the first of the displacement equations for the η component is   \begin{eqnarray} \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty & & i^n \Big [ b_{mn}^{(L)} L_{\eta mn}^{(1)}(c_{\alpha ^{\prime }}) + b_{mn}^{(N)} N_{\eta mn}^{(1)}(c_{\beta ^{\prime }}) + b_{mn}^{(M)} M_{\eta mn}^{(1)}(c_{\beta ^{\prime }}) \nonumber \\ & & -\,\, a_{mn}^{(L)} L_{\eta mn}^{(4)}(c_{\alpha }) - a_{mn}^{(N)} N_{\eta mn}^{(4)}(c_{\beta }) - a_{mn}^{(M)} M_{\eta mn}^{(4)}(c_{\beta }) \Big ] \nonumber \\ &=& \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \Big [ s_{mn}^{(L)} L_{\eta mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} N_{\eta mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} M_{\eta mn}^{(1)}(c_{\beta }) \Big ] , \end{eqnarray} (22)with similar equations for the ξ and ϕ components. The boundary equation for the η component of traction is   \begin{eqnarray} \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty & & i^n \Big [ b_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(1)}(c_{\alpha ^{\prime }}) \big \rbrace + b_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(1)}(c_{\beta ^{\prime }}) \big \rbrace + b_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(1)}(c_{\beta ^{\prime }}) \big \rbrace \nonumber \\ & & - a_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(4)}(c_{\alpha }) \big \rbrace - a_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(4)}(c_{\beta }) \big \rbrace - a_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(4)}(c_{\beta }) \big \rbrace \Big ] \nonumber \\ &=& \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \Big [ s_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(1)}(c_{\alpha }) \big \rbrace + s_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(1)}(c_{\beta }) \big \rbrace + s_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(1)}(c_{\beta }) \big \rbrace \Big ] , \end{eqnarray} (23)with similar equations for the ξ and ϕ components. At this point the six boundary equations that have to be satisfied are eq. (22) and similar displacement equations for the ξ and ϕ components together with eq. (23) and similar traction equations for the ξ and ϕ components. First note that the frequency ω, which enters through the variables cα and cβ (eq. 9) is fixed in these equations. Next note that (Appendix A and eq. 8) each of the components $$L_{.mn}^{(j)}$$, $$N_{.mn}^{(j)}$$, $$M_{.mn}^{(j)}$$, $$T_. \big \lbrace {\bf L}_{mn}^{(j)} \big \rbrace$$, $$T_. \big \lbrace {\bf N}_{mn}^{(j)} \big \rbrace$$ and $$T_. \big \lbrace {\bf M}_{mn}^{(j)} \big \rbrace$$ that appear in the boundary equations is a function of the coordinates (η, ξ, ϕ). The equations have to be satisfied for the single boundary value ξ = ξB and for all possible values of η and ϕ. Finally note that there is a double sum on each side of the equality sign in these equations. In the case of a spherical inclusion it is possible to use the orthogonality of the wave functions to remove from the boundary equations the dependence on η and ϕ and also the summations over the separation indexes m and n. In the case of a spheroidal inclusion all of this is also possible, except it is not possible to remove the summation over n and it is necessary to add another index l. This important difference from the spherical problem means that there is coupling between modes having different values of n and the boundary equations have to be solved simultaneously for all values of this index. The manipulations required to put the boundary equations in a form suitable for solution are given in Appendix B. After the application of the process described in Appendix B the boundary eq. (22) for the η component of displacement becomes   \begin{eqnarray} \sum _{n=|m|}^\infty & & b_{mn}^{(L)} E_{\eta mnl}^{(L1)}(\alpha ^{\prime }) + b_{mn}^{(N)} E_{\eta mnl}^{(N1)}(\beta ^{\prime }) + b_{mn}^{(M)} E_{\eta mnl}^{(M1)}(\beta ^{\prime }) \nonumber \\ & & - a_{mn}^{(L)} E_{\eta mnl}^{(L4)}(\alpha ) - a_{mn}^{(N)} E_{\eta mnl}^{(N4)}(\beta ) - a_{mn}^{(M)} E_{\eta mnl}^{(M4)}(\beta ) \nonumber \\ & = &\sum _{n=|m|}^\infty s_{mn}^{(L)} E_{\eta mnl}^{(L1)}(\alpha ) + s_{mn}^{(N)} E_{\eta mnl}^{(N1)}(\beta ) + s_{mn}^{(M)} E_{\eta mnl}^{(M1)}(\beta ) , \end{eqnarray} (24)and eq. (23) for the η component of traction becomes   \begin{eqnarray} \sum _{n=|m|}^\infty & & b_{mn}^{(L)} F_{\eta mnl}^{(L1)}(\alpha ^{\prime }) + b_{mn}^{(N)} F_{\eta mnl}^{(N1)}(\beta ^{\prime }) + b_{mn}^{(M)} F_{\eta mnl}^{(M1)}(\beta ^{\prime }) \nonumber \\ & & - a_{mn}^{(L)} F_{\eta mnl}^{(L4)}(\alpha ) - a_{mn}^{(N)} F_{\eta mnl}^{(N4)}(\beta ) - a_{mn}^{(M)} F_{\eta mnl}^{(M4)}(\beta ) \nonumber \\ & & = \sum _{n=|m|}^\infty s_{mn}^{(L)} F_{\eta mnl}^{(L1)}(\alpha ) + s_{mn}^{(N)} F_{\eta mnl}^{(N1)}(\beta ) + s_{mn}^{(M)} F_{\eta mnl}^{(M1)}(\beta ) , \end{eqnarray} (25)with similar equations for the ξ and ϕ components of both displacement and traction at the boundary. Expressions for $$E_{.mnl}^{(.j)}(v)$$ and $$F_{.mnl}^{(.j)}(v)$$ are given in Appendix C and, while they are functions of the boundary value ξ = ξB, they are not functions of η or ϕ. Also note the index m is fixed in these boundary equations, so a set of equations of this form have to be solved independently for each desired value of m. For given values of m and n the six boundary equations can be written in the matrix form   $${\bf C}_{mnl} \, {\bf w}_{mn} = {\bf S}_{mnl} \, {\bf s}_{mn} = {\bf D}_{mnl} ,$$ (26)where   $${\bf C}_{mnl} = \left[ \begin{array}{cccccc}E_{\eta mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\eta mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\eta mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\eta mnl}^{(L4)}(\alpha ) &\quad -E_{\eta mnl}^{(N4)}(\beta ) &\quad -E_{\eta mnl}^{(M4)}(\beta ) \\ E_{\xi mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\xi mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\xi mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\xi mnl}^{(L4)}(\alpha ) &\quad -E_{\xi mnl}^{(N4)}(\beta ) &\quad -E_{\xi mnl}^{(M4)}(\beta ) \\ E_{\phi mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\phi mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\phi mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\phi mnl}^{(L4)}(\alpha ) &\quad -E_{\phi mnl}^{(N4)}(\beta ) &\quad -E_{\phi mnl}^{(M4)}(\beta ) \\ F_{\eta mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\eta mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\eta mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\eta mnl}^{(L4)}(\alpha ) &\quad -F_{\eta mnl}^{(N4)}(\beta ) &\quad -F_{\eta mnl}^{(M4)}(\beta ) \\ F_{\xi mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\xi mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\xi mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\xi mnl}^{(L4)}(\alpha ) &\quad -F_{\xi mnl}^{(N4)}(\beta ) &\quad -F_{\xi mnl}^{(M4)}(\beta ) \\ F_{\phi mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\phi mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\phi mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\phi mnl}^{(L4)}(\alpha ) &\quad -F_{\phi mnl}^{(N4)}(\beta ) &\quad -F_{\phi mnl}^{(M4)}(\beta ) \\ \end{array} \right] ,$$ (27a)  $${\bf w}_{mn} = \left[ \begin{array}{c}b_{mn}^{(L)} \\ b_{mn}^{(N)} \\ b_{mn}^{(M)} \\ a_{mn}^{(L)} \\ a_{mn}^{(N)} \\ a_{mn}^{(M)} \end{array} \right] ,$$ (27b)  $${\bf S}_{mnl} = \left[ \begin{array}{ccc}E_{\eta mnl}^{(L1)}(\alpha ) &\quad E_{\eta mnl}^{(N1)}(\beta ) &\quad E_{\eta mnl}^{(M1)}(\beta ) \\ E_{\xi mnl}^{(L1)}(\alpha ) &\quad E_{\xi mnl}^{(N1)}(\beta ) &\quad E_{\xi mnl}^{(M1)}(\beta ) \\ E_{\phi mnl}^{(L1)}(\alpha ) &\quad E_{\phi mnl}^{(N1)}(\beta ) &\quad E_{\phi mnl}^{(M1)}(\beta ) \\ F_{\eta mnl}^{(L1)}(\alpha ) &\quad F_{\eta mnl}^{(N1)}(\beta ) &\quad F_{\eta mnl}^{(M1)}(\beta ) \\ F_{\xi mnl}^{(L1)}(\alpha ) &\quad F_{\xi mnl}^{(N1)}(\beta ) &\quad F_{\xi mnl}^{(M1)}(\beta ) \\ F_{\phi mnl}^{(L1)}(\alpha ) &\quad F_{\phi mnl}^{(N1)}(\beta ) &\quad F_{\phi mnl}^{(M1)}(\beta ) \\ \end{array} \right] , \ \ {\bf s}_{mn} = \left[ \begin{array}{c}s_{mn}^{(L)} \\ s_{mn}^{(N)} \\ s_{mn}^{(M)} \end{array} \right] ,$$ (27c)  $${\bf D}_{mnl} = \left[ \begin{array}{c}E_{\eta mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\eta mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\eta mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ E_{\xi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\xi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\xi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ E_{\phi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\phi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\phi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\eta mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\eta mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\eta mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\xi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\xi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\xi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\phi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\phi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\phi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ \end{array} \right] .$$ (27d) Using the results from eq. (19) the source vector smn is selected as one of the following   $$\left[ {\bf s}_{mn}^{P} \ {\bf s}_{mn}^{SV} \ {\bf s}_{mn}^{SH} \right] = \left[ \begin{array}{ccc}-i f_{mn}^{(0)}(\zeta ) &\quad 0 &\quad 0 \\ 0 &\quad -i f_{mn}^{(2)}(\zeta ) &\quad - f_{mn}^{(1)}(\zeta ) \\ 0 &\quad -i f_{mn}^{(1)}(\zeta ) &\quad - f_{mn}^{(2)}(\zeta ) \\ \end{array} \right] .$$ (28) The boundary equations for given values of m and n lead to a system of six equations with six unknowns. However, because of coupling between modes it is necessary to simultaneously solve for all values of n. Assuming that the terms in the equations become negligibly small for n > m + J, the combined system for n = m to n = m + J can be written as   $$\left[ \begin{array}{ccccc}\bf C_{mml} &\quad {\bf C}_{m(m+1)l} &\quad {\bf C}_{m(m+2)l} &\quad \cdots &\quad {\bf C}_{m(m+J)l} \end{array} \right] \, \left[ \begin{array}{l}\bf w_{mm} \\ {\bf w}_{m(m+1)} \\ {\bf w}_{m(m+2)} \\ \vdots \\ {\bf w}_{m(m+J)} \end{array} \right] = \sum _{n=m}^{m+J} {\bf D}_{mnl} .$$ (29)This is a system of 6 equations with 6(J + 1) unknowns. This discrepancy is removed by recognizing that this system of equations holds for any value of l. Letting l vary from l = 0 to l = J leads to a square system of 6(J + 1) equations with 6(J + 1) unknowns:   $$\left[ \begin{array}{ccccc}\bf C_{mm0} &\quad {\bf C}_{m(m+1)0} &\quad {\bf C}_{m(m+2)0} &\quad \cdots &\quad {\bf C}_{m(m+J)0} \\ {\bf C}_{mm1} &\quad {\bf C}_{m(m+1)1} &\quad {\bf C}_{m(m+2)1} &\quad \cdots &\quad {\bf C}_{m(m+J)1} \\ {\bf C}_{mm2} &\quad {\bf C}_{m(m+1)2} &\quad {\bf C}_{m(m+2)2} &\quad \cdots &\quad {\bf C}_{m(m+J)2} \\ \vdots &\quad \vdots &\quad \vdots &\quad \cdots &\quad \vdots \\ {\bf C}_{mmJ} &\quad {\bf C}_{m(m+1)J} &\quad {\bf C}_{m(m+2)J} &\quad \cdots &\quad {\bf C}_{m(m+J)J} \\ \end{array} \right] \, \left[ \begin{array}{l}\bf w_{mm} \\ {\bf w}_{m(m+1)} \\ {\bf w}_{m(m+2)} \\ \vdots \\ {\bf w}_{m(m+J)} \end{array} \right] = \left[ \begin{array}{l}\sum _{n=m}^{m+J} {\bf D}_{mn0} \\ \sum _{n=m}^{m+J} {\bf D}_{mn1} \\ \sum _{n=m}^{m+J} {\bf D}_{mn2} \\ \vdots \\ \sum _{n=m}^{m+J} {\bf D}_{mnJ} \\ \end{array} \right] .$$ (30) In the special case of m = 0 several terms of the coefficient matrix Cmnl are identically zero and the 6 × 6 system of eq. (26) splits into two independent systems, a 4 × 4 system involving the L and N potentials and a 2 × 2 system involving the M potential. For the 4 × 4 system the individual matrices are   $${\bf C}_{0nl} = \left[ \begin{array}{cccc}E_{\eta 0nl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\eta 0nl}^{(N1)}(\beta ^{\prime }) &\quad - E_{\eta 0nl}^{(L4)}(\alpha ) &\quad -E_{\eta 0nl}^{(N4)}(\beta ) \\ E_{\xi 0nl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\xi 0nl}^{(N1)}(\beta ^{\prime }) &\quad - E_{\xi 0nl}^{(L4)}(\alpha ) &\quad -E_{\xi 0nl}^{(N4)}(\beta ) \\ F_{\eta 0nl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\eta 0nl}^{(N1)}(\beta ^{\prime }) &\quad - F_{\eta 0nl}^{(L4)}(\alpha ) &\quad -F_{\eta 0nl}^{(N4)}(\beta ) \\ F_{\xi 0nl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\xi 0nl}^{(N1)}(\beta ^{\prime }) &\quad - F_{\xi 0nl}^{(L4)}(\alpha ) &\quad -F_{\xi 0nl}^{(N4)}(\beta ) \\ \end{array} \right] ,$$ (31a)  $${\bf w}_{0n} = \left[ \begin{array}{c}b_{0n}^{(L)} \\ b_{0n}^{(N)} \\ a_{0n}^{(L)} \\ a_{0n}^{(N)} \end{array} \right] ,$$ (31b)  $${\bf D}_{0nl} = \left[ \begin{array}{c}E_{\eta 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + E_{\eta 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ E_{\xi 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + E_{\xi 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ F_{\eta 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + F_{\eta 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ F_{\xi 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + F_{\xi 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ \end{array} \right] ,$$ (31c) and for the 2 × 2 system they are   $${\bf C}_{0nl} = \left[ \begin{array}{cc}E_{\phi 0nl}^{(M1)}(\beta ^{\prime }) &\quad -E_{\phi 0nl}^{(M4)}(\beta ) \\ F_{\phi 0nl}^{(M1)}(\beta ^{\prime }) &\quad -F_{\phi 0nl}^{(M4)}(\beta ) \\ \end{array} \right] ,$$ (32a)  $${\bf w}_{0n} = \left[ \begin{array}{c}b_{0n}^{(M)} \\ a_{0n}^{(M)} \end{array} \right] ,$$ (32b)  $${\bf D}_{0nl} = \left[ \begin{array}{c}E_{\phi 0nl}^{(M1)}(\beta ) s_{0n}^{(M)} \\ F_{\phi 0nl}^{(M1)}(\beta ) s_{0n}^{(M)} \\ \end{array} \right] .$$ (32c) Although not necessary, the dimension of this system of equations can be further reduced by taking advantage of the fact pointed out in Appendix C that many of the terms in the coefficient matrices of eqs (27a), (31a) and (32a) are zero, depending upon whether n − m + l is even or odd. For m ≠ 0, this means that the 6(J + 1) × 6(J + 1) general system can be split into two 3(J + 1) × 3(J + 1) systems. Let p = 0, 1, 2, 3, 4, ⋅⋅⋅ be an integer. Then the first system involves the coefficients   $$b_{m,m+2p+1}^{(L)} , b_{m,m+2p+1}^{(N)} , b_{m,m+2p}^{(M)} , a_{m,m+2p+1}^{(L)} , a_{m,m+2p+1}^{(N)} , a_{m,m+2p}^{(M)} ,$$ (33a)while the second involves   $$b_{m,m+2p}^{(L)} , b_{m,m+2p}^{(N)} , b_{m,m+2p+1}^{(M)} , a_{m,m+2p}^{(L)} , a_{m,m+2p}^{(N)} , a_{m,m+2p+1}^{(M)} .$$ (33b) For m = 0 the general system can be split into four systems. The first is a 2(J + 1) × 2(J + 1) system that involves   $$b_{0,2p}^{(L)} , b_{0,2p}^{(N)} , a_{0,2p}^{(L)} , a_{0,2p}^{(N)} .$$ (34a)The second is a 2(J + 1) × 2(J + 1) system that involves   $$b_{0,2p+1}^{(L)} , b_{0,2p+1}^{(N)} , a_{0,2p+1}^{(L)} , a_{0,2p+1}^{(N)} .$$ (34b)The third is a (J + 1) × (J + 1) system that involves   $$b_{0,2p}^{(M)} , a_{0,2p}^{(M)} .$$ (34c)The fourth is a (J + 1) × (J + 1) system that involves   $$b_{0,2p+1}^{(M)} , a_{0,2p+1}^{(M)} .$$ (34d) What is clear from these results is that for any of the expansion coefficients $$b_{mn}^{(.)}$$ or $$a_{mn}^{(.)}$$ there is no coupling between coefficients having different values of m. However, there is coupling between all coefficients having even values of n and between all coefficients having odd values of n, but no coupling between even and odd values. Furthermore, for m ≠ 0 there is coupling between coefficients having superscripts L, N and M. For m = 0 there is coupling between coefficients having superscripts L and N, but these are decoupled from coefficients having superscript M. The solution to the scattering problem is now complete. Given the geometry of the problem and the material properties of the inclusion and background, it is possible to calculate the spheroidal wave functions Smn(cv, η) and $$R_{mn}^{(j)}(c_v,\xi )$$. The requirement of continuity of displacement and traction on the boundary of the inclusion leads to equations in the form of (22) and (23), and applying the methods of Appendix B to these equations results in the boundary coefficients $$E_{.mnl}^{(.j)}$$ and $$F_{.mnl}^{(.j)}$$ (listed in Appendix C) that comprise the coefficient matrix of eq. (30). Selection of a source type from eq. (28) determines the right hand side of eq. (30). Solving the square linear system of equations in eq. (30) leads to the solution wmn, and then eq. (17) determines the part of the displacement due to scattering both outside and inside the inclusion. Note that this solution process is repeated for each value of the index m and each value of the frequency ω. The preceding development is for the prolate spheroidal coordinate system, but, as pointed out by Flammer (1957) and Asano & Yamamoto (1975), all of the results can be converted to the oblate spheroidal system through the substitutions   $$\xi \rightarrow i \xi , \ \ c \rightarrow - i c .$$ (35) 4 NUMERICAL CONSIDERATIONS In order to solve the boundary equations and obtain numerical values for the scattering coefficients it is necessary to calculate the eigenvalues Λmn(c), the expansion coefficients $$d_r^{mn}(c)$$, the angular functions Smn(c, η) and the radial functions $$R_{mn}^{(j)}(c,\xi )$$. In these calculations it is possible to take advantage of numerous resources available in the literature that were developed primarily for the scattering of electromagnetic waves. Morse & Feshbach (1953), Erdélyi et al. (1955), Stratton et al. (1956), Flammer (1957) and Abramowitz & Stegun (1964) all provide useful discussions of functions and algorithms that can be used in the calculations. Additional discussions of this type along with computer programs can be found in Zhang & Jin (1996), Thompson (1997), Li et al. (2002) and Falloon et al. (2003). More detailed discussions of specific aspects of the calculations can be found in Hodge (1970), Sinha et al. (1973), Do-Nhat & MacPhie (1996, 1997) and Do-Nhat (1999, 2001). While in theory the same algorithm can be used for both prolate and oblate calculations with the insertion of an i in the appropriate places (eq. 35), additional modifications may be necessary in order to obtain accurate results. The calculation of the spheroidal radial functions is particularly difficult, with different algorithms generally used for prolate and oblate functions and for large and small values of cξ. A practical problem involves the fact that the limits of the sums in the representations of incident and scattered fields in eq. (17) are infinite. In numerical calculations these limits have to be replaced by finite numbers, so the incident field (eq. 17c) becomes   $${\bf u}^{\text{inc}} = \sum _{m=-m_{\text{max}}}^{m_{\text{max}}} \sum _{n=|m|}^{n_{\text{max}}} i^n \left[ s_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta }) \right] ,$$ (36)with similar modifications to the scattered and interior fields. This same problem is discussed by Korneev & Johnson (1993a) and a similar approach is followed here. The basic goal is to choose mmax and nmax large enough so that convergence of the series is achieved to some desired precision, but at the same time not so large that the terms in the series cannot be calculated to this same precision. Experience has shown that this goal is best obtained by considering the incident field, as this field converges more slowly than the scattered and interior fields. The calculations in the present study use   $$m_{\text{max}} = 10 + \frac{20}{\pi } \zeta + k_{\beta } r_s ,$$ (37a)  $$n_{\text{max}} = m + 10 + \frac{1}{2} \Big ( e + \frac{h}{r_s} \Big ) k_{\beta } r_s,$$ (37b) where ζ is given in units of radians and $$r_s = (r_p r_e^2)^{1/3}$$ is the spherical radius for a volume equal to that of the spheroid. Note that in the special case where ζ = 0, mmax = 1. The application of the boundary conditions leads to a set of linear equations in the form of eq. (30) that must be solved to obtain the solution coefficients wmm...wm(m + J). Formally this amounts to inverting a square matrix of dimension 6(J + 1) in the case where m > 0 and a smaller dimension when m = 0. In practice this is best achieved with an LU factorization followed by back substitution, which is performed using subroutines from LAPACK in the present study. 5 NUMERICAL RESULTS This section presents a few sample calculations that illustrate some of the properties of the elastic waves scattered by a spheroidal inclusion. The calculations are similar to some of those presented by Korneev & Johnson (1993a, 1996) for elastic waves scattered by a spherical inclusion and make use of the geometry shown in Figs 1 and 2. The model of material properties used for the calculations is listed in Table 2 and is a spheroidal inclusion that has velocities and density about 25 per cent less than the surrounding background material. This model is identical to model 2 of Korneev & Johnson (1996). Table 2. Model used for calculations.   Density  Compressional velocity  Shear velocity  Background  ρ = 2700 kgm−3  α = 6000 ms−1  β = 3500 ms−1  Inclusion  ρ΄ = 2300 kgm−3  α΄ = 4500 ms−1  β΄ = 2600 ms−1    Density  Compressional velocity  Shear velocity  Background  ρ = 2700 kgm−3  α = 6000 ms−1  β = 3500 ms−1  Inclusion  ρ΄ = 2300 kgm−3  α΄ = 4500 ms−1  β΄ = 2600 ms−1  View Large As shown in Figs 1 and 2, a plane wave is incident upon the spheroidal inclusion and solutions for the three components of the vector displacement usca(x, t) associated with the scattered field are calculated on a plane displaced a distance zo from the centre of the spheroid. The coordinate system (x, y, z) shown in these figures is rotated by an angle ζ with respect to the (x1, x2, x3) system of the spheroid given in eq. (3). Unlike the case of a spherical inclusion, the angle ζ between the direction of propagation of the incident wave and the polar axis of the spheroid is now an important parameter. It is necessary to distinguish three different types of incident waves and these were listed above in Table 1. In the geometry shown in Figs 1 and 2 the incident P wave has displacement in the z-direction, the incident SV wave in the x-direction, and the incident SH wave in the y-direction. In the case where ζ = 0 there is no distinction between the SV and SH incident waves, which is similar to the case of a spherical inclusion. Both prolate and oblate spheroidal inclusions are considered with both having an aspect ratio of 2 and both having a volume equal to that of a sphere having a radius of 1 m. Thus for the prolate inclusion   $$r_p = 1.588 \, m , \ \ r_e = 0.794 \, m ,$$ (38a)and for the oblate spheroid   $$r_p = 0.630 \, m , \ \ r_e = 1.260 \, m .$$ (38b)These dimensions give the boundary of the spheroid as ξB = 1.15 for the prolate case and ξB = 0.58 for the oblate case. The angle of incidence ζ and offset distance zo used for the calculations are   $$\zeta = 30^{\circ} , \ \ z_o = 2 \, m .$$ (38c) One method of characterizing the scattering of waves by an inclusion is to calculate and display scattering cross-sections as a function of frequency. The approach of Korneev & Johnson (1993a) is used to derive normalized total scattering cross-sections for the scattering of elastic waves by a spheroidal inclusion. Such cross-sections for incident P, SV and SH waves are   $$\sigma _{P} = - \frac{4\pi }{A_g(\zeta ) k_{\alpha }^2} \Im \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } a_{mn}^{(L)} W_{mn}^{(1)}(\zeta ) \right] ,$$ (39a)  $$\sigma _{SV} = \frac{4\pi }{A_g(\zeta ) k_{\beta }^2} \Im \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } \left( a_{mn}^{(N)} W_{mn}^{(3)}(\zeta ) + a_{mn}^{(M)} W_{mn}^{(2)}(\zeta ) \right) \right],$$ (39b)  $$\sigma _{SH} = - \frac{4\pi }{A_g(\zeta ) k_{\beta }^2} \Re \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } \left( a_{mn}^{(M)} W_{mn}^{(3)}(\zeta ) + a_{mn}^{(N)} W_{mn}^{(2)}(\zeta ) \right) \right] ,$$ (39c) where   $$W_{mn}^{(1)}(\zeta ) = P_n^m(\cos \zeta ) ,$$ (40a)  $$W_{mn}^{(2)}(\zeta ) = m \frac{P_n^m(\cos \zeta )}{\sin \zeta } ,$$ (40b)  $$W_{mn}^{(3)}(\zeta ) = \frac{d}{d \zeta } P_n^m(\cos \zeta ) ,$$ (40c) and the geometrical shadow of the inclusion is   $$A_g (\zeta ) = \pi r_e^2 \left[ 1 + \left( (r_p/r_e)^2 - 1 \right) \sin ^2(\zeta ) \right]^{1/2} .$$ (41)This geometrical shadow is the area of the plane incident wave that intersects the inclusion and is a function of the angle of incidence ζ. Normalizing the scattering cross-sections by this area allows them to be interpreted as the ratio of total energy scattered by the inclusion to total energy incident on the inclusion. The scattering cross-sections for incident P, SV and SH waves are plotted in Figs 3– 5, respectively, for the dimensions and angle of incidence listed in eq. (38). Each plot contains results for a prolate and oblate spheroid, both having an aspect ratio of 2, and a sphere, with all three inclusions having the same volume. The abscissa in each plot is a normalized frequency that is the wavenumber of the background compressional wave kα times the radius of a sphere having the same volume as the inclusion, which is rs = 1 m for Figs 3–5. This wavenumber is that of the incident wave in Fig. 3 but not in Figs 4 and 5 so that the three scattering cross-sections can be compared with the same horizontal scale. Figure 3. View largeDownload slide Normalized total scattering cross-sections for a plane P wave incident at an angle of 30° with respect to the polar axis of the inclusion. The solid line is for a sphere, the long dash line for a prolate spheroid and the short dash line for an oblate spheroid, where all three inclusions have the same volume of a sphere with a 1 m radius. Figure 3. View largeDownload slide Normalized total scattering cross-sections for a plane P wave incident at an angle of 30° with respect to the polar axis of the inclusion. The solid line is for a sphere, the long dash line for a prolate spheroid and the short dash line for an oblate spheroid, where all three inclusions have the same volume of a sphere with a 1 m radius. Figure 4. View largeDownload slide Similar to Fig. 3 for an incident plane SV wave. Figure 4. View largeDownload slide Similar to Fig. 3 for an incident plane SV wave. Figure 5. View largeDownload slide Similar to Fig. 3 for an incident plane SH wave. Figure 5. View largeDownload slide Similar to Fig. 3 for an incident plane SH wave. The scattering cross-sections in Figs 3–5 are all similar in that they exhibit characteristic Rayleigh type scattering at low frequencies where the amplitude is proportional to ω4 (Korneev & Johnson 1996). However, in this range the amplitudes for both the prolate and oblate spheroids are always less than for the sphere, regardless of the type of incident wave. This low frequency behaviour merges into large long oscillations that eventually tend to a constant value of 2.0 at high frequencies. Superimposed on these large long oscillations are numerous small short oscillations that are referred to as ripple. The basic features of the scattering cross-sections in Figs 3–5 are qualitatively similar to scattering cross-sections for electromagnetic waves incident upon a spherical inclusion, so interpretations developed for the electromagnetic problem can be borrowed for the elastic problem. Consider the scattering cross-section for incident P waves shown in Fig. 3. van de Hulst (1957) explains the large long oscillations as an interference pattern between waves transmitted through the inclusion and waves diffracted around it for electromagnetic waves. This interpretation applied to the elastic problem does reasonably well in predicting the peaks in the large long oscillations of the scattering cross-section for the sphere in Fig. 3 and can be extended to the cross-sections for the spheroids. The basic idea is that when waves passing through the inclusion are in-phase with those traveling around the inclusion there will be constructive interference and peaks in scattering cross-sections. Similarly, when these two types of waves are out-of-phase there will be destructive interference and troughs in the scattering cross-sections. The phase difference between waves passing through and around the inclusion depends upon frequency, velocity contrast of the inclusion, and dimension of the inclusion parallel to the direction of the incident wave. Noting that the projection of the inclusion onto the z-axis is 2.00 m for the sphere, 2.86 m for the prolate spheroid, and 1.67 m for the oblate spheroid (see Figs 1 and 2), the interpretation of van de Hulst (1957) predicts that the peaks in the large oscillations should be shifted to lower frequencies than the sphere for the prolate spheroid and to higher frequencies for the oblate spheroid, which agrees with the plots in Fig. 3. Also note that the first peak of the cross-section is larger for the prolate spheroid compared to the sphere by a factor of 1.08 and smaller for the oblate spheroid by a factor of 0.90. For cross-sections of electromagnetic waves incident upon a spherical inclusion Bohren & Huffman (1983) give an interpretation for the ripple in terms of the resonant modes of the inclusion. It seems likely that such an interpretation could apply to the case of elastic waves incident upon a spherical inclusion and also to the cases of prolate and oblate spheroids. It is important to note that the appearance of the ripple is influenced by the frequency sampling rate used in the calculations. The scattering cross-sections for incident SV and SH waves are shown in Figs 4 and 5, respectively. The general appearance of the cross-sections are similar to those of P waves but with all features shifted to lower frequencies. The velocities of S waves within the inclusion are about 0.7 times those of P waves. Again applying the interpretation of van de Hulst (1957), one would expect the peaks in the large oscillations to be shifted to lower frequencies by about a factor of about 0.7, and this is approximately the case. Similar to the case of an incident P wave, the first peak in the cross-section is shifted to lower frequencies for the prolate spheroid compared to the sphere and to higher frequencies for the oblate spheroid. In all cases the amplitude of the first peak of the scattering cross-section is larger for incident S waves than for incident P waves. For the prolate case the ratio of SV amplitude to P amplitude is 1.08 and the ratio for SH amplitude to P amplitude is 1.06. For the oblate case these ratios are both 1.01 The differences between the scattering cross-sections for an incident SV wave in Fig. 4 and an incident SH wave in Fig. 5 are small. This just reflects the fact that effects related to the polarization of the incident wave that might appear in differential cross-sections get averaged out when the integration is performed to obtain the total cross-sections. The scattering cross-sections in Figs 3–5 are for an angle of incidence ζ = 30° and an aspect ratio Ra = 2. A complete investigation of variations in these parameters as well as the velocity contrast of the inclusion is beyond the scope of the present study, but some of the effects caused by such variations can be anticipated. Asano (1979) considers a range in velocity contrasts, aspect ratios, and angles of incidence for the electromagnetic problem, and it appears that some of these results can be extended to the elastic problem, at least in a qualitative sense. As a check on this possibility consider the effect on the scattering cross-sections of Figs 3–5 when the angle of incidence is changed from ζ = 30° to ζ = 0°. The results for the sphere do not change but the peaks in the scattering cross sections shift to lower frequencies for the prolate spheroid and to higher frequencies for the oblate spheroid, which is consistent with the results of Asano (1979). These results are also consistent with the van de Hulst (1957) interpretation in terms of an interference pattern, as the projection of the inclusion onto the direction of the incident wave has now changed from 2.86 to 3.17 m for the prolate spheroid, suggesting a shift to lower frequencies, and from 1.67 to 1.26 m for the oblate spheroid, suggesting a shift to higher frequencies. This change in angle of incidence also causes the amplitude of the peaks to increase for the prolate spheroid and decrease for the oblate spheroid, which is also consistent with the findings of Asano (1979). The amount of ripple shows little change for this change in the angle of incidence. Given the general similarities between scattering cross-sections for electromagnetic and elastic waves, some of the other observations of Asano (1979) might also apply to scattering of elastic waves by spheroidal inclusions. A good test of any scattering solution is the ability to generate waveforms of the scattered fields, as this requires both modulus and phase to be accurately calculated over a broad frequency range. Figs 6– 15 display waveforms of scattered waves that are generated exterior to the inclusion for P, SV and SH types of incident waves. The waveforms are recorded along a profile that extends along the x-axis at a distance zo from the centre of the inclusion with the y coordinate set to zero (see Figs 1 and 2). The incident wave is a plane wave propagating in a direction that makes an angle ζ = 30° with the rp axis of the spheroidal inclusion, which has an aspect ratio of 2. Figs 6–10 are for a prolate spheroid and Figs 11– 15 are for an oblate spheroid. The calculations are performed with a Nyquist frequency of 12.5 KHz, and the incident wave is a delta function in time that is passed through a low-pass causal filter having a corner frequency of 6.25 KHz and a high frequency roll-off of 4. Figure 6. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 6. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 7. View largeDownload slide Similar to Fig. 6 for the scattered field with displacements in the z-direction. Figure 7. View largeDownload slide Similar to Fig. 6 for the scattered field with displacements in the z-direction. Figure 8. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 8. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 9. View largeDownload slide Similar to Fig. 8 for the scattered field with displacements in the z-direction. Figure 9. View largeDownload slide Similar to Fig. 8 for the scattered field with displacements in the z-direction. Figure 10. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 10. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 11. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 11. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 12. View largeDownload slide Similar to Fig. 11 for the scattered field with displacements in the z-direction. Figure 12. View largeDownload slide Similar to Fig. 11 for the scattered field with displacements in the z-direction. Figure 13. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 13. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 14. View largeDownload slide Similar to Fig. 13 for the scattered field with displacements in the z-direction. Figure 14. View largeDownload slide Similar to Fig. 13 for the scattered field with displacements in the z-direction. Figure 15. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 15. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. In viewing the waveforms in Figs 6–9 and Figs 11–14 it is important to keep in mind that the scattering process causes conversions between P and SV waves and in general the converted waves are not propagating parallel to the incident wave. Thus in plots with either of these waves as the incident wave, it is not always obvious whether a particular waveform represents a P wave, an SV wave, or a combination of the two. Figs 6 and 7 show waveforms with displacements in the x and z directions, respectively, for a plane P wave incident on a prolate spheroid with velocities lower than the surrounding background material (Table 2). Only the scattered field outside the inclusion is shown, so all of these waveforms would have zero amplitude if the inclusion were not present. Because of the lower velocity within the inclusion, the waves passing through the inclusion become focused and have a fairly complicated pattern of triplications in the traveltime curves, caustics, and diffractions. The scattered waves extend well beyond the geometrical shadow of the inclusion, which has limits of ±1.05 m along the x coordinate in this case. For x positions inside the geometrical shadow the z component of displacement, which is that of the incident wave, is larger than the x component, but for x positions outside the geometrical shadow this situation can be reversed. Because of the tilt of the inclusion with ζ = 30° the waveforms are not symmetrical with respect to the x coordinate, with larger amplitudes and more extended durations for negative values of the x coordinate where the inclusion is in closer proximity to the observation point (see Fig. 1). For comparison, note that in the case of a spherical inclusion or a spheroidal inclusion with ζ = 0 the waveforms are symmetrical with respect to the x and y coordinates of the observation point. Figs 8 and 9 are similar to Figs 6 and 7 except that the incident wave is a plane SV wave with displacement in the x-direction. In many respects the waveforms are similar to the case of the incident P wave except the roles of the x and z displacement are reversed and the first arriving waves are delayed in time. The asymmetry with respect to the x-axis is now somewhat larger. Also note that in Fig. 9 prior to the time of the incident SV wave there are arrivals caused by P waves converted from the incident wave by the scattering process. The waveforms generated by an incident SH wave on a prolate inclusion are shown in Fig. 10 for a profile along the x-axis. The waveforms are for displacement in the y-direction, the same direction as the displacement of the incident wave. Because of the mirror symmetry of the inclusion about the xz plane, the displacements of the scattered waves in the x and z-directions are zero along this profile. The waveforms for this incident SH wave are quite similar in shape to those of the incident SV wave in Fig. 8, but the amplitudes are generally larger. This is because part of the energy of the incident SV wave is converted to displacement in the z-direction, as shown in Fig. 9. The waveforms in Figs 11–15 contain the same incident waves as those of Figs 6–10 except the inclusion is now oblate instead of prolate (compare Fig. 2 to Fig. 1). The limits of the geometrical shadow are now ±1.14 m along the x-axis. While there is a general similarity of the scattered waves for the prolate and oblate inclusions, there are also some specific differences. The amplitudes of the displacements for the oblate case are generally less than for the prolate case, which is consistent with the lower peak values of the scattering cross-sections in Figs 3–5. An exception is the y component of displacement for an incident SH wave in Fig. 15 where strong focussing causes amplitudes that are more than 3 times the amplitude of the incident wave. The waveforms are still asymmetrical with respect to the x coordinate, but now the larger amplitudes occur for positive values of the x coordinate where the inclusion is in closer proximity to the observation point (see Fig. 2). Fig. 14 is similar to Fig. 9 in that it shows converted P waves arriving before the incident SV wave. It is important to point out that the plots of the scattered waveforms in Figs 6–15 are only part of the realizable field that consists of the sum of the incident field and scattered field. In Fig. 16, the scattered field for the incident plane SV wave shown in Fig. 8 is combined with the incident field for that case to give the total field. Comparing Fig. 8 with Fig. 16 demonstrates the fact that one of the functions of the scattered field is to cancel out part of the incident field, as at small distances and times extending beyond the arrival time of the incident wave there is a shadow zone. Also note that outside the geometrical shadow of the inclusion the incident wave has a reduced amplitude, so the effect of the inclusion extends well beyond its geometrical shadow. Figure 16. View largeDownload slide Waveforms of the total field that consists of the scattered field of Fig. 8 plus the incident field for displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 16. View largeDownload slide Waveforms of the total field that consists of the scattered field of Fig. 8 plus the incident field for displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. The calculations of this section used a maximum value for the wavenumber in the background material of kβ = 22 m−1. With ζ = 30° this required mmax = 35 (eq. 37a) and nmax = 107 (eq. 37b). The corresponding maximums in the spheroidal parameter c are cβ = 31 for the prolate case and cβ = 24 for the oblate case. When either the maximum frequency or aspect ratio of the spheroid is increased so that larger values of c are required, problems are encountered with the numerical accuracy of the radial wave functions and alternative methods of calculation may be required (Do-Nhat 1999, 2001; Kirby 2006, 2010). 6 DISCUSSION AND CONCLUSIONS One method of characterizing the solution for the scattering of elastic waves by a spheroidal inclusion is to compare it to the corresponding solution for a sphere. The spheroidal problem transforms to the spherical problem in the limit as   $$r_p \rightarrow r_e = r_s , \ h = 2 \sqrt{ r_p^2 - r_e^2} \rightarrow 0 , \ c_v = \frac{h}{2} k_v \rightarrow 0 ,$$ (42a)where v can be either the compressional velocity α or shear velocity β and kv = ω/v is a wavenumber. In this limit the spheroidal coordinates transform as   $$\eta \rightarrow \cos \theta , \ \xi \rightarrow \infty , \ \phi \rightarrow \phi ,$$ (42b)with the condition that   $$c_v \xi \rightarrow k_v r .$$ (42c)The rather complicated spheroidal eigenvalues transform to the simple form   $$\Lambda _{mn} (c_v) \rightarrow n (n+1) .$$ (42d)The angular spheroidal wave functions transform to associated Legendre functions   $$S_{mn} (c_v,\eta ) \rightarrow P_n^m (\cos \theta ) ,$$ (42e)and the radial spheroidal wave functions transform to spherical Bessel and Hankel functions   $$R_{mn}^{(1)} (c_v,\xi ) \rightarrow j_n (k_v r) , \ R_{mn}^{(2)} (c_v,\xi ) \rightarrow y_n (k_v r) , \ R_{mn}^{(4)} (c_v,\xi ) \rightarrow k_n^{(2)} (k_v r) .$$ (42f)In the general case ζ ≠ 0 the range of the separation index m is greatly reduced   $$0 \le m \le \infty \rightarrow m = 0, 1 ,$$ (42g)while the range of the separation index n is the same in the spheroidal and spherical problems. An obvious difference between the spheroidal and spherical problem is that the spheroidal problem actually consists of two different types of inclusions, prolate and oblate, with different computational methods leading to rather different scattered fields for the two types. This difference also represents a major asset of the solution to the spheroidal problem, the ability to model inclusions having a large variety in both shape and aspect ratio. Another very important difference between the spheroidal and spherical problems in a computational sense is the fact that the boundary equations remain coupled with respect to the separation index n in the case of a spheroid but are fully uncoupled in the case of a sphere. A few sample calculations are presented to illustrate some of the properties of the solution. Normalized total scattering cross-sections for incident P, SV and SH waves are compared for spherical, prolate, and oblate inclusions having velocities less than the surrounding background material. While the general features of the cross-sections are similar for these three cases, there are also some systematic differences in the amplitudes and frequency dependences. The interpretation of the scattering cross-sections can be guided by methods that have been developed for scattering cross sections of electromagnetic waves. Waveforms are calculated for the scattered field exterior to the inclusion for obliquely incident P, SV and SH waves. These calculations illustrate the importance of the angle ζ between the direction of propagation of the incident wave and the polar axis of the spheroid. In the case where ζ ≠ 0 some features of the spheroidal problem not found in the spherical problem include scattered fields that change as ζ is changed, a difference in the scattered fields for an incident SV wave and an incident SH wave, and scattered fields that do not have symmetry about the centre of the spheroid. Other more general scattering effects contained in the calculated waveforms are scattered fields that extend well outside the geometrical shadow of the inclusion, strong focussing for inclusions with velocities less than the background, and conversions between P and SV waves. While the solution presented in this study is complete in the sense that it is valid for any prolate or oblate spheroid, the application of this solution still has limits imposed by current capabilities to efficiently calculate spheroidal eigenvalues and wave functions for arbitrary combinations of material properties, aspect ratios, and frequency. Further progress on the numerical aspects of these calculations could help extend the utility of the spheroidal solution. Given the exact analytical solution presented in this study, it should be possible to test various approximate solutions such as those mentioned in the introduction, or possibly develop new approximate solutions (see Korneev & Johnson 1993b, for such approximations in the case of a spherical inclusion). For certain selected ranges of the variables and parameters, such approximate solutions may be completely acceptable and more computationally feasible than the exact solution. ACKNOWLEDGEMENTS This research is dedicated to the memory of Valeri Korneev and his contributions to problems involving the scattering of elastic waves by a spherical inclusion. REFERENCES Abramowitz M., Stegun I.A. (Eds), 1964. Handbook of Mathematical Functions  AMS 55, National Bureau of Standards, 1046 pp. Asano S., 1979. Light scattering properties of spheroidal particles, Appl. Opt.  18 712– 723. Google Scholar CrossRef Search ADS PubMed  Asano S., Yamamoto G., 1975. Light scattering by a spheroidal particle, Appl. Opt.  14 29– 49. Google Scholar CrossRef Search ADS PubMed  Bohren C.F., Huffman D.R., 1983. Absorption and Scattering of Light by Small Particles  Wiley, 530 pp. Dassios G., Kleinman R., 2000. Low Frequency Scattering  Clarendon Press, 297 pp. Datta S.K., 1977. Diffraction of plane elastic waves by ellipsoidal inclusions, J. acoust. Soc. Am.  61 1432– 1437. Google Scholar CrossRef Search ADS   Do-Nhat T., 1999. Asymptotic expansions of the Mathieu and prolate spheroidal eigenvalues for large parameter c, Can. J. Phys.  77 635– 652. Google Scholar CrossRef Search ADS   Do-Nhat T., 2001. Asymptotic expansions of the oblate spheroidal eigenvalues and wave functions for large parameter c, Can. J. Phys.  79 813– 831. Google Scholar CrossRef Search ADS   Do-Nhat T., MacPhie R.H., 1996. On the accurate computation of prolate spheroidal radial functions of the second kind, Q. Appl. Math.  55 677– 685. Google Scholar CrossRef Search ADS   Do-Nhat T., MacPhie R.H., 1997. Accurate values of prolate spheroidal radial functions of the second kind, Can. J. Phys.  75 671– 675. Google Scholar CrossRef Search ADS   Erdélyi A. (ed), Magnus W., Oberhettinger F., Tricomi F.G., 1953. Higher Transcendental Functions  Vol. 1, McGraw-Hill, 302 pp. Erdélyi A. (ed), Magnus W., Oberhettinger F., Tricomi F. G., 1955. Higher Transcendental Functions  Vol. 3, McGraw-Hill, 292 pp. Falloon P.E., Abbot P.C., Wang J.B., 2003. Theory and computation of spheroidal wave functions, J. Phys. A: Math. Gen.  36 5477– 5495. Google Scholar CrossRef Search ADS   Flammer C., 1957. Spheroidal Wave Functions  Stanford University Press, 220 pp. Gubernatis J.E., Domany E., Krumhansl J.A., 1977a. Formal aspects of the theory of the scattering of ultrasound by flaws in elastic materials, J. Appl. Phys.  48 2804– 2811. Google Scholar CrossRef Search ADS   Gubernatis J.E., Domany E., Krumhansl J.A., 1977b. The Born approximation in the theory of the scattering of elastic waves by flaws, J. Appl. Phys.  48 2812– 2819. Google Scholar CrossRef Search ADS   Hackman R.H., 1984. The transition matrix for acoustic and elastic wave scattering in prolate spheroidal coordinates, J. acoust. Soc. Am.  75 35– 45. Google Scholar CrossRef Search ADS   Hackman R.H., Lim R., 1994. Development and application of the spheroidal coordinate based T matrix solution to elastic wave scattering, Radio Sci.  29 1035– 1049. Google Scholar CrossRef Search ADS   Hodge D.P., 1970. Eigenvalues and eigenfunctions of the spheroidal wave equation, J. Math. Phys.  11 2308– 2312. Google Scholar CrossRef Search ADS   Kirby P., 2006. Calculation of spheroidal wave functions, Comput. Phys. Commun.  175 465– 472. Google Scholar CrossRef Search ADS   Kirby P., 2010. Calculation of radial prolate spheroidal wave functions of the second kind, Comput. Phys. Commun.  181 514– 519. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1993a. Scattering of elastic waves by a spherical inclusion—I. Theory and numerical results, Geophys. J. Int.  115 230– 250. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1993b. Scattering of elastic waves by a spherical inclusion—II. Limitations of asymptotic solutions, Geophys. J. Int.  115 251– 263. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1996. Scattering of P and S waves by a spherically symmetric inclusion, Pure appl. Geophys.  147 675– 718. Google Scholar CrossRef Search ADS   Li L.-W., Kang X.-K., Leong M.-S., 2002. Spheroidal Wave Functions in Electromagnetic Theory  Wiley, 295 pp. Miles J.W., 1960. Scattering of elastic waves by small inhomogeneities, Geophysics  25 642– 648. Google Scholar CrossRef Search ADS   Mishchenko M.I., Hovenier J.W., Travis L.D., 2000. Light Scattering by Nonspherical Particles  Academic Press, 690 pp. Morse P.M., Feshbach H., 1953. Methods of Theoretical Physics  McGraw-Hill, 1978 pp. Sinha B.P., MacPhie R.H., Prasad T., 1973. Accurate computation of eigenvalues for prolate spheroidal wave functions, IEEE Trans. Antennas Propag.  21 406– 407. Google Scholar CrossRef Search ADS   Stratton J.A., Morse P.M., Chu L.J., Little J.D.S., Corbató F.J., 1956. Spheroidal Wave Functions  Wiley, 613 pp. Thompson W.J., 1997. Atlas for Computing Mathematical Functions  Wiley, 888 pp. van de Hulst H.C., 1957. Light Scattering by Small Particles  Wiley, 470 pp. Wu R., Aki K., 1985. Scattering characteristics of elastic waves by an elastic heterogeneity, Geophysics  50 582– 595. Google Scholar CrossRef Search ADS   Zhang S., Jin J., 1996. Computation of Special Functions  Wiley, 717 pp. APPENDIX A: DISPLACEMENTS AND TRACTIONS IN SPHEROIDAL COORDINATES Listed below are the components of the three displacement vectors L, M, N in spheroidal coordinates. The dependence on (η, ξ, ϕ, ω) as well as the superscript (j) in these vectors and the potential functions ψ, χ are omitted to improve clarity, and the standard notation of a subscript following a comma is used to denote a partial spatial derivative. Denoting the components of the displacement vectors in the form   $${\bf L} = L_{\eta } \hat{\bf a}_{\eta } + L_{\xi } \hat{\bf a}_{\xi } + L_{\phi } \hat{\bf a}_{\phi } ,$$ (A1)with similar equations for M and N, the individual components are listed below:   $$L_{\eta } = \frac{1}{c_{\alpha }} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}} \, \psi _{,\eta } ,$$ (A2a)  $$L_{\xi } = \frac{1}{c_{\alpha }} \sqrt{\frac{\xi ^2-1}{\xi ^2-\eta ^2}} \, \psi _{,\xi } ,$$ (A2b)  $$L_{\phi } = \frac{1}{c_{\alpha }}\frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}} \, \psi _{,\phi } ,$$ (A2c)  $$M_{\eta } = - \frac{\xi }{\sqrt{(1-\eta ^2)(\xi ^2-\eta ^2)}} \, \chi _{,\phi } ,$$ (A3a)  $$M_{\xi } = \frac{\eta }{\sqrt{(\xi ^2-1)(\xi ^2-\eta ^2)}} \, \chi _{,\phi } ,$$ (A3b)  $$M_{\phi } = \frac{\sqrt{(1-\eta ^2)(\xi ^2-1)}}{\xi ^2-\eta ^2} \, (\xi \chi _{,\eta }-\eta \chi _{,\xi }) ,$$ (A3c)  \begin{eqnarray} N_{\eta } & = & \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}} \Bigg [ \Big ( (\xi ^2-1) + 2 \xi ^2 \frac{1-\eta ^2}{\xi ^2-\eta ^2} \Big ) \chi _{,\eta } + 2 \eta \xi \frac{\xi ^2-1}{\xi ^2-\eta ^2} \chi _{,\xi } \nonumber \\ & &+\,\, \xi (\xi ^2-1)\chi _{,\eta \xi } - \eta \Big ( \Lambda _{mn}(c_\beta ) - m^2 \frac{1}{1-\eta ^2} - c_\beta ^2 \xi ^2 \Big ) \chi \Bigg ] , \end{eqnarray} (A4a)  \begin{eqnarray} N_{\xi } & = & \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \sqrt{\frac{\xi ^2-1}{\xi ^2-\eta ^2}} \Bigg [ -2 \eta \xi \frac{1-\eta ^2}{\xi ^2-\eta ^2} \chi _{,\eta } + \Big ( (1-\eta ^2) - 2\eta ^2\frac{\xi ^2-1)}{\xi ^2-\eta ^2} \Big ) \chi _{,\xi } \nonumber \\ & &+\,\, \eta (1-\eta ^2)\chi _{,\eta \xi } + \xi \Big ( \Lambda _{mn}(c_\beta ) + m^2 \frac{1}{\xi ^2-1} - c_\beta ^2 \eta ^2 \Big ) \chi \Bigg ] , \end{eqnarray} (A4b)  \begin{eqnarray} N_{\phi }= \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}} \left[ \eta (1-\eta ^2) \chi _{,\phi \eta } + \xi (\xi ^2-1) \chi _{,\phi \xi} +\,\, (\xi ^2-\eta ^2) \chi _{,\phi }) \right]. \end{eqnarray} (A4c) Next consider the tractions on a surface of constant ξ associated with the vector L written in component form as   $${\bf T} \big \lbrace {\bf L} \big \rbrace = T_{\eta } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\eta } + T_{\xi } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\xi } + T_{\phi } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\phi } ,$$ (A5)with similar expressions for T{M} and T{N}. The individual components of the tractions are listed below:   $$T_{\eta } \big \lbrace {\bf L} \big \rbrace = \frac{2\mu }{c_{\alpha }} \frac{2}{h} \frac{\sqrt{(1-\eta ^2)(\xi ^2-1)}}{\xi ^2-\eta ^2} \bigg [ -\frac{\xi }{\xi ^2-\eta ^2}\psi _{,\eta } +\frac{\eta }{\xi ^2-\eta ^2}\psi _{,\xi } + \psi _{,\eta \xi } \bigg ] ,$$ (A6a)  \begin{eqnarray} T_{\xi } \big \lbrace {\bf L} \big \rbrace & = & \frac{\kappa }{c_\alpha } \frac{2}{h} \Bigg \lbrace \left(\frac{2\mu }{3\kappa } -1\right) c_\alpha ^2 \psi + \frac{2\mu }{\kappa } \frac{1}{\xi ^2-\eta ^2} \bigg [ -\eta \frac{1-\eta ^2}{\xi ^2-\eta ^2}\psi _{,\eta } \nonumber \\ & & -\,\,\xi \frac{2(\xi ^2-1) + (1-\eta ^2)}{\xi ^2-\eta ^2}\psi _{,\xi } + \left( \Lambda _{mn}(c_\alpha ) + \frac{m^2}{\xi ^2-1} - c_\alpha ^2 \xi ^2 \right) \psi \bigg ] \Bigg \rbrace , \end{eqnarray} (A6b)  $$T_{\phi } \big \lbrace {\bf L} \big \rbrace = \frac{2\mu }{c_{\alpha }} \frac{2}{h} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-\eta ^2)}} \bigg [ \psi _{,\phi \xi } - \frac{\xi }{\xi ^2-1} \psi _{,\phi } \bigg ] ,$$ (A6c)  \begin{eqnarray} T_{\eta } \left\lbrace {\bf M} \right\rbrace & = & \mu \frac{2}{h} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}}\frac{1}{\xi ^2-\eta ^2} \left[ \eta (1-\eta ^2) \chi _{,\eta \phi } - \xi (\xi ^2-1) \chi _{,\xi \phi }+\,\, \left( (\xi ^2+1) - (1-\eta ^2) \right) \chi _{,\phi } \right], \end{eqnarray} (A7a)  $$T_{\xi } \big \lbrace {\bf M} \big \rbrace = 2 \mu \frac{2}{h} \frac{\eta }{\xi ^2-\eta ^2} \bigg [ \chi _{,\xi \phi } - \frac{\xi }{\xi ^2-1} \chi _{,\phi } \bigg ] ,$$ (A7b)  \begin{eqnarray} T_{\phi } \left\lbrace {\bf M} \right\rbrace & = & \mu \frac{2}{h} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}}\frac{1}{\xi ^2-\eta ^2} \bigg [ - \frac{(\xi ^2-1)[(\xi ^2+1)-(1-\eta ^2)]}{\xi ^2-\eta ^2} \chi _{,\eta } \nonumber \\ & & +\,\, 2 \eta \xi \frac{2(\xi ^2-1)+(1-\eta ^2)}{\xi ^2-\eta ^2} \chi _{,\xi } + \xi (\xi ^2-1) \chi _{,\eta \xi } \nonumber \\ & & -\,\, \eta \left( \Lambda _{mn}(c_\beta ) + m^2 [ \frac{2}{\xi ^2-1} + \frac{1}{1-\eta ^2} ] - c_\beta ^2 \xi ^2 \right) \chi \bigg ] , \end{eqnarray} (A7c)  \begin{eqnarray} T_{\eta } \left\lbrace {\bf N} \right\rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{\sqrt{(\xi ^2-1)(1-\eta ^2)}}{(\xi ^2-\eta ^2)^2} \bigg [ \xi \bigg ( \Lambda _{mn}(c_\beta ) -1 + \frac{m^2}{\xi ^2-1} - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \nonumber \\ & & +\,\, 9 \frac{1-\eta ^2}{\xi ^2-\eta ^2} - 12 \xi ^2 \frac{1-\eta ^2}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\eta } \nonumber \\ & & -\,\, \eta \bigg ( \Lambda _{mn}(c_\beta ) -1 - \frac{m^2}{1-\eta ^2} -\frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \nonumber \\ & & -\,\, 3 \frac{\xi ^2-1}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\xi } \nonumber \\ & & +\,\, \bigg ( 5\xi ^2 - 2 - (1-\eta ^2) - 6\xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \bigg ) \chi _{,\eta \xi } \nonumber \\ & & +\,\, \frac{6 \eta \xi }{\xi ^2-\eta ^2} \left( \Lambda _{mn}(c_\beta ) + \frac{m^2}{2} \left[ \frac{1}{\xi ^2-1} - \frac{1}{1-\eta ^2} \right] - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \right) \chi \bigg ] , \end{eqnarray} (A8a)  \begin{eqnarray} T_{\xi } \big \lbrace {\bf N} \big \rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{1}{(\xi ^2-\eta ^2)^2} \bigg [ \eta (1-\eta ^2) \bigg ( \Lambda _{mn}(c_\beta ) + \frac{m^2}{\xi ^2-1} - c_\beta ^2 \xi ^2 \nonumber \\ & & -\,\, \frac{7\xi ^2-3}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\eta } \nonumber \\ & & +\,\, \xi (\xi ^2-1) \bigg ( \Lambda _{mn}(c_\beta ) + 1 + \frac{m^2}{\xi ^2-1} - c_\beta ^2 + (1-\eta ^2)c_\beta ^2 \nonumber \\ & & -\,\, \frac{1-\eta ^2}{\xi ^2-1} - \frac{13\xi ^2-9}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\xi } \nonumber \\ & & -\,\, \xi \eta (1-\eta ^2) \left( 1 + 6 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right) \chi _{,\eta \xi } \nonumber \\ & & +\,\, \bigg ( \Lambda _{mn}(c_\beta ) \left[ 5\xi ^2 - 3 - 6 \xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right] + m^2 \left[ 3 - \frac{1}{\xi ^2-1} + \frac{1-\eta ^2}{\xi ^2-1} - 6 \frac{\xi ^2}{\xi ^2-\eta ^2} \right] \nonumber \\ & & -\,\, c_\beta ^2 \left[ 3\xi ^2 (2\xi ^2-1) - 1 - (2\xi ^2-1)(1-\eta ^2) - 6 \xi ^4 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right] \bigg ) \chi \bigg ] , \end{eqnarray} (A8b)  \begin{eqnarray} T_{\phi } \big \lbrace {\bf N} \big \rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{1}{\sqrt{1-\eta ^2}} \frac{1}{(\xi ^2-\eta ^2)^{3/2}} \bigg [ - \xi \eta (1-\eta ^2) \bigg ( \frac{1}{\xi ^2-1}+\frac{2}{\xi ^2-\eta ^2} \bigg ) \chi _{,\eta \phi } \nonumber \\ & & + \,\,\left( \xi ^2 - 2 + (1-\eta ^2) -2 \xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right) \chi _{,\xi \phi } + \eta (1-\eta ^2) \chi _{,\eta \xi \phi } \nonumber \\ & & + \,\,\xi \left( \Lambda _{mn}(c_\beta ) -1 + \frac{m^2}{\xi ^2-1} - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] - \frac{1-\eta ^2}{\xi ^2-1} \right) \chi _{,\phi } \bigg] . \end{eqnarray} (A8c) APPENDIX B: REMOVAL OF ϕ AND η DEPENDENCE FROM BOUNDARY EQUATIONS When written out in component form the boundary equations of eqs (22) and (23) consist of six scalar equations that contain the various components of displacement and traction listed in Appendix A that in turn contain the potential functions of eq. (8). These six equations are functions of the spheroidal coordinates (η, ξ, ϕ), and they are to be evaluated on the boundary where ξ = ξB. Furthermore, each term of these boundary equations involves infinite sums over the two indices m and n. The task of this appendix is to show how the dependence on η and ϕ in these equations can be removed. The removal of the ϕ dependence is straightforward and follows the same approach used in spherical problems. The only ϕ dependence in the boundary equations is the factor eimϕ that occurs in every term. Multiplying all of the boundary equations by $$e^{-im^{\prime } \phi } / 2 \pi$$ and then integrating over ϕ from 0 to 2π removes all of the ϕ dependence from the equations and leaves behind in each term the Kronecker delta function $$\delta _{m m^{\prime }}$$. Thus the infinite sums over m reduce to a single term where m = m΄. This means that the boundary equations separate with respect to the index m so that a separate set of equations exists for every value of m. The removal of the η dependence is more complicated and the general approach of Asano & Yamamoto (1975) is followed. The first step is to multiply the displacement boundary equations by the factors   $$u_\eta \ :\ \frac{(\xi ^2-\eta ^2)^{5/2}}{\xi (\xi ^2-1)^2} ,$$ (B1a)  $$u_\xi \ :\ \frac{(\xi ^2-\eta ^2)^{5/2}}{\xi ^2 (\xi ^2-1)^{3/2}} ,$$ (B1b)  $$u_\phi \ :\ \frac{(\xi ^2-\eta ^2)}{\xi (\xi ^2-1)^{1/2}} ,$$ (B1c)and to multiply the traction boundary equations by the factors   $$T_\eta \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^4}{(\xi ^2-1)^{7/2}} ,$$ (B1d)  $$T_\xi \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^4}{\xi (\xi ^2-1)^3} ,$$ (B1e)  $$T_\phi \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^{5/2}}{(\xi ^2-1)^2} .$$ (B1f) After this step all of the η dependence of the boundary equations is isolated into four types of factors   $$Q_{mn}^{(k)}(c,\eta ) = (1-\eta ^2)^{k/2} S_{mn}(c,\eta ) , \ \ -1 \le k \le 8 ,$$ (B2a)  $$\hat{Q}_{mn}^{(k)}(c,\eta ) = \eta Q_{mn}^{(k)}(c,\eta ) , \ \ -1 \le k \le 7 ,$$ (B2b)  $$Q_{mn}^{^{\prime }(k)}(c,\eta ) = (1-\eta ^2)^{k/2} \frac{d}{d\eta } S_{mn}(c,\eta ) , \ \ 1 \le k \le 7 ,$$ (B2c)  $$\hat{Q}_{mn}^{^{\prime }(k)}(c,\eta ) = \eta Q_{mn}^{^{\prime }(k)}(c,\eta ) , \ \ 1 \le k \le 7 ,$$ (B2d) where k is an integer in the specified range. A fundamental difficulty is now apparent. Although the wave functions Smn(c, η) are orthogonal with respect to the variable η (see eq. 13), the four types of factors all contain additional dependence upon η that prevents the use of this orthogonality to obtain a separate set of boundary equations for each value of n. While it is not possible to decouple the boundary equations with respect to the separation index n, it is possible to remove the η dependence from these equations. The basic approach is to substitute for the spheroidal angle functions Smn(c, η) a series expansion in terms of associated Legendre functions (eq. 12), and then use the completeness, orthogonality, and other properties of the associated Legendre functions to eliminate the η dependence. The approach varies slightly, depending upon the values of k and m. In four of the six boundary equations, the η and ϕ components of displacement and traction, the factors of eq. (B2) only appear with k as an odd integer. So first consider the case of k an odd integer and m not equal to zero. Let $$P_{m-1+l}^{m-1}(\eta )$$ be an appropriate selection function with normalization   $$N_{ml} = \int _{-1}^{+1} [ P_{m-1+l}^{m-1}(\eta ) ]^2 \, d \eta = \frac{2(2m-2+l)!}{(2m-1+2l)(l)!} ,$$ (B3)and note that l is an unspecified integer index at this point. Then the η dependence can be removed from the factors of eq. (B2) by performing the integrations   $$q_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} Q_{mn}^{(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4a)  $$\hat{q}_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} \hat{Q}_{mn}^{(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4b)  $$q_{mnl}^{^{\prime }(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} Q_{mn}^{^{\prime }(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} \hat{Q}_{mn}^{^{\prime }(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta .$$ (B4d) The remaining task is to evaluate the integrals in these expressions. Let the spheroidal angular function be expanded in associated Legendre functions (eq. 12) and then eq. (B4a) becomes   $$q_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \, d_r^{mn}(c) \int _{-1}^{+1} (1-\eta ^2)^{k/2} P_{m+r}^m(\eta ) P_{m-1+l}^{m-1}(\eta ) \, \text{d} \eta ,$$ (B5)with similar expressions for eqs (B4b)–(B4d). The integrals in these expressions can be evaluated by using properties of the associated Legendre functions (Erdélyi et al.1953). For the lowest required values of k the results are   $$q_{mnl}^{(-1)} (c) = - (2m-1+2l) \sum _{r=l}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{even} ,$$ (B6a)  $$\hat{q}_{mnl}^{(-1)} (c) = -l d_{l-1}^{mn}(c) - (2m-1+2l) \sum _{r=l+1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{odd} ,$$ (B6b)  $$q_{mnl}^{^{\prime }(1)}(c) = (m-1+l) l d_{l-1}^{mn}(c) - m(2m-1+2l) \sum _{r=l+1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{odd} ,$$ (B6c)  \begin{eqnarray} \hat{q}_{mnl}^{^{\prime }(1)}(c) & = & \frac{(m-2+l)(l-1)l}{2m-3+2l} d_{l-2}^{mn}(c) + \frac{1}{2} \left[ (l-1)l + \frac{(2m-1+l)(2m+l)}{2m+1+2l} \right] d_l^{mn}(c) \nonumber \\ & & -\,\, m(2m-1+2l) \sum _{r=l+2}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{even} . \end{eqnarray} (B6d) Regarded as a function of the index l, half of the integrals will be zero because they involve the product of even and odd functions of η, so in the above expressions the condition on whether (n − m + l) is even or odd is given to denote the non-zero values. This practice of listing only the non-zero values will apply in what follows. The approach outlined in the preceding paragraphs, which will be termed the direct approach, must be applied for the lowest values of k and can be applied for all values of k. However, for large values of k it becomes rather cumbersome and tedious. A more efficient and more systematic approach is to use the direct approach for the smallest values of k and then use recurrence relationships for the larger values of k. The recurrence process makes use of the following five functions:   $$R^{(-2)}(m,l) = - \frac{(2m-1+l)(2m+l)}{(2m-1+2l)(2m+1+2l)} ,$$ (B7a)  $$R^{(-1)}(m,l) = \frac{2m+l}{2m+1+2l} ,$$ (B7b)  $$R^{(0)}(m,l) = - \frac{1}{2m+3+2l} \left[ 1 + 2 \frac{(2m+l)l)}{2m-1+2l} \right] ,$$ (B7c)  $$R^{(1)}(m,l) = \frac{l+1}{2m+1+2l} ,$$ (B7d)  $$R^{(2)}(m,l) = - \frac{(l+1)(l+2)}{(2m+1+2l)(2m+3+2l)} .$$ (B7e) For the present case of k odd and m not equal to zero, the desired results can be written   $$q_{mnl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{even} ,$$ (B8a)  $$\hat{q}_{mnl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B8b)  $$q_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k-1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B8c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{even} .$$ (B8d) Assuming the expansion coefficients $$d_r^{mn}(c)$$ are already known, what remains is to determine the expansion coefficients $$p_j^{(k)}(m,r)$$. They are determined recursively by   $$p_j^{(k)}(m,l) = p_j^{(k-2)}(m,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{(k-2)}(m,l) ,$$ (B9a)  $$\hat{p}_j^{(k)}(m,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{(k)}(m,l) ,$$ (B9b)  $$p_j^{^{\prime }(k)}(m,l) = p_j^{^{\prime }(k-2)}(m,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{^{\prime }(k-2)}(m,l) ,$$ (B9c)  $$\hat{p}_j^{^{\prime }(k)}(m,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{^{\prime }(k)}(m,l) .$$ (B9d) In using these recursive equations values of Ri(m, l) are taken to be zero when the index i is not in the range −2 ≤ i ≤ 2. The first terms that are needed to start the process are   $$p_{-1}^{(1)}(m,l) = -\frac{(2m-1+l)(2m+l)}{2m+1+2l} ,$$ (B10a)  $$p_{1}^{(1)}(m,l) = \frac{(l+1)(l+2)}{2m+1+2l} ,$$ (B10b)  $$\hat{p}_{-2}^{^{\prime }(3)}(m,l) = - \frac{(2m-2+l)(2m-1+l)(2m+l)(m+1+l)}{(2m-l+2l)(2m+1+2l)} ,$$ (B10c)  $$\hat{p}_{0}^{^{\prime }(3)}(m,l) = \frac{(2m+l)(l+1)[2(m+l)(m+1+l)-3m]}{(2m-l+2l)(2m+3+2l)} ,$$ (B10d)  $$\hat{p}_{2}^{^{\prime }(3)}(m,l) = - \frac{(m+l)(l+1)(l+2)(l+3)}{(2m+l+2l)(2m+3+2l)} .$$ (B10e) With these initial values and the recurrence expressions of eq. (B9) the coefficients $$p_j^{k}(m,l)$$ can be determined for all values of (k, m, l) that are needed, and then the $$q_{mnl}^{(k)}(c)$$ are determined from eq. (B8). For example, for k = 1 eq. (B8a) gives   $$q_{mnl}^{(1)}(c) = p_{-1}^{(1)}(m,l) d_{l}^{mn}(c) + p_{1}^{(1)}(m,l-2) d_{l-2}^{mn}(c) ,$$ (B11a)with the coefficients $$p_{-1}^{(1)}(m,l)$$ and $$p_{1}^{(1)}(m,l)$$ given in eqs (B10a) and (B10b). Then eq. (B8b) gives   $$\hat{q}_{mnl}^{(1)}(c) = \hat{p}_{-2}^{(1)}(m,l+1) d_{l+1}^{mn}(c) + \hat{p}_{0}^{(1)}(m,l-1) d_{l-1}^{mn}(c) + \hat{p}_{2}^{(1)}(m,l-3) d_{l-3}^{mn}(c) ,$$ (B11b) with eq. (B9b) combined with eqs (B10a) and (B10b) to obtain   $$\hat{p}_{-2}^{(1)}(m,l) = R^{(-1)}(m-1,l) p_{-1}^{(1)}(m,l) ,$$ (B12a)  $$\hat{p}_{0}^{(1)}(m,l) = R^{(1)}(m-1,l) p_{-1}^{(1)}(m,l) + R^{(-1)}(m-1,l+2) p_{1}^{(1)}(m,l) ,$$ (B12b)  $$\hat{p}_{2}^{(1)}(m,l) = R^{(1)}(m-1,l+2) p_{1}^{(1)}(m,l) .$$ (B12c) Following this general process, all of the $$q_{mnl}^{\cdot (k)}(c)$$ functions listed in eq. (B4) can be evaluated for larger values of k. A second case that is only slightly different is that of m equal to zero and k an odd integer. The selection function is now taken to be $$P_{l+1}^{1}(\eta )$$ with the normalization factor   $$N_{0l} = \int _{-1}^{+1} [ P_{l+1}^{1}(\eta ) ]^2 \, \text{d} \eta = \frac{2(l+2)!}{(2l+3)(l)!} .$$ (B13)The direct approach produces   $$q_{0nl}^{(-1)} (c) = \frac{2l+3}{(l+1)(l+2)} \sum _{r=l+2}^{\infty \, ^{\prime \prime }} d_r^{0n} , \ (n+l) \, \text{even} ,$$ (B14a)  $$\hat{q}_{0nl}^{(-1)} (c) = \frac{1}{l+1} d_{l+1}^{0n} + \frac{2l+3}{(l+1)(l+2)} \sum _{r=l+3}^{\infty \, ^{\prime \prime }} d_r^{0n} , \ (n+l) \, \text{odd} ,$$ (B14b)  $$q_{0nl}^{^{\prime }(1)} (c) = - d_{l+1}^{0n}(c) , \ (n+l) \, \text{odd} .$$ (B14c)  $$\hat{q}_{0nl}^{^{\prime }(1)} (c) = - \frac{1}{2l+1} d_{l}^{0n}(c) - \frac{l+3}{2l+5} d_{l+2}^{0n}(c) , \ (n+l) \, \text{even} .$$ (B14d) With the recursion approach for larger values of k the results are   $$q_{0nl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15a)  $$\hat{q}_{0nl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15b)  $$q_{0nl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k-1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15c)  $$\hat{q}_{0nl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) .$$ (B15d) The expansion coefficients $$p_j^{(k)}(0,l)$$ now satisfy   $$p_j^{(k)}(0,l) = p_j^{(k-2)}(0,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{(k-2)}(0,l) ,$$ (B16a)  $$\hat{p}_j^{(k)}(0,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{(k)}(0,l) ,$$ (B16b)  $$p_j^{^{\prime }(k)}(0,l) = p_j^{^{\prime }(k-2)}(0,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{^{\prime }(k-2)}(0,l) ,$$ (B16c)  $$\hat{p}_j^{^{\prime }(k)}(0,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{^{\prime }(k)}(0,l) ,$$ (B16d) with starting values   $$p_{-1}^{(1)}(0,l) = \frac{1}{2l+1} ,$$ (B17a)  $$p_{1}^{(1)}(0,l) = - \frac{1}{2l+1} ,$$ (B17b)  $$p_{-2}^{^{\prime }(3)}(0,l) = \frac{l(l+l)}{(2l-1)(2l+1)} ,$$ (B17c)  $$p_{0}^{^{\prime }(3)}(0,l) = - 2 \frac{l(l+1)}{(2l-1)(2l+3)} ,$$ (B17d)  $$p_{2}^{^{\prime }(3)}(0,l) = \frac{l(l+1)}{(2l+1)(2l+3)} .$$ (B17e) In two of the six boundary equations, the ξ component of displacement and traction, the factors of eq. (B2) only appear with k as an even integer. The selection function is now $$P_{m+l}^{m}(\eta )$$ with the normalization factor   $$N_{ml} = \int _{-1}^{+1} [ P_{m+l}^{m}(\eta ) ]^2 \, \text{d} \eta = \frac{2(2m+l)!}{(2m+1+2l)(l)!} .$$ (B18)In this case the results for all values of k can be written   $$q_{mnl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{even} ,$$ (B19a)  $$\hat{q}_{mnl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B19b)  $$q_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k+1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B19c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{even} .$$ (B19d) The expansion coefficients $$p_j^{(k)}(m,l)$$ satisfy   $$p_j^{(k)}(m,l) = p_j^{(k-2)}(m,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{(k-2)}(m,l) ,$$ (B20a)  $$\hat{p}_j^{(k)}(m,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{(k)}(m,l) ,$$ (B20b)  $$p_j^{^{\prime }(k)}(m,l) = p_j^{^{\prime }(k-2)}(m,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{^{\prime }(k-2)}(m,l) ,$$ (B20c)  $$\hat{p}_j^{^{\prime }(k)}(m,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{^{\prime }(k)}(m,l) .$$ (B20d) The process begins with   $$p_0^{(0)}(m,l) = 1 ,$$ (B21a)  $$p_{-1}^{^{\prime }(2)}(m,l) = \frac{(2m+l)(m+1+l)}{(2m+1+2l)} ,$$ (B21b)  $$p_{1}^{^{\prime }(2)}(m,l) = - \frac{(m+l)(l+1)}{(2m+1+2l)} .$$ (B21c) APPENDIX C: COEFFICIENTS OF THE BOUNDARY EQUATIONS The boundary equations such as eqs (22) and (23) are functions of the vectors L, M and N (eq. 7) and each of these is a function of the spheroidal coordinates η, ξ and ϕ, as listed in Appendix A in component form. The process described in Appendix B is used to remove the η and ϕ dependence from the vector components and this leads to a transformation of the form   $$i^n L_{\eta mn}^{(j)} (\eta ,\xi ,\phi ,c_v) \rightarrow E_{\eta mnl}^{(Lj)}(v) ,$$ (C1a)  $$i^n T_{\eta } \big \lbrace {\bf L}_{mn}^{(j)}(\eta ,\xi ,\phi ,c_v) \big \rbrace \rightarrow F_{\eta mnl}^{(Lj)}(v) ,$$ (C1b) with similar transformations for the ξ and ϕ components and for the M and N vectors. With this transformation the boundary equations (22) and (23) become the boundary equations (24) and (25), with similar results for the ξ and ϕ components. The 18 coefficients of the transformed boundary equations are listed in this appendix. Note that in these coefficients the velocity v can take on any of the values α, β, α΄, or β΄, and the coordinate ξ is set equal to its boundary value ξ = ξB.   $$E_{\eta mnl}^{(Lj)}(v) = i^n \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } ,$$ (C2a)  \begin{eqnarray} E_{\eta mnl}^{(Nj)}(v) & = & i^n \Bigg [ \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{(3\xi ^2-1)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 \Big ) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} \Big ( 2 \hat{q}_{mnl}^{(1)}(c_v) + q_{mnl}^{^{\prime }(3)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C2b)  $$E_{\eta mnl}^{(Mj)}(v) = - i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) ,$$ (C2c)  $$E_{\xi mnl}^{(Lj)}(v) = i^n \frac{(\xi ^2-1)}{\xi ^2} \bigg \lbrace q_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} ,$$ (C3a)  \begin{eqnarray} E_{\xi mnl}^{(Nj)}(v) & = & i^n \Bigg [ \bigg \lbrace \Big [ \Lambda _{mn}(c_v) + \frac{m^2}{\xi ^2-1} \Big ] \Big [ q_{mnl}^{(0)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \Big ] \nonumber \\ & & -\,\, c_v^2 \Big [ q_{mnl}^{(0)}(c_v) - \frac{(\xi ^2-2)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) - \frac{1}{(\xi ^2-1)} q_{mnl}^{(4)}(c_c) \Big ] \nonumber \\ & & -\,\, \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(2)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v \xi } \nonumber \\ & & +\,\, \frac{1}{\xi ^2} \bigg \lbrace -2 q_{mnl}^{(0)}(c_v) + 3 q_{mnl}^{(2)}(c_v) + \hat{q}_{mnl}^{^{\prime }(2)}(c_v) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)} \Big [ q_{mnl}^{(4)}(c_v) + \hat{q}_{mnl}^{^{\prime }(4)}(c_v) \Big ] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C3b)  $$E_{\xi mnl}^{(Mj)}(v) = i^{n+1} \frac{m}{\xi ^2} \bigg \lbrace \hat{q}_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(4)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) ,$$ (C3c)  $$E_{\phi mnl}^{(Lj)}(v) = i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } ,$$ (C4a)  \begin{eqnarray} E_{\phi mnl}^{(Nj)}(v) & = & i^{n+1} m \left[ \left\lbrace q_{mn}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} \left[ q_{mn}^{(1)}(c_v) + \hat{q}_{mn}^{^{\prime }(1)}(c_v) \right] \right\rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } + q_{mnl}^{(-1)}(c_v) \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \right] , \end{eqnarray} (C4b)  $$E_{\phi mnl}^{(Mj)}(v) = i^n \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) R_{mn}^{(j)}(c_v,\xi ) - \hat{q}_{mnl}^{(1)}(c_v) \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{\xi } \bigg \rbrace ,$$ (C4c)  \begin{eqnarray} F_{\eta mnl}^{(Lj)}(v) & = & 2 \mu i^n \Bigg [ - \frac{\xi ^2}{(\xi ^2-1)}\bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)}[\hat{q}_{mnl}^{(1)}(c_v)+3 q_{mnl}^{^{\prime }(3)}(c_v)] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [2 \hat{q}_{mnl}^{(3)}(c_v) + 3 q_{mnl}^{^{\prime }(5)}(c_v)] + \frac{1}{(\xi ^2-1)^3} [\hat{q}_{mnl}^{(5)}(c_v) + q_{mnl}^{^{\prime }(7)}(c_v) ] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5a)  \begin{eqnarray} F_{\eta mnl}^{(Nj)}(v) & = & \frac{2 \mu i^n}{(\xi ^2-1)} \Bigg [ \xi ^2 \bigg \lbrace \Lambda _{mn}(c_v) \Big ( q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} [3 \hat{q}_{mnl}^{(1)}(c_v) + q_{mnl}^{^{\prime }(3)}(c_v)] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [6 \hat{q}_{mnl}^{(3)}(c_v) + q_{mnl}^{^{\prime }(5)}(c_v)] \Big ) \nonumber \\ & & -\,\, \frac{c_v^2}{2} \Big ( (\xi ^2+1) q_{mnl}^{^{\prime }(1)}(c_v) +6 \frac{\xi ^2+1}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{(\xi ^2+3)}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{(\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) + \frac{1}{(\xi ^2-1)^2} [12 \hat{q}_{mnl}^{(3)}(c_v)-6\hat{q}_{mnl}^{(5)}(c_v)-q_{mnl}^{^{\prime }(7)}(c_v)] \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( -3 \hat{q}_{mnl}^{(-1)}(c_v) + q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & + \,\,\frac{1}{(\xi ^2-1)^2} [ 3 \hat{q}_{mnl}^{(3)}(c_v) + q_{mnl}^{^{\prime }(5)}(c_v) ] \Big ) \nonumber \\ & & -\,\, q_{mnl}^{^{\prime }(1)}(c_v) - \frac{(5\xi ^2+7)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{8}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace - \Lambda _{mn}(c_v) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(5)}(c_v) \Big ) \nonumber \\ & & +\,\,\frac{c_v^2}{2} \Big ( (\xi ^2+1) \hat{q}_{mnl}^{(1)}(c_v) + \frac{(\xi ^2+3)}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \nonumber \\ & & -\,\, \frac{(\xi ^2-3)}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(5)}(c_v) - \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(7)}(c_v) \Big ) \nonumber \\ & & +\,\, m^2 \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & -\,\, (\xi ^2+2)q_{mnl}^{^{\prime }(1)}(c_v) + 3 q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} [ 4 (2\xi ^2+1) \hat{q}_{mnl}^{(1)}(c_v) - 5 \hat{q}_{mnl}^{(3)}(c_v) ] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [ 3 \xi ^2 q_{mnl}^{^{\prime }(5)}(c_v) + \hat{q}_{mnl}^{(5)}(c_v) - q_{mnl}^{^{\prime }(7)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5b)  \begin{eqnarray} F_{\eta mnl}^{(Mj)}(v) & = & \mu i^{n+1} m \Bigg [ \frac{1}{(\xi ^2-1)} \bigg \lbrace (\xi ^2+1) q_{mnl}^{(-1)}(c_v) + 2 \frac{(\xi ^2+2)}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \hat{q}_{mnl}^{^{\prime }(1)}(c_v) \nonumber \\ & & +\,\, \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(3)}(c_v) + \frac{3}{(\xi ^2-1)^2} [ 2q_{mnl}^{(3)}(c_v) + \hat{q}_{mnl}^{^{\prime }(5)}(c_v) ] \nonumber \\ & & -\,\, 2 \frac{(\xi ^2-2)}{(\xi ^2-1)^3} q_{mnl}^{(5)}(c_v) - \frac{1}{(\xi ^2-1)^3} [ q_{mnl}^{(7)}(c_v) - \hat{q}_{mnl}^{^{\prime }(7)}(c_v) ] \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & -\,\, c_v \xi \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{3}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v)+ \frac{3}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(5)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5c)  \begin{eqnarray} F_{\xi mnl}^{(Lj)}(v) & = & (\frac{2}{3} \mu - \kappa ) i^n c_v \xi \frac{(\xi ^2-1)}{\xi ^2} \Big ( q_{mnl}^{(0)}(c_v) + \frac{4}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{6}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{4}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) + \frac{1}{(\xi ^2-1)^4} q_{mnl}^{(8)}(c_v) \Big ) R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & +\,\, 2\mu i^n \Bigg [ \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 + \frac{m^2}{\xi ^2-1} \Big ) \Big ( q_{mnl}^{(0)}(c_v) + \frac{3}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{3}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} \Big ( \hat{q}_{mnl}^{^{\prime }(2)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(6)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v \xi } \nonumber \\ & & -\,\, \bigg \lbrace 2 q_{mnl}^{(0)}(c_v) + \frac{5}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{4}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C6a)  \begin{eqnarray} F_{\xi mnl}^{(Nj)}(v) & = & 2 \mu i^n \Bigg [ \frac{1}{(\xi ^2-1)} \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 + \frac{m^2}{\xi ^2-1} \Big ) \cdot \nonumber \\ & & \cdot \Big ( \hat{q}_{mnl}^{^{\prime }(2)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \Lambda _{mn}(c_v) \Big (- (\xi ^2+3) q_{mnl}^{(0)}(c_v) + 2 \frac{(2\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{(5\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \Big ) \nonumber \\ & & + c_v^2 \Big ( (3\xi ^2+1) q_{mnl}^{(0)}(c_v) - \frac{(4\xi ^4-3\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & -\,\, \frac{(2\xi ^4+3\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{(2\xi ^2-1)}{(\xi ^2-1)^2} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( -(3\xi ^2+4) q_{mnl}^{(0)}(c_v) + \frac{(\xi ^2-9)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\,\frac{5\xi ^2-6}{(\xi ^2-1)^2}q_{mnl}^{(4)}(c_v) +\frac{1}{(\xi ^2-1)^2} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & + \,\,\frac{(5\xi ^2+3)}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(2)}(c_v) - \frac{(7\xi ^2-3)}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & + \,\,\bigg \lbrace \Big ( \Lambda _{mn}(c_v) + \frac{m^2}{(\xi ^2-1)} \Big ) \Big ( q_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \Big ) \nonumber \\ & & -\,\, c_v^2 \Big ( q_{mnl}^{(0)}(c_v) - \frac{(\xi ^2-3)}{(\xi ^2-1}) q_{mnl}^{(2)}(c_v) - \frac{(2\xi ^2-3)}{(\xi ^2-1)^2}q_{mnl}^{(4)}(c_v)-\frac{1}{(\xi ^2-1)^2}q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)} \Big ( 8 q_{mnl}^{(0)}(c_v) - 4 \frac{(3\xi ^2-2)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & -\,\, 7 \hat{q}_{mnl}^{^{\prime }(2)}(c_v) -\frac{1}{(\xi ^2-1)} [ q_{mnl}^{(4)}(c_v) + 8 \hat{q}_{mnl}^{^{\prime }(4)}(c_v) ] \Big ) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)^2} [q_{mnl}^{(6)}(c_v) + \hat{q}_{mnl}^{^{\prime }(6)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Big ] , \end{eqnarray} (C6b)  \begin{eqnarray} F_{\xi mnl}^{(Mj)}(v) & = & 2 \mu i^{n+1} m \Bigg [ \hat{q}_{mnl}^{(0)}(c_v) + \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{(2)}(c_v) + \frac{3}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(4)}(c_v) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^3} \hat{q}_{mnl}^{(6)}(c_v) \Bigg ] \Bigg [ \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{\xi } - \frac{1}{(\xi ^2-1)} R_{mn}^{(j)}(c_v,\xi ) \Bigg ] , \end{eqnarray} (C6c)  \begin{eqnarray} F_{\phi mnl}^{(Lj)}(v) & = & 2\mu i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) \bigg \rbrace \nonumber \\ & & \bigg \lbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} - \frac{\xi ^2}{(\xi ^2-1)} \frac{ R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \bigg \rbrace , \end{eqnarray} (C7a)  \begin{eqnarray} F_{\phi mnl}^{(Nj)}(v) & = & \frac{2\mu i^{n+1} m}{(\xi ^2-1)} \Bigg [ \xi ^2 \bigg \lbrace \Big ( \Lambda _{mn}(c_v) + \frac{m^2}{(\xi ^2-1)} \Big ) \Big ( q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) \Big ) \nonumber \\ & & -\,\, \frac{c_v^2}{2} \Big ( (\xi ^2+1) q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} [ 2 q_{mnl}^{(1)}(c_v) - q_{mnl}^{(3)}(c_v) ] \Big ) \nonumber \\ & & -\,\, q_{mnl}^{(-1)}(c_v) - \frac{1}{(\xi ^2-1)} [ 2q_{mnl}^{(1)}(c_v)+3\hat{q}_{mnl}^{^{\prime }(1)}(c_v)] \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)^2} [q_{mnl}^{(3)}(c_v) + \hat{q}_{mnl}^{^{\prime }(3)}(c_v)] \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace -(\xi ^2+2) q_{mnl}^{(-1)}(c_v) + \frac{(2\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(1)}{c_v} \nonumber \\ & & +\,\, \hat{q}_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} [q_{mnl}^{(3)}{c_v} + \hat{q}_{mnl}^{^{\prime }(3)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] . \end{eqnarray} (C7b)  \begin{eqnarray} F_{\phi mnl}^{(Mj)}(v) & = & \mu i^n \Bigg [ - \frac{1}{(\xi ^2-1)} \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 \Big ) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, m^2 \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{2}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, (\xi ^2+1) q_{mnl}^{^{\prime }(1)}(c_v) - q_{mnl}^{^{\prime }(3)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & +\,\, c_v \xi \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} [ 4\hat{q}_{mnl}^{(1)}{c_v}+q_{mnl}^{^{\prime }(3)}(c_v) ] \nonumber \\ & & +\,\, \frac{2}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}{c_v} \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ]. \end{eqnarray} (C7c) Similar to the $$q_{mnl}^{(k)}$$ functions of Appendix B, many of the $$E_{.mnl}^{(.j)}$$ and $$F_{.mnl}^{(.j)}$$ functions are zero. The non-zero values are for (n − m + l) even   $$E_{\eta mnl}^{(Mj)}, \, E_{\xi mnl}^{(Lj)}, \, E_{\xi mnl}^{(Nj)}, \, E_{\phi mnl}^{(Lj)}, \, E_{\phi mnl}^{(Nj)}, \, F_{\eta mnl}^{(Mj)}, \, F_{\xi mnl}^{(Lj)}, \, F_{\xi mnl}^{(Nj)}, \, F_{\phi mnl}^{(Lj)}, \, F_{\phi mnl}^{(Nj)} ,$$ (C8a)and for (n − m + l) odd   $$E_{\eta mnl}^{(Lj)}, \, E_{\eta mnl}^{(Nj)}, \, E_{\xi mnl}^{(Mj)}, \, E_{\phi mnl}^{(Mj)}, \, F_{\eta mnl}^{(Lj)}, \, F_{\eta mnl}^{(Nj)}, \, F_{\xi mnl}^{(Mj)}, \, F_{\phi mnl}^{(Mj)} .$$ (C8b) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Scattering of elastic waves by a spheroidal inclusion

, Volume 212 (3) – Mar 1, 2018
30 pages

/lp/ou_press/scattering-of-elastic-waves-by-a-spheroidal-inclusion-0od9rAIMdp
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx482
Publisher site
See Article on Publisher Site

### Abstract

Summary An analytical solution is presented for scattering of elastic waves by prolate and oblate spheroidal inclusions. The problem is solved in the frequency domain where separation of variables leads to a solution involving spheroidal wave functions of the angular and radial kind. Unlike the spherical problem, the boundary equations remain coupled with respect to one of the separation indices. Expanding the angular spheroidal wave functions in terms of associated Legendre functions and using their orthogonality properties leads to a set of linear equations that can be solved to simultaneously obtain solutions for all coupled modes of both scattered and interior fields. To illustrate some of the properties of the spheroidal solution, total scattering cross-sections for P, SV and SH plane waves incident at an oblique angle on a prolate spheroid, an oblate spheroid and a sphere are compared. The waveforms of the scattered field exterior to the inclusion are calculated for these same incident waves. The waveforms scattered by a spheroid are strongly dependent upon the angle of incidence, are different for incident SV and SH waves and are asymmetrical about the centre of the spheroid with the asymmetry different for prolate and oblate spheroids. Wave propagation, Wave scattering and diffraction 1 INTRODUCTION The situation of an inclusion having material properties different from the surrounding medium is commonly encountered in studies of elastic wave propagation. When the shape of the inclusion has all dimensions approximately equal, the known solutions for scattering of elastic waves by a sphere are useful (Korneev & Johnson 1993a, 1996). However, in many cases, such as markedly elongated or flattened inclusions, the approximation by a sphere is questionable and more appropriate solutions are needed. The next step up from a sphere in terms of complexity is a spheroidal inclusion, which has the advantage of being a good approximation to many elongated or flattened shapes while still leading to differential equations that are separable. It appears that a complete solution for the scattering of elastic waves by a spheroidal inclusion has not yet appeared in the literature, so that is the purpose of the present study. Unlike the case of elastic waves, scattering of electromagnetic waves by a spheroidal inclusion has a well-developed literature (see Mishchenko et al. 2000, for a review). The elastic wave problem is more complicated than the electromagnetic wave problem in that it involves two velocities instead of one, requires two potential functions instead of one, allows three types of incident waves instead of two, and the relationship between displacement and stress fields lacks the symmetry of the relationship between electric and magnetic fields. However, essentially all of the methods and mathematics that have been developed for the electromagnetic problem can be applied to the elastic problem. Of the various methods that have been developed for electromagnetic problems, the present study follows most closely that of Asano & Yamamoto (1975). The approach followed in this study is to write the differential equations for elastic displacement in spheroidal coordinates and construct separate solutions to these equations for the incident field, the scattered field outside the inclusion and the interior field within the inclusion. The arbitrary constants in these fields are determined by requiring that boundary conditions on the surface of the inclusion are satisfied. This leads to a set of matrix equations than can be solved to yield an exact analytical solution to the scattering problem. There are a number of alternative approaches to the elastic scattering problem for a spheroidal inclusion that yield approximate solutions. One such approach is the transition matrix method that has been used by Hackman (1984) and Hackman & Lim (1994) to obtain the scattering solutions for inclusions of various shapes including a prolate spheroid. Another possible approach is to formulate the scattering problem in terms of an equivalent integral equation. While the integral formulation does not generally lead to an exact analytical solution, it is useful as the starting point for various approximations such as the Born series, low frequency, low contrast in material properties and far field. Examples of such approximate solutions can be found in Miles (1960), Gubernatis et al. (1977a,b), Datta (1977), Wu & Aki (1985) and Dassios & Kleinman (2000). 2 BASIC EQUATIONS IN SPHEROIDAL COORDINATES Consider a spheroidal inclusion with a polar axis rp having rotational symmetry and a perpendicular equatorial axis re. The degree to which the shape of the spheroid departs from that of a sphere is characterized by the interfocal distance h defined as   $$h = 2 \sqrt{r_p^2 - r_e^2} .$$ (1)It is also convenient to use the aspect ratio defined as   $$R_a = \frac{\text{max} \lbrace r_p , r_e\rbrace }{ \text{min} \lbrace r_p , r_e \rbrace } .$$ (2)When rp > re the spheroid is prolate and when rp < re it is oblate. Prolate and oblate spheroids are sketched in Figs 1 and 2, respectively, along with the geometry that is used later in numerical calculations. Figure 1. View largeDownload slide 2-Dview of a 3-D prolate spheroid having a polar axis rp and an equatorial axis re. An incident plane wave propagates from below in a positive z-direction that makes an angle ζ with the rp axis. Observations of the scattered waves are registered in the x, y and z directions on a plane that is displaced a distance zo from the centre of the spheroid. Figure 1. View largeDownload slide 2-Dview of a 3-D prolate spheroid having a polar axis rp and an equatorial axis re. An incident plane wave propagates from below in a positive z-direction that makes an angle ζ with the rp axis. Observations of the scattered waves are registered in the x, y and z directions on a plane that is displaced a distance zo from the centre of the spheroid. Figure 2. View largeDownload slide A sketch similar to Fig. 1 for an oblate spheroid. Figure 2. View largeDownload slide A sketch similar to Fig. 1 for an oblate spheroid. The spheroidal inclusion is imbedded in a homogeneous background having material properties density ρ, bulk modulus κ and shear modulus μ. The elastic compressional velocity is $$\alpha = [ (\kappa +\frac{4}{3}\mu )/\rho ]^{1/2}$$ and the shear velocity is β = [μ/ρ]1/2. Material properties within the inclusion will in general be different from the background and are denoted with a prime, so they are given by ρ΄, κ΄, μ΄, α΄, β΄. Considering first the prolate case, the spheroidal coordinates (η, ξ, ϕ) can be defined as an orthogonal right-handed system with their relationship to rectangular coordinates (x1, x2, x3) given by (Morse & Feshbach 1953; Flammer 1957)   $$x_1 = \frac{h}{2} \left[(\xi ^2-1)(1-\eta ^2)\right]^{1/2} \cos \phi ,$$ (3a)  $$x_2 = \frac{h}{2} \left[(\xi ^2-1)(1-\eta ^2)\right]^{1/2} \sin \phi ,$$ (3b)  $$x_3 = \frac{h}{2} \eta \xi .$$ (3c) with x3 parallel to the axis of rotational symmetry. These spheroidal coordinates have the ranges −1 ≤ η ≤ 1, 1 ≤ ξ ≤ ∞ and 0 ≤ ϕ ≤ 2π. The particle displacement is the vector u(x, t) where x is the position vector and t is time. In terms of angular frequency ω, this can be represented as   $${\bf u}({\bf x},t) = {\bf u}({\bf x},\omega ) \, e^{i \omega t} .$$ (4)The problem will be solved in the frequency domain, so the objective is to determine u(x, ω). In a source-free region conservation of linear momentum requires that   $$\alpha ^2 \nabla \left[ \nabla \cdot {\bf u}({\bf x},\omega ) \right] - \beta ^2 \nabla \times \left[ \nabla \times {\bf u}({\bf u},\omega ) \right] + \omega ^2 {\bf u}({\bf x},\omega ) = 0 .$$ (5)A general solution of this equation can be written in terms of three vectors (L, M, N)   $${\bf u}({\bf x},\omega ) = a^{(L)} {\bf L}({\bf x},\omega ) + a^{(M)} {\bf M}({\bf x},\omega ) + a^{(N)} {\bf N}({\bf x},\omega ) ,$$ (6)where a(L), a(M) and a(N) are constants. These vectors are related to two scalar potential functions ψ and χ that satisfy Helmholtz equations by   $${\bf L}({\bf x},\omega ) = \frac{1}{k_{\alpha }} \nabla \psi ({\bf x},\omega ), \ \ \left[ \nabla ^2 + k_{\alpha }^2 \right] \psi ({\bf x},\omega ) = 0 ,$$ (7a)  $${\bf M}({\bf x},\omega ) = \nabla \times [\chi ({\bf x},\omega ) {\bf r}], \ \ \left[ \nabla ^2 + k_{\beta }^2 \right] \chi ({\bf x},\omega ) = 0 ,$$ (7b)  $${\bf N}({\bf x},\omega ) = \frac{1}{k_{\beta }} \nabla \times {\bf M}({\bf x},\omega ) = \frac{1}{k_{\beta }} \nabla \times ( \nabla \times [ \chi ({\bf x},\omega ) {\bf r} ] ) ,$$ (7c) with wavenumbers kα = ω/α, kβ = ω/β and r the radial vector. Using this representation of u allows separation of variables to be used on eq. (5) in spheroidal coordinates (Morse & Feshbach 1953), so solutions can be written (Flammer 1957) in terms of angular wave functions eimϕ, angular wave functions Smn(c, η), and radial wave functions $$R_{mn}^{(j)}(c,\xi )$$, where m and n are separation indices. Solutions for the compressional and shear potential functions are   $$\psi _{mn}^{(j)}({\bf x},\omega ) = S_{mn}(c_\alpha ,\eta ) R_{mn}^{(j)}(c_\alpha ,\xi ) \, e^{im\phi } ,$$ (8a)and   $$\chi _{mn}^{(j)}({\bf x},\omega ) = S_{mn}(c_\beta ,\eta ) R_{mn}^{(j)}(c_\beta ,\xi ) \, e^{im\phi } ,$$ (8b) where   $$c_\alpha = \frac{h}{2} k_\alpha , \ \ c_\beta = \frac{h}{2} k_\beta .$$ (9)While the ϕ dependence is identical to the case of spherical coordinates, the other two wave functions satisfy the ordinary differential equations   $$\left[ \frac{\text{d}}{\text{d}\eta } (1-\eta ^2)\frac{\text{d}}{\text{d}\eta } - \frac{m^2}{1-\eta ^2} - c^2 \eta ^2 \right] S_{mn}(c,\eta ) = - \Lambda _{mn}(c) S_{mn}(c,\eta ) ,$$ (10a)  $$\left[ \frac{\text{d}}{\text{d}\xi } (\xi ^2-1)\frac{\text{d}}{\text{d}\xi } - \frac{m^2}{\xi ^2-1} + c^2 \xi ^2 \right] R_{mn}^{(j)}(c,\xi ) = \Lambda _{mn}(c) R_{mn}^{(j)}(c,\xi ) ,$$ (10b) where Λmn(c) is an eigenvalue. The two independent radial wave functions are denoted $$R_{mn}^{(1)}$$ and $$R_{mn}^{(2)}$$, where, similar to spherical Bessel and Neumann functions, the former is convergent at the origin and the latter is divergent. Similar to the spherical Hankel function of the second kind, a fourth radial wave function is defined as   $$R_{mn}^{(4)}(c,\xi ) = R_{mn}^{(1)}(c,\xi ) - i R_{mn}^{(2)}(c,\xi ) ,$$ (11)so that the combination of this with eq. (4) corresponds to an outward propagating wave at large ξ. The spheroidal angular function can be related to associated Legendre functions $$P_n^m(.)$$ by   $$S_{mn}(c,\eta ) = \sum _{r=0,1}^{\infty \, ^{\prime \prime }} d_r^{mn} (c) P_{m+r}^m(\eta ) ,$$ (12)where the $$d_r^{mn}(c)$$ are connection coefficients described by Flammer (1957). The double quotes over the summation sign indicates that when n − m is even the sum starts with r = 0 and increments by 2 so that only even values are included, whereas when n − m is odd it starts with r = 1 and increments by 2 so that only odd values are included. This convention will be used throughout this study. The normalization of these angular functions comes from the equation   $$\int _{-1}^{1} S_{mn}(c,\eta ) S_{ml}(c,\eta ) \, \text{d}\eta = N_{mn}^S(c) \delta _{nl} ,$$ (13)where   $$N_{mn}^S(c) = 2 \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{(2m+r)!}{(2m+2r+1)r!} \left[d_r^{mn}(c)\right]^2 .$$ (14) The stress tensor $$\tau$$ generated by a displacement u is (Morse & Feshbach 1953)   $$\tau \left\lbrace {\bf u} \right\rbrace = \left(\kappa - \frac{2}{3} \mu \right) (\nabla \cdot {\bf u}) {\bf I} + \mu \left[ \nabla {\bf u}+{\bf u}\nabla \right].$$ (15)What is actually needed for the boundary conditions is the traction on a surface of constant ξ having unit normal $$\hat{\bf a}_{\xi }$$ and this is the vector   $${\bf T} \left\lbrace {\bf u} \right\rbrace = \tau \left\lbrace {\bf u} \right\rbrace \cdot \hat{\bf a}_{\xi } .$$ (16)Specific expressions in spheroidal coordinates for the displacement vectors L, M, N and the associated tractions are given in Appendix A. Consider the boundary value problem for a general background material having material properties (ρ, α, β) surrounding an inclusion having material properties (ρ΄, α΄, β΄). The displacement of the scattered field outside the inclusion is   $${\bf u}^{\text{sca}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ a_{mn}^{(L)} {\bf L}_{mn}^{(4)}(c_{\alpha }) + a_{mn}^{(N)} {\bf N}_{mn}^{(4)}(c_{\beta }) + a_{mn}^{(M)} {\bf M}_{mn}^{(4)}(c_{\beta }) \right] ,$$ (17a)that of the field interior to the inclusion is   $${\bf u}^{\text{int}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ b_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha ^{\prime }}) + b_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta ^{\prime }}) + b_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta ^{\prime }}) \right] ,$$ (17b)and that of the incident field is   $${\bf u}^{\text{inc}} = \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \left[ s_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta }) \right] ,$$ (17c) The coefficients $$a_{mn}^{(.)}$$ and $$b_{mn}^{(.)}$$ are determined by the boundary conditions. The incident field is taken to be a plane wave parallel to the xy plane that propagates in a positive z-direction that makes an angle ζ with the rp axis of the spheroid (see Figs 1 and 2). First define three functions of the incident angle ζ as   $$f_{mn}^{(0)}(\zeta ) = \frac{2}{N_{mn}^S(c_\alpha )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c_\alpha ) P_{m+r}^m(\cos \zeta ) ,$$ (18a)  $$f_{mn}^{(1)}(\zeta ) = \frac{2m}{N_{mn}^S(c_\beta )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{d_r^{mn}(c_\beta )}{(r+m)(r+m+1)} \frac{P_{m+r}^m(\cos \zeta )}{\sin \zeta } ,$$ (18b)  $$f_{mn}^{(2)}(\zeta ) = \frac{2}{N_{mn}^S(c_\beta )} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \frac{d_r^{mn}(c_\beta )}{(r+m)(r+m+1)} \frac{\text{d}}{\text{d}\zeta } P_{m+r}^m(\cos \zeta ) .$$ (18c) Three different types of incident waves are possible with their properties listed in Table 1. For an incident P wave the non-zero expansion coefficients are   $$s_{mn}^{(L)}(\zeta ) = - i f_{mn}^{(0)}(\zeta ) .$$ (19a)For an incident SV wave the non-zero expansion coefficients are   $$s_{mn}^{(N)}(\zeta ) = - i f_{mn}^{(2)}(\zeta ) , \ \ s_{mn}^{(M)}(\zeta ) = - i f_{mn}^{(1)}(\zeta ) .$$ (19b)For an incident SH wave the non-zero expansion coefficients are   $$s_{mn}^{(N)}(\zeta ) = - f_{mn}^{(1)}(\zeta ) , \ \ s_{mn}^{(M)}(\zeta ) = - f_{mn}^{(2)}(\zeta ) .$$ (19c) Table 1. Types of incident waves.   Velocity  Displacement direction  Displacement plane  P  α  ∥ to direction of propagation  in plane containing ζ  SV  β  ⊥ to direction of propagation  in plane containing ζ  SH  β  ⊥ to direction of propagation  ⊥ to plane containing ζ    Velocity  Displacement direction  Displacement plane  P  α  ∥ to direction of propagation  in plane containing ζ  SV  β  ⊥ to direction of propagation  in plane containing ζ  SH  β  ⊥ to direction of propagation  ⊥ to plane containing ζ  View Large 3 BOUNDARY EQUATIONS The boundary conditions are that both displacement and traction should be continuous on the boundary of the inclusion. Letting the boundary of the inclusion be ξ = ξB the displacement boundary equation is   $${\bf u}^{\text{inc}}(\eta ,\xi _B,\phi ,\omega ) + {\bf u}^{\text{sca}}(\eta ,\xi _B,\phi ,\omega ) = {\bf u}^{\text{int}}(\eta ,\xi _B,\phi ,\omega ) ,$$ (20)and the traction boundary equation is   $${\bf T} \left\lbrace {\bf u}^{\text{inc}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace + {\bf T} \left\lbrace {\bf u}^{\text{sca}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace = {\bf T} \left\lbrace {\bf u}^{\text{int}}(\eta ,\xi _B,\phi ,\omega ) \right\rbrace .$$ (21) Written in component form, the first of the displacement equations for the η component is   \begin{eqnarray} \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty & & i^n \Big [ b_{mn}^{(L)} L_{\eta mn}^{(1)}(c_{\alpha ^{\prime }}) + b_{mn}^{(N)} N_{\eta mn}^{(1)}(c_{\beta ^{\prime }}) + b_{mn}^{(M)} M_{\eta mn}^{(1)}(c_{\beta ^{\prime }}) \nonumber \\ & & -\,\, a_{mn}^{(L)} L_{\eta mn}^{(4)}(c_{\alpha }) - a_{mn}^{(N)} N_{\eta mn}^{(4)}(c_{\beta }) - a_{mn}^{(M)} M_{\eta mn}^{(4)}(c_{\beta }) \Big ] \nonumber \\ &=& \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \Big [ s_{mn}^{(L)} L_{\eta mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} N_{\eta mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} M_{\eta mn}^{(1)}(c_{\beta }) \Big ] , \end{eqnarray} (22)with similar equations for the ξ and ϕ components. The boundary equation for the η component of traction is   \begin{eqnarray} \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty & & i^n \Big [ b_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(1)}(c_{\alpha ^{\prime }}) \big \rbrace + b_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(1)}(c_{\beta ^{\prime }}) \big \rbrace + b_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(1)}(c_{\beta ^{\prime }}) \big \rbrace \nonumber \\ & & - a_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(4)}(c_{\alpha }) \big \rbrace - a_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(4)}(c_{\beta }) \big \rbrace - a_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(4)}(c_{\beta }) \big \rbrace \Big ] \nonumber \\ &=& \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^\infty i^n \Big [ s_{mn}^{(L)} T_{\eta } \big \lbrace {\bf L}_{mn}^{(1)}(c_{\alpha }) \big \rbrace + s_{mn}^{(N)} T_{\eta } \big \lbrace {\bf N}_{mn}^{(1)}(c_{\beta }) \big \rbrace + s_{mn}^{(M)} T_{\eta } \big \lbrace {\bf M}_{mn}^{(1)}(c_{\beta }) \big \rbrace \Big ] , \end{eqnarray} (23)with similar equations for the ξ and ϕ components. At this point the six boundary equations that have to be satisfied are eq. (22) and similar displacement equations for the ξ and ϕ components together with eq. (23) and similar traction equations for the ξ and ϕ components. First note that the frequency ω, which enters through the variables cα and cβ (eq. 9) is fixed in these equations. Next note that (Appendix A and eq. 8) each of the components $$L_{.mn}^{(j)}$$, $$N_{.mn}^{(j)}$$, $$M_{.mn}^{(j)}$$, $$T_. \big \lbrace {\bf L}_{mn}^{(j)} \big \rbrace$$, $$T_. \big \lbrace {\bf N}_{mn}^{(j)} \big \rbrace$$ and $$T_. \big \lbrace {\bf M}_{mn}^{(j)} \big \rbrace$$ that appear in the boundary equations is a function of the coordinates (η, ξ, ϕ). The equations have to be satisfied for the single boundary value ξ = ξB and for all possible values of η and ϕ. Finally note that there is a double sum on each side of the equality sign in these equations. In the case of a spherical inclusion it is possible to use the orthogonality of the wave functions to remove from the boundary equations the dependence on η and ϕ and also the summations over the separation indexes m and n. In the case of a spheroidal inclusion all of this is also possible, except it is not possible to remove the summation over n and it is necessary to add another index l. This important difference from the spherical problem means that there is coupling between modes having different values of n and the boundary equations have to be solved simultaneously for all values of this index. The manipulations required to put the boundary equations in a form suitable for solution are given in Appendix B. After the application of the process described in Appendix B the boundary eq. (22) for the η component of displacement becomes   \begin{eqnarray} \sum _{n=|m|}^\infty & & b_{mn}^{(L)} E_{\eta mnl}^{(L1)}(\alpha ^{\prime }) + b_{mn}^{(N)} E_{\eta mnl}^{(N1)}(\beta ^{\prime }) + b_{mn}^{(M)} E_{\eta mnl}^{(M1)}(\beta ^{\prime }) \nonumber \\ & & - a_{mn}^{(L)} E_{\eta mnl}^{(L4)}(\alpha ) - a_{mn}^{(N)} E_{\eta mnl}^{(N4)}(\beta ) - a_{mn}^{(M)} E_{\eta mnl}^{(M4)}(\beta ) \nonumber \\ & = &\sum _{n=|m|}^\infty s_{mn}^{(L)} E_{\eta mnl}^{(L1)}(\alpha ) + s_{mn}^{(N)} E_{\eta mnl}^{(N1)}(\beta ) + s_{mn}^{(M)} E_{\eta mnl}^{(M1)}(\beta ) , \end{eqnarray} (24)and eq. (23) for the η component of traction becomes   \begin{eqnarray} \sum _{n=|m|}^\infty & & b_{mn}^{(L)} F_{\eta mnl}^{(L1)}(\alpha ^{\prime }) + b_{mn}^{(N)} F_{\eta mnl}^{(N1)}(\beta ^{\prime }) + b_{mn}^{(M)} F_{\eta mnl}^{(M1)}(\beta ^{\prime }) \nonumber \\ & & - a_{mn}^{(L)} F_{\eta mnl}^{(L4)}(\alpha ) - a_{mn}^{(N)} F_{\eta mnl}^{(N4)}(\beta ) - a_{mn}^{(M)} F_{\eta mnl}^{(M4)}(\beta ) \nonumber \\ & & = \sum _{n=|m|}^\infty s_{mn}^{(L)} F_{\eta mnl}^{(L1)}(\alpha ) + s_{mn}^{(N)} F_{\eta mnl}^{(N1)}(\beta ) + s_{mn}^{(M)} F_{\eta mnl}^{(M1)}(\beta ) , \end{eqnarray} (25)with similar equations for the ξ and ϕ components of both displacement and traction at the boundary. Expressions for $$E_{.mnl}^{(.j)}(v)$$ and $$F_{.mnl}^{(.j)}(v)$$ are given in Appendix C and, while they are functions of the boundary value ξ = ξB, they are not functions of η or ϕ. Also note the index m is fixed in these boundary equations, so a set of equations of this form have to be solved independently for each desired value of m. For given values of m and n the six boundary equations can be written in the matrix form   $${\bf C}_{mnl} \, {\bf w}_{mn} = {\bf S}_{mnl} \, {\bf s}_{mn} = {\bf D}_{mnl} ,$$ (26)where   $${\bf C}_{mnl} = \left[ \begin{array}{cccccc}E_{\eta mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\eta mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\eta mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\eta mnl}^{(L4)}(\alpha ) &\quad -E_{\eta mnl}^{(N4)}(\beta ) &\quad -E_{\eta mnl}^{(M4)}(\beta ) \\ E_{\xi mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\xi mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\xi mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\xi mnl}^{(L4)}(\alpha ) &\quad -E_{\xi mnl}^{(N4)}(\beta ) &\quad -E_{\xi mnl}^{(M4)}(\beta ) \\ E_{\phi mnl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\phi mnl}^{(N1)}(\beta ^{\prime }) &\quad E_{\phi mnl}^{(M1)}(\beta ^{\prime }) &\quad - E_{\phi mnl}^{(L4)}(\alpha ) &\quad -E_{\phi mnl}^{(N4)}(\beta ) &\quad -E_{\phi mnl}^{(M4)}(\beta ) \\ F_{\eta mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\eta mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\eta mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\eta mnl}^{(L4)}(\alpha ) &\quad -F_{\eta mnl}^{(N4)}(\beta ) &\quad -F_{\eta mnl}^{(M4)}(\beta ) \\ F_{\xi mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\xi mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\xi mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\xi mnl}^{(L4)}(\alpha ) &\quad -F_{\xi mnl}^{(N4)}(\beta ) &\quad -F_{\xi mnl}^{(M4)}(\beta ) \\ F_{\phi mnl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\phi mnl}^{(N1)}(\beta ^{\prime }) &\quad F_{\phi mnl}^{(M1)}(\beta ^{\prime }) &\quad - F_{\phi mnl}^{(L4)}(\alpha ) &\quad -F_{\phi mnl}^{(N4)}(\beta ) &\quad -F_{\phi mnl}^{(M4)}(\beta ) \\ \end{array} \right] ,$$ (27a)  $${\bf w}_{mn} = \left[ \begin{array}{c}b_{mn}^{(L)} \\ b_{mn}^{(N)} \\ b_{mn}^{(M)} \\ a_{mn}^{(L)} \\ a_{mn}^{(N)} \\ a_{mn}^{(M)} \end{array} \right] ,$$ (27b)  $${\bf S}_{mnl} = \left[ \begin{array}{ccc}E_{\eta mnl}^{(L1)}(\alpha ) &\quad E_{\eta mnl}^{(N1)}(\beta ) &\quad E_{\eta mnl}^{(M1)}(\beta ) \\ E_{\xi mnl}^{(L1)}(\alpha ) &\quad E_{\xi mnl}^{(N1)}(\beta ) &\quad E_{\xi mnl}^{(M1)}(\beta ) \\ E_{\phi mnl}^{(L1)}(\alpha ) &\quad E_{\phi mnl}^{(N1)}(\beta ) &\quad E_{\phi mnl}^{(M1)}(\beta ) \\ F_{\eta mnl}^{(L1)}(\alpha ) &\quad F_{\eta mnl}^{(N1)}(\beta ) &\quad F_{\eta mnl}^{(M1)}(\beta ) \\ F_{\xi mnl}^{(L1)}(\alpha ) &\quad F_{\xi mnl}^{(N1)}(\beta ) &\quad F_{\xi mnl}^{(M1)}(\beta ) \\ F_{\phi mnl}^{(L1)}(\alpha ) &\quad F_{\phi mnl}^{(N1)}(\beta ) &\quad F_{\phi mnl}^{(M1)}(\beta ) \\ \end{array} \right] , \ \ {\bf s}_{mn} = \left[ \begin{array}{c}s_{mn}^{(L)} \\ s_{mn}^{(N)} \\ s_{mn}^{(M)} \end{array} \right] ,$$ (27c)  $${\bf D}_{mnl} = \left[ \begin{array}{c}E_{\eta mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\eta mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\eta mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ E_{\xi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\xi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\xi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ E_{\phi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + E_{\phi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + E_{\phi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\eta mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\eta mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\eta mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\xi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\xi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\xi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ F_{\phi mnl}^{(L1)}(\alpha ) s_{mn}^{(L)} + F_{\phi mnl}^{(N1)}(\beta ) s_{mn}^{(N)} + F_{\phi mnl}^{(M1)}(\beta ) s_{mn}^{(M)} \\ \end{array} \right] .$$ (27d) Using the results from eq. (19) the source vector smn is selected as one of the following   $$\left[ {\bf s}_{mn}^{P} \ {\bf s}_{mn}^{SV} \ {\bf s}_{mn}^{SH} \right] = \left[ \begin{array}{ccc}-i f_{mn}^{(0)}(\zeta ) &\quad 0 &\quad 0 \\ 0 &\quad -i f_{mn}^{(2)}(\zeta ) &\quad - f_{mn}^{(1)}(\zeta ) \\ 0 &\quad -i f_{mn}^{(1)}(\zeta ) &\quad - f_{mn}^{(2)}(\zeta ) \\ \end{array} \right] .$$ (28) The boundary equations for given values of m and n lead to a system of six equations with six unknowns. However, because of coupling between modes it is necessary to simultaneously solve for all values of n. Assuming that the terms in the equations become negligibly small for n > m + J, the combined system for n = m to n = m + J can be written as   $$\left[ \begin{array}{ccccc}\bf C_{mml} &\quad {\bf C}_{m(m+1)l} &\quad {\bf C}_{m(m+2)l} &\quad \cdots &\quad {\bf C}_{m(m+J)l} \end{array} \right] \, \left[ \begin{array}{l}\bf w_{mm} \\ {\bf w}_{m(m+1)} \\ {\bf w}_{m(m+2)} \\ \vdots \\ {\bf w}_{m(m+J)} \end{array} \right] = \sum _{n=m}^{m+J} {\bf D}_{mnl} .$$ (29)This is a system of 6 equations with 6(J + 1) unknowns. This discrepancy is removed by recognizing that this system of equations holds for any value of l. Letting l vary from l = 0 to l = J leads to a square system of 6(J + 1) equations with 6(J + 1) unknowns:   $$\left[ \begin{array}{ccccc}\bf C_{mm0} &\quad {\bf C}_{m(m+1)0} &\quad {\bf C}_{m(m+2)0} &\quad \cdots &\quad {\bf C}_{m(m+J)0} \\ {\bf C}_{mm1} &\quad {\bf C}_{m(m+1)1} &\quad {\bf C}_{m(m+2)1} &\quad \cdots &\quad {\bf C}_{m(m+J)1} \\ {\bf C}_{mm2} &\quad {\bf C}_{m(m+1)2} &\quad {\bf C}_{m(m+2)2} &\quad \cdots &\quad {\bf C}_{m(m+J)2} \\ \vdots &\quad \vdots &\quad \vdots &\quad \cdots &\quad \vdots \\ {\bf C}_{mmJ} &\quad {\bf C}_{m(m+1)J} &\quad {\bf C}_{m(m+2)J} &\quad \cdots &\quad {\bf C}_{m(m+J)J} \\ \end{array} \right] \, \left[ \begin{array}{l}\bf w_{mm} \\ {\bf w}_{m(m+1)} \\ {\bf w}_{m(m+2)} \\ \vdots \\ {\bf w}_{m(m+J)} \end{array} \right] = \left[ \begin{array}{l}\sum _{n=m}^{m+J} {\bf D}_{mn0} \\ \sum _{n=m}^{m+J} {\bf D}_{mn1} \\ \sum _{n=m}^{m+J} {\bf D}_{mn2} \\ \vdots \\ \sum _{n=m}^{m+J} {\bf D}_{mnJ} \\ \end{array} \right] .$$ (30) In the special case of m = 0 several terms of the coefficient matrix Cmnl are identically zero and the 6 × 6 system of eq. (26) splits into two independent systems, a 4 × 4 system involving the L and N potentials and a 2 × 2 system involving the M potential. For the 4 × 4 system the individual matrices are   $${\bf C}_{0nl} = \left[ \begin{array}{cccc}E_{\eta 0nl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\eta 0nl}^{(N1)}(\beta ^{\prime }) &\quad - E_{\eta 0nl}^{(L4)}(\alpha ) &\quad -E_{\eta 0nl}^{(N4)}(\beta ) \\ E_{\xi 0nl}^{(L1)}(\alpha ^{\prime }) &\quad E_{\xi 0nl}^{(N1)}(\beta ^{\prime }) &\quad - E_{\xi 0nl}^{(L4)}(\alpha ) &\quad -E_{\xi 0nl}^{(N4)}(\beta ) \\ F_{\eta 0nl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\eta 0nl}^{(N1)}(\beta ^{\prime }) &\quad - F_{\eta 0nl}^{(L4)}(\alpha ) &\quad -F_{\eta 0nl}^{(N4)}(\beta ) \\ F_{\xi 0nl}^{(L1)}(\alpha ^{\prime }) &\quad F_{\xi 0nl}^{(N1)}(\beta ^{\prime }) &\quad - F_{\xi 0nl}^{(L4)}(\alpha ) &\quad -F_{\xi 0nl}^{(N4)}(\beta ) \\ \end{array} \right] ,$$ (31a)  $${\bf w}_{0n} = \left[ \begin{array}{c}b_{0n}^{(L)} \\ b_{0n}^{(N)} \\ a_{0n}^{(L)} \\ a_{0n}^{(N)} \end{array} \right] ,$$ (31b)  $${\bf D}_{0nl} = \left[ \begin{array}{c}E_{\eta 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + E_{\eta 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ E_{\xi 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + E_{\xi 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ F_{\eta 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + F_{\eta 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ F_{\xi 0nl}^{(L1)}(\alpha ) s_{0n}^{(L)} + F_{\xi 0nl}^{(N1)}(\beta ) s_{0n}^{(N)} \\ \end{array} \right] ,$$ (31c) and for the 2 × 2 system they are   $${\bf C}_{0nl} = \left[ \begin{array}{cc}E_{\phi 0nl}^{(M1)}(\beta ^{\prime }) &\quad -E_{\phi 0nl}^{(M4)}(\beta ) \\ F_{\phi 0nl}^{(M1)}(\beta ^{\prime }) &\quad -F_{\phi 0nl}^{(M4)}(\beta ) \\ \end{array} \right] ,$$ (32a)  $${\bf w}_{0n} = \left[ \begin{array}{c}b_{0n}^{(M)} \\ a_{0n}^{(M)} \end{array} \right] ,$$ (32b)  $${\bf D}_{0nl} = \left[ \begin{array}{c}E_{\phi 0nl}^{(M1)}(\beta ) s_{0n}^{(M)} \\ F_{\phi 0nl}^{(M1)}(\beta ) s_{0n}^{(M)} \\ \end{array} \right] .$$ (32c) Although not necessary, the dimension of this system of equations can be further reduced by taking advantage of the fact pointed out in Appendix C that many of the terms in the coefficient matrices of eqs (27a), (31a) and (32a) are zero, depending upon whether n − m + l is even or odd. For m ≠ 0, this means that the 6(J + 1) × 6(J + 1) general system can be split into two 3(J + 1) × 3(J + 1) systems. Let p = 0, 1, 2, 3, 4, ⋅⋅⋅ be an integer. Then the first system involves the coefficients   $$b_{m,m+2p+1}^{(L)} , b_{m,m+2p+1}^{(N)} , b_{m,m+2p}^{(M)} , a_{m,m+2p+1}^{(L)} , a_{m,m+2p+1}^{(N)} , a_{m,m+2p}^{(M)} ,$$ (33a)while the second involves   $$b_{m,m+2p}^{(L)} , b_{m,m+2p}^{(N)} , b_{m,m+2p+1}^{(M)} , a_{m,m+2p}^{(L)} , a_{m,m+2p}^{(N)} , a_{m,m+2p+1}^{(M)} .$$ (33b) For m = 0 the general system can be split into four systems. The first is a 2(J + 1) × 2(J + 1) system that involves   $$b_{0,2p}^{(L)} , b_{0,2p}^{(N)} , a_{0,2p}^{(L)} , a_{0,2p}^{(N)} .$$ (34a)The second is a 2(J + 1) × 2(J + 1) system that involves   $$b_{0,2p+1}^{(L)} , b_{0,2p+1}^{(N)} , a_{0,2p+1}^{(L)} , a_{0,2p+1}^{(N)} .$$ (34b)The third is a (J + 1) × (J + 1) system that involves   $$b_{0,2p}^{(M)} , a_{0,2p}^{(M)} .$$ (34c)The fourth is a (J + 1) × (J + 1) system that involves   $$b_{0,2p+1}^{(M)} , a_{0,2p+1}^{(M)} .$$ (34d) What is clear from these results is that for any of the expansion coefficients $$b_{mn}^{(.)}$$ or $$a_{mn}^{(.)}$$ there is no coupling between coefficients having different values of m. However, there is coupling between all coefficients having even values of n and between all coefficients having odd values of n, but no coupling between even and odd values. Furthermore, for m ≠ 0 there is coupling between coefficients having superscripts L, N and M. For m = 0 there is coupling between coefficients having superscripts L and N, but these are decoupled from coefficients having superscript M. The solution to the scattering problem is now complete. Given the geometry of the problem and the material properties of the inclusion and background, it is possible to calculate the spheroidal wave functions Smn(cv, η) and $$R_{mn}^{(j)}(c_v,\xi )$$. The requirement of continuity of displacement and traction on the boundary of the inclusion leads to equations in the form of (22) and (23), and applying the methods of Appendix B to these equations results in the boundary coefficients $$E_{.mnl}^{(.j)}$$ and $$F_{.mnl}^{(.j)}$$ (listed in Appendix C) that comprise the coefficient matrix of eq. (30). Selection of a source type from eq. (28) determines the right hand side of eq. (30). Solving the square linear system of equations in eq. (30) leads to the solution wmn, and then eq. (17) determines the part of the displacement due to scattering both outside and inside the inclusion. Note that this solution process is repeated for each value of the index m and each value of the frequency ω. The preceding development is for the prolate spheroidal coordinate system, but, as pointed out by Flammer (1957) and Asano & Yamamoto (1975), all of the results can be converted to the oblate spheroidal system through the substitutions   $$\xi \rightarrow i \xi , \ \ c \rightarrow - i c .$$ (35) 4 NUMERICAL CONSIDERATIONS In order to solve the boundary equations and obtain numerical values for the scattering coefficients it is necessary to calculate the eigenvalues Λmn(c), the expansion coefficients $$d_r^{mn}(c)$$, the angular functions Smn(c, η) and the radial functions $$R_{mn}^{(j)}(c,\xi )$$. In these calculations it is possible to take advantage of numerous resources available in the literature that were developed primarily for the scattering of electromagnetic waves. Morse & Feshbach (1953), Erdélyi et al. (1955), Stratton et al. (1956), Flammer (1957) and Abramowitz & Stegun (1964) all provide useful discussions of functions and algorithms that can be used in the calculations. Additional discussions of this type along with computer programs can be found in Zhang & Jin (1996), Thompson (1997), Li et al. (2002) and Falloon et al. (2003). More detailed discussions of specific aspects of the calculations can be found in Hodge (1970), Sinha et al. (1973), Do-Nhat & MacPhie (1996, 1997) and Do-Nhat (1999, 2001). While in theory the same algorithm can be used for both prolate and oblate calculations with the insertion of an i in the appropriate places (eq. 35), additional modifications may be necessary in order to obtain accurate results. The calculation of the spheroidal radial functions is particularly difficult, with different algorithms generally used for prolate and oblate functions and for large and small values of cξ. A practical problem involves the fact that the limits of the sums in the representations of incident and scattered fields in eq. (17) are infinite. In numerical calculations these limits have to be replaced by finite numbers, so the incident field (eq. 17c) becomes   $${\bf u}^{\text{inc}} = \sum _{m=-m_{\text{max}}}^{m_{\text{max}}} \sum _{n=|m|}^{n_{\text{max}}} i^n \left[ s_{mn}^{(L)} {\bf L}_{mn}^{(1)}(c_{\alpha }) + s_{mn}^{(N)} {\bf N}_{mn}^{(1)}(c_{\beta }) + s_{mn}^{(M)} {\bf M}_{mn}^{(1)}(c_{\beta }) \right] ,$$ (36)with similar modifications to the scattered and interior fields. This same problem is discussed by Korneev & Johnson (1993a) and a similar approach is followed here. The basic goal is to choose mmax and nmax large enough so that convergence of the series is achieved to some desired precision, but at the same time not so large that the terms in the series cannot be calculated to this same precision. Experience has shown that this goal is best obtained by considering the incident field, as this field converges more slowly than the scattered and interior fields. The calculations in the present study use   $$m_{\text{max}} = 10 + \frac{20}{\pi } \zeta + k_{\beta } r_s ,$$ (37a)  $$n_{\text{max}} = m + 10 + \frac{1}{2} \Big ( e + \frac{h}{r_s} \Big ) k_{\beta } r_s,$$ (37b) where ζ is given in units of radians and $$r_s = (r_p r_e^2)^{1/3}$$ is the spherical radius for a volume equal to that of the spheroid. Note that in the special case where ζ = 0, mmax = 1. The application of the boundary conditions leads to a set of linear equations in the form of eq. (30) that must be solved to obtain the solution coefficients wmm...wm(m + J). Formally this amounts to inverting a square matrix of dimension 6(J + 1) in the case where m > 0 and a smaller dimension when m = 0. In practice this is best achieved with an LU factorization followed by back substitution, which is performed using subroutines from LAPACK in the present study. 5 NUMERICAL RESULTS This section presents a few sample calculations that illustrate some of the properties of the elastic waves scattered by a spheroidal inclusion. The calculations are similar to some of those presented by Korneev & Johnson (1993a, 1996) for elastic waves scattered by a spherical inclusion and make use of the geometry shown in Figs 1 and 2. The model of material properties used for the calculations is listed in Table 2 and is a spheroidal inclusion that has velocities and density about 25 per cent less than the surrounding background material. This model is identical to model 2 of Korneev & Johnson (1996). Table 2. Model used for calculations.   Density  Compressional velocity  Shear velocity  Background  ρ = 2700 kgm−3  α = 6000 ms−1  β = 3500 ms−1  Inclusion  ρ΄ = 2300 kgm−3  α΄ = 4500 ms−1  β΄ = 2600 ms−1    Density  Compressional velocity  Shear velocity  Background  ρ = 2700 kgm−3  α = 6000 ms−1  β = 3500 ms−1  Inclusion  ρ΄ = 2300 kgm−3  α΄ = 4500 ms−1  β΄ = 2600 ms−1  View Large As shown in Figs 1 and 2, a plane wave is incident upon the spheroidal inclusion and solutions for the three components of the vector displacement usca(x, t) associated with the scattered field are calculated on a plane displaced a distance zo from the centre of the spheroid. The coordinate system (x, y, z) shown in these figures is rotated by an angle ζ with respect to the (x1, x2, x3) system of the spheroid given in eq. (3). Unlike the case of a spherical inclusion, the angle ζ between the direction of propagation of the incident wave and the polar axis of the spheroid is now an important parameter. It is necessary to distinguish three different types of incident waves and these were listed above in Table 1. In the geometry shown in Figs 1 and 2 the incident P wave has displacement in the z-direction, the incident SV wave in the x-direction, and the incident SH wave in the y-direction. In the case where ζ = 0 there is no distinction between the SV and SH incident waves, which is similar to the case of a spherical inclusion. Both prolate and oblate spheroidal inclusions are considered with both having an aspect ratio of 2 and both having a volume equal to that of a sphere having a radius of 1 m. Thus for the prolate inclusion   $$r_p = 1.588 \, m , \ \ r_e = 0.794 \, m ,$$ (38a)and for the oblate spheroid   $$r_p = 0.630 \, m , \ \ r_e = 1.260 \, m .$$ (38b)These dimensions give the boundary of the spheroid as ξB = 1.15 for the prolate case and ξB = 0.58 for the oblate case. The angle of incidence ζ and offset distance zo used for the calculations are   $$\zeta = 30^{\circ} , \ \ z_o = 2 \, m .$$ (38c) One method of characterizing the scattering of waves by an inclusion is to calculate and display scattering cross-sections as a function of frequency. The approach of Korneev & Johnson (1993a) is used to derive normalized total scattering cross-sections for the scattering of elastic waves by a spheroidal inclusion. Such cross-sections for incident P, SV and SH waves are   $$\sigma _{P} = - \frac{4\pi }{A_g(\zeta ) k_{\alpha }^2} \Im \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } a_{mn}^{(L)} W_{mn}^{(1)}(\zeta ) \right] ,$$ (39a)  $$\sigma _{SV} = \frac{4\pi }{A_g(\zeta ) k_{\beta }^2} \Im \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } \left( a_{mn}^{(N)} W_{mn}^{(3)}(\zeta ) + a_{mn}^{(M)} W_{mn}^{(2)}(\zeta ) \right) \right],$$ (39b)  $$\sigma _{SH} = - \frac{4\pi }{A_g(\zeta ) k_{\beta }^2} \Re \left[ \sum _{m=-\infty }^{\infty } \sum _{n=|m|}^{\infty } \left( a_{mn}^{(M)} W_{mn}^{(3)}(\zeta ) + a_{mn}^{(N)} W_{mn}^{(2)}(\zeta ) \right) \right] ,$$ (39c) where   $$W_{mn}^{(1)}(\zeta ) = P_n^m(\cos \zeta ) ,$$ (40a)  $$W_{mn}^{(2)}(\zeta ) = m \frac{P_n^m(\cos \zeta )}{\sin \zeta } ,$$ (40b)  $$W_{mn}^{(3)}(\zeta ) = \frac{d}{d \zeta } P_n^m(\cos \zeta ) ,$$ (40c) and the geometrical shadow of the inclusion is   $$A_g (\zeta ) = \pi r_e^2 \left[ 1 + \left( (r_p/r_e)^2 - 1 \right) \sin ^2(\zeta ) \right]^{1/2} .$$ (41)This geometrical shadow is the area of the plane incident wave that intersects the inclusion and is a function of the angle of incidence ζ. Normalizing the scattering cross-sections by this area allows them to be interpreted as the ratio of total energy scattered by the inclusion to total energy incident on the inclusion. The scattering cross-sections for incident P, SV and SH waves are plotted in Figs 3– 5, respectively, for the dimensions and angle of incidence listed in eq. (38). Each plot contains results for a prolate and oblate spheroid, both having an aspect ratio of 2, and a sphere, with all three inclusions having the same volume. The abscissa in each plot is a normalized frequency that is the wavenumber of the background compressional wave kα times the radius of a sphere having the same volume as the inclusion, which is rs = 1 m for Figs 3–5. This wavenumber is that of the incident wave in Fig. 3 but not in Figs 4 and 5 so that the three scattering cross-sections can be compared with the same horizontal scale. Figure 3. View largeDownload slide Normalized total scattering cross-sections for a plane P wave incident at an angle of 30° with respect to the polar axis of the inclusion. The solid line is for a sphere, the long dash line for a prolate spheroid and the short dash line for an oblate spheroid, where all three inclusions have the same volume of a sphere with a 1 m radius. Figure 3. View largeDownload slide Normalized total scattering cross-sections for a plane P wave incident at an angle of 30° with respect to the polar axis of the inclusion. The solid line is for a sphere, the long dash line for a prolate spheroid and the short dash line for an oblate spheroid, where all three inclusions have the same volume of a sphere with a 1 m radius. Figure 4. View largeDownload slide Similar to Fig. 3 for an incident plane SV wave. Figure 4. View largeDownload slide Similar to Fig. 3 for an incident plane SV wave. Figure 5. View largeDownload slide Similar to Fig. 3 for an incident plane SH wave. Figure 5. View largeDownload slide Similar to Fig. 3 for an incident plane SH wave. The scattering cross-sections in Figs 3–5 are all similar in that they exhibit characteristic Rayleigh type scattering at low frequencies where the amplitude is proportional to ω4 (Korneev & Johnson 1996). However, in this range the amplitudes for both the prolate and oblate spheroids are always less than for the sphere, regardless of the type of incident wave. This low frequency behaviour merges into large long oscillations that eventually tend to a constant value of 2.0 at high frequencies. Superimposed on these large long oscillations are numerous small short oscillations that are referred to as ripple. The basic features of the scattering cross-sections in Figs 3–5 are qualitatively similar to scattering cross-sections for electromagnetic waves incident upon a spherical inclusion, so interpretations developed for the electromagnetic problem can be borrowed for the elastic problem. Consider the scattering cross-section for incident P waves shown in Fig. 3. van de Hulst (1957) explains the large long oscillations as an interference pattern between waves transmitted through the inclusion and waves diffracted around it for electromagnetic waves. This interpretation applied to the elastic problem does reasonably well in predicting the peaks in the large long oscillations of the scattering cross-section for the sphere in Fig. 3 and can be extended to the cross-sections for the spheroids. The basic idea is that when waves passing through the inclusion are in-phase with those traveling around the inclusion there will be constructive interference and peaks in scattering cross-sections. Similarly, when these two types of waves are out-of-phase there will be destructive interference and troughs in the scattering cross-sections. The phase difference between waves passing through and around the inclusion depends upon frequency, velocity contrast of the inclusion, and dimension of the inclusion parallel to the direction of the incident wave. Noting that the projection of the inclusion onto the z-axis is 2.00 m for the sphere, 2.86 m for the prolate spheroid, and 1.67 m for the oblate spheroid (see Figs 1 and 2), the interpretation of van de Hulst (1957) predicts that the peaks in the large oscillations should be shifted to lower frequencies than the sphere for the prolate spheroid and to higher frequencies for the oblate spheroid, which agrees with the plots in Fig. 3. Also note that the first peak of the cross-section is larger for the prolate spheroid compared to the sphere by a factor of 1.08 and smaller for the oblate spheroid by a factor of 0.90. For cross-sections of electromagnetic waves incident upon a spherical inclusion Bohren & Huffman (1983) give an interpretation for the ripple in terms of the resonant modes of the inclusion. It seems likely that such an interpretation could apply to the case of elastic waves incident upon a spherical inclusion and also to the cases of prolate and oblate spheroids. It is important to note that the appearance of the ripple is influenced by the frequency sampling rate used in the calculations. The scattering cross-sections for incident SV and SH waves are shown in Figs 4 and 5, respectively. The general appearance of the cross-sections are similar to those of P waves but with all features shifted to lower frequencies. The velocities of S waves within the inclusion are about 0.7 times those of P waves. Again applying the interpretation of van de Hulst (1957), one would expect the peaks in the large oscillations to be shifted to lower frequencies by about a factor of about 0.7, and this is approximately the case. Similar to the case of an incident P wave, the first peak in the cross-section is shifted to lower frequencies for the prolate spheroid compared to the sphere and to higher frequencies for the oblate spheroid. In all cases the amplitude of the first peak of the scattering cross-section is larger for incident S waves than for incident P waves. For the prolate case the ratio of SV amplitude to P amplitude is 1.08 and the ratio for SH amplitude to P amplitude is 1.06. For the oblate case these ratios are both 1.01 The differences between the scattering cross-sections for an incident SV wave in Fig. 4 and an incident SH wave in Fig. 5 are small. This just reflects the fact that effects related to the polarization of the incident wave that might appear in differential cross-sections get averaged out when the integration is performed to obtain the total cross-sections. The scattering cross-sections in Figs 3–5 are for an angle of incidence ζ = 30° and an aspect ratio Ra = 2. A complete investigation of variations in these parameters as well as the velocity contrast of the inclusion is beyond the scope of the present study, but some of the effects caused by such variations can be anticipated. Asano (1979) considers a range in velocity contrasts, aspect ratios, and angles of incidence for the electromagnetic problem, and it appears that some of these results can be extended to the elastic problem, at least in a qualitative sense. As a check on this possibility consider the effect on the scattering cross-sections of Figs 3–5 when the angle of incidence is changed from ζ = 30° to ζ = 0°. The results for the sphere do not change but the peaks in the scattering cross sections shift to lower frequencies for the prolate spheroid and to higher frequencies for the oblate spheroid, which is consistent with the results of Asano (1979). These results are also consistent with the van de Hulst (1957) interpretation in terms of an interference pattern, as the projection of the inclusion onto the direction of the incident wave has now changed from 2.86 to 3.17 m for the prolate spheroid, suggesting a shift to lower frequencies, and from 1.67 to 1.26 m for the oblate spheroid, suggesting a shift to higher frequencies. This change in angle of incidence also causes the amplitude of the peaks to increase for the prolate spheroid and decrease for the oblate spheroid, which is also consistent with the findings of Asano (1979). The amount of ripple shows little change for this change in the angle of incidence. Given the general similarities between scattering cross-sections for electromagnetic and elastic waves, some of the other observations of Asano (1979) might also apply to scattering of elastic waves by spheroidal inclusions. A good test of any scattering solution is the ability to generate waveforms of the scattered fields, as this requires both modulus and phase to be accurately calculated over a broad frequency range. Figs 6– 15 display waveforms of scattered waves that are generated exterior to the inclusion for P, SV and SH types of incident waves. The waveforms are recorded along a profile that extends along the x-axis at a distance zo from the centre of the inclusion with the y coordinate set to zero (see Figs 1 and 2). The incident wave is a plane wave propagating in a direction that makes an angle ζ = 30° with the rp axis of the spheroidal inclusion, which has an aspect ratio of 2. Figs 6–10 are for a prolate spheroid and Figs 11– 15 are for an oblate spheroid. The calculations are performed with a Nyquist frequency of 12.5 KHz, and the incident wave is a delta function in time that is passed through a low-pass causal filter having a corner frequency of 6.25 KHz and a high frequency roll-off of 4. Figure 6. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 6. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 7. View largeDownload slide Similar to Fig. 6 for the scattered field with displacements in the z-direction. Figure 7. View largeDownload slide Similar to Fig. 6 for the scattered field with displacements in the z-direction. Figure 8. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 8. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 9. View largeDownload slide Similar to Fig. 8 for the scattered field with displacements in the z-direction. Figure 9. View largeDownload slide Similar to Fig. 8 for the scattered field with displacements in the z-direction. Figure 10. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 10. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 11. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 11. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane P wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 12. View largeDownload slide Similar to Fig. 11 for the scattered field with displacements in the z-direction. Figure 12. View largeDownload slide Similar to Fig. 11 for the scattered field with displacements in the z-direction. Figure 13. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 13. View largeDownload slide Waveforms of the scattered field with displacements in the x-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 14. View largeDownload slide Similar to Fig. 13 for the scattered field with displacements in the z-direction. Figure 14. View largeDownload slide Similar to Fig. 13 for the scattered field with displacements in the z-direction. Figure 15. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 15. View largeDownload slide Waveforms of the scattered field with displacements in the y-direction observed at points along the x-axis in Fig. 2 at an offset distance zo = 2 m for a plane SH wave incident at an angle ζ = 30° with respect to the polar axis of a oblate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in metres and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. In viewing the waveforms in Figs 6–9 and Figs 11–14 it is important to keep in mind that the scattering process causes conversions between P and SV waves and in general the converted waves are not propagating parallel to the incident wave. Thus in plots with either of these waves as the incident wave, it is not always obvious whether a particular waveform represents a P wave, an SV wave, or a combination of the two. Figs 6 and 7 show waveforms with displacements in the x and z directions, respectively, for a plane P wave incident on a prolate spheroid with velocities lower than the surrounding background material (Table 2). Only the scattered field outside the inclusion is shown, so all of these waveforms would have zero amplitude if the inclusion were not present. Because of the lower velocity within the inclusion, the waves passing through the inclusion become focused and have a fairly complicated pattern of triplications in the traveltime curves, caustics, and diffractions. The scattered waves extend well beyond the geometrical shadow of the inclusion, which has limits of ±1.05 m along the x coordinate in this case. For x positions inside the geometrical shadow the z component of displacement, which is that of the incident wave, is larger than the x component, but for x positions outside the geometrical shadow this situation can be reversed. Because of the tilt of the inclusion with ζ = 30° the waveforms are not symmetrical with respect to the x coordinate, with larger amplitudes and more extended durations for negative values of the x coordinate where the inclusion is in closer proximity to the observation point (see Fig. 1). For comparison, note that in the case of a spherical inclusion or a spheroidal inclusion with ζ = 0 the waveforms are symmetrical with respect to the x and y coordinates of the observation point. Figs 8 and 9 are similar to Figs 6 and 7 except that the incident wave is a plane SV wave with displacement in the x-direction. In many respects the waveforms are similar to the case of the incident P wave except the roles of the x and z displacement are reversed and the first arriving waves are delayed in time. The asymmetry with respect to the x-axis is now somewhat larger. Also note that in Fig. 9 prior to the time of the incident SV wave there are arrivals caused by P waves converted from the incident wave by the scattering process. The waveforms generated by an incident SH wave on a prolate inclusion are shown in Fig. 10 for a profile along the x-axis. The waveforms are for displacement in the y-direction, the same direction as the displacement of the incident wave. Because of the mirror symmetry of the inclusion about the xz plane, the displacements of the scattered waves in the x and z-directions are zero along this profile. The waveforms for this incident SH wave are quite similar in shape to those of the incident SV wave in Fig. 8, but the amplitudes are generally larger. This is because part of the energy of the incident SV wave is converted to displacement in the z-direction, as shown in Fig. 9. The waveforms in Figs 11–15 contain the same incident waves as those of Figs 6–10 except the inclusion is now oblate instead of prolate (compare Fig. 2 to Fig. 1). The limits of the geometrical shadow are now ±1.14 m along the x-axis. While there is a general similarity of the scattered waves for the prolate and oblate inclusions, there are also some specific differences. The amplitudes of the displacements for the oblate case are generally less than for the prolate case, which is consistent with the lower peak values of the scattering cross-sections in Figs 3–5. An exception is the y component of displacement for an incident SH wave in Fig. 15 where strong focussing causes amplitudes that are more than 3 times the amplitude of the incident wave. The waveforms are still asymmetrical with respect to the x coordinate, but now the larger amplitudes occur for positive values of the x coordinate where the inclusion is in closer proximity to the observation point (see Fig. 2). Fig. 14 is similar to Fig. 9 in that it shows converted P waves arriving before the incident SV wave. It is important to point out that the plots of the scattered waveforms in Figs 6–15 are only part of the realizable field that consists of the sum of the incident field and scattered field. In Fig. 16, the scattered field for the incident plane SV wave shown in Fig. 8 is combined with the incident field for that case to give the total field. Comparing Fig. 8 with Fig. 16 demonstrates the fact that one of the functions of the scattered field is to cancel out part of the incident field, as at small distances and times extending beyond the arrival time of the incident wave there is a shadow zone. Also note that outside the geometrical shadow of the inclusion the incident wave has a reduced amplitude, so the effect of the inclusion extends well beyond its geometrical shadow. Figure 16. View largeDownload slide Waveforms of the total field that consists of the scattered field of Fig. 8 plus the incident field for displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. Figure 16. View largeDownload slide Waveforms of the total field that consists of the scattered field of Fig. 8 plus the incident field for displacements in the x-direction observed at points along the x-axis in Fig. 1 at an offset distance zo = 2 m for a plane SV wave incident at an angle ζ = 30° with respect to the polar axis of a prolate inclusion having an aspect ratio of 2. The long dash vertical line denotes the arrival time of the incident wave and the two short dash horizontal lines denote the limits of the geometrical shadow. The number on the left is the x coordinate of the observation point in meters and the number on the right is the maximum amplitude of the waveform divided by the maximum amplitude of the incident wave. The calculations of this section used a maximum value for the wavenumber in the background material of kβ = 22 m−1. With ζ = 30° this required mmax = 35 (eq. 37a) and nmax = 107 (eq. 37b). The corresponding maximums in the spheroidal parameter c are cβ = 31 for the prolate case and cβ = 24 for the oblate case. When either the maximum frequency or aspect ratio of the spheroid is increased so that larger values of c are required, problems are encountered with the numerical accuracy of the radial wave functions and alternative methods of calculation may be required (Do-Nhat 1999, 2001; Kirby 2006, 2010). 6 DISCUSSION AND CONCLUSIONS One method of characterizing the solution for the scattering of elastic waves by a spheroidal inclusion is to compare it to the corresponding solution for a sphere. The spheroidal problem transforms to the spherical problem in the limit as   $$r_p \rightarrow r_e = r_s , \ h = 2 \sqrt{ r_p^2 - r_e^2} \rightarrow 0 , \ c_v = \frac{h}{2} k_v \rightarrow 0 ,$$ (42a)where v can be either the compressional velocity α or shear velocity β and kv = ω/v is a wavenumber. In this limit the spheroidal coordinates transform as   $$\eta \rightarrow \cos \theta , \ \xi \rightarrow \infty , \ \phi \rightarrow \phi ,$$ (42b)with the condition that   $$c_v \xi \rightarrow k_v r .$$ (42c)The rather complicated spheroidal eigenvalues transform to the simple form   $$\Lambda _{mn} (c_v) \rightarrow n (n+1) .$$ (42d)The angular spheroidal wave functions transform to associated Legendre functions   $$S_{mn} (c_v,\eta ) \rightarrow P_n^m (\cos \theta ) ,$$ (42e)and the radial spheroidal wave functions transform to spherical Bessel and Hankel functions   $$R_{mn}^{(1)} (c_v,\xi ) \rightarrow j_n (k_v r) , \ R_{mn}^{(2)} (c_v,\xi ) \rightarrow y_n (k_v r) , \ R_{mn}^{(4)} (c_v,\xi ) \rightarrow k_n^{(2)} (k_v r) .$$ (42f)In the general case ζ ≠ 0 the range of the separation index m is greatly reduced   $$0 \le m \le \infty \rightarrow m = 0, 1 ,$$ (42g)while the range of the separation index n is the same in the spheroidal and spherical problems. An obvious difference between the spheroidal and spherical problem is that the spheroidal problem actually consists of two different types of inclusions, prolate and oblate, with different computational methods leading to rather different scattered fields for the two types. This difference also represents a major asset of the solution to the spheroidal problem, the ability to model inclusions having a large variety in both shape and aspect ratio. Another very important difference between the spheroidal and spherical problems in a computational sense is the fact that the boundary equations remain coupled with respect to the separation index n in the case of a spheroid but are fully uncoupled in the case of a sphere. A few sample calculations are presented to illustrate some of the properties of the solution. Normalized total scattering cross-sections for incident P, SV and SH waves are compared for spherical, prolate, and oblate inclusions having velocities less than the surrounding background material. While the general features of the cross-sections are similar for these three cases, there are also some systematic differences in the amplitudes and frequency dependences. The interpretation of the scattering cross-sections can be guided by methods that have been developed for scattering cross sections of electromagnetic waves. Waveforms are calculated for the scattered field exterior to the inclusion for obliquely incident P, SV and SH waves. These calculations illustrate the importance of the angle ζ between the direction of propagation of the incident wave and the polar axis of the spheroid. In the case where ζ ≠ 0 some features of the spheroidal problem not found in the spherical problem include scattered fields that change as ζ is changed, a difference in the scattered fields for an incident SV wave and an incident SH wave, and scattered fields that do not have symmetry about the centre of the spheroid. Other more general scattering effects contained in the calculated waveforms are scattered fields that extend well outside the geometrical shadow of the inclusion, strong focussing for inclusions with velocities less than the background, and conversions between P and SV waves. While the solution presented in this study is complete in the sense that it is valid for any prolate or oblate spheroid, the application of this solution still has limits imposed by current capabilities to efficiently calculate spheroidal eigenvalues and wave functions for arbitrary combinations of material properties, aspect ratios, and frequency. Further progress on the numerical aspects of these calculations could help extend the utility of the spheroidal solution. Given the exact analytical solution presented in this study, it should be possible to test various approximate solutions such as those mentioned in the introduction, or possibly develop new approximate solutions (see Korneev & Johnson 1993b, for such approximations in the case of a spherical inclusion). For certain selected ranges of the variables and parameters, such approximate solutions may be completely acceptable and more computationally feasible than the exact solution. ACKNOWLEDGEMENTS This research is dedicated to the memory of Valeri Korneev and his contributions to problems involving the scattering of elastic waves by a spherical inclusion. REFERENCES Abramowitz M., Stegun I.A. (Eds), 1964. Handbook of Mathematical Functions  AMS 55, National Bureau of Standards, 1046 pp. Asano S., 1979. Light scattering properties of spheroidal particles, Appl. Opt.  18 712– 723. Google Scholar CrossRef Search ADS PubMed  Asano S., Yamamoto G., 1975. Light scattering by a spheroidal particle, Appl. Opt.  14 29– 49. Google Scholar CrossRef Search ADS PubMed  Bohren C.F., Huffman D.R., 1983. Absorption and Scattering of Light by Small Particles  Wiley, 530 pp. Dassios G., Kleinman R., 2000. Low Frequency Scattering  Clarendon Press, 297 pp. Datta S.K., 1977. Diffraction of plane elastic waves by ellipsoidal inclusions, J. acoust. Soc. Am.  61 1432– 1437. Google Scholar CrossRef Search ADS   Do-Nhat T., 1999. Asymptotic expansions of the Mathieu and prolate spheroidal eigenvalues for large parameter c, Can. J. Phys.  77 635– 652. Google Scholar CrossRef Search ADS   Do-Nhat T., 2001. Asymptotic expansions of the oblate spheroidal eigenvalues and wave functions for large parameter c, Can. J. Phys.  79 813– 831. Google Scholar CrossRef Search ADS   Do-Nhat T., MacPhie R.H., 1996. On the accurate computation of prolate spheroidal radial functions of the second kind, Q. Appl. Math.  55 677– 685. Google Scholar CrossRef Search ADS   Do-Nhat T., MacPhie R.H., 1997. Accurate values of prolate spheroidal radial functions of the second kind, Can. J. Phys.  75 671– 675. Google Scholar CrossRef Search ADS   Erdélyi A. (ed), Magnus W., Oberhettinger F., Tricomi F.G., 1953. Higher Transcendental Functions  Vol. 1, McGraw-Hill, 302 pp. Erdélyi A. (ed), Magnus W., Oberhettinger F., Tricomi F. G., 1955. Higher Transcendental Functions  Vol. 3, McGraw-Hill, 292 pp. Falloon P.E., Abbot P.C., Wang J.B., 2003. Theory and computation of spheroidal wave functions, J. Phys. A: Math. Gen.  36 5477– 5495. Google Scholar CrossRef Search ADS   Flammer C., 1957. Spheroidal Wave Functions  Stanford University Press, 220 pp. Gubernatis J.E., Domany E., Krumhansl J.A., 1977a. Formal aspects of the theory of the scattering of ultrasound by flaws in elastic materials, J. Appl. Phys.  48 2804– 2811. Google Scholar CrossRef Search ADS   Gubernatis J.E., Domany E., Krumhansl J.A., 1977b. The Born approximation in the theory of the scattering of elastic waves by flaws, J. Appl. Phys.  48 2812– 2819. Google Scholar CrossRef Search ADS   Hackman R.H., 1984. The transition matrix for acoustic and elastic wave scattering in prolate spheroidal coordinates, J. acoust. Soc. Am.  75 35– 45. Google Scholar CrossRef Search ADS   Hackman R.H., Lim R., 1994. Development and application of the spheroidal coordinate based T matrix solution to elastic wave scattering, Radio Sci.  29 1035– 1049. Google Scholar CrossRef Search ADS   Hodge D.P., 1970. Eigenvalues and eigenfunctions of the spheroidal wave equation, J. Math. Phys.  11 2308– 2312. Google Scholar CrossRef Search ADS   Kirby P., 2006. Calculation of spheroidal wave functions, Comput. Phys. Commun.  175 465– 472. Google Scholar CrossRef Search ADS   Kirby P., 2010. Calculation of radial prolate spheroidal wave functions of the second kind, Comput. Phys. Commun.  181 514– 519. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1993a. Scattering of elastic waves by a spherical inclusion—I. Theory and numerical results, Geophys. J. Int.  115 230– 250. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1993b. Scattering of elastic waves by a spherical inclusion—II. Limitations of asymptotic solutions, Geophys. J. Int.  115 251– 263. Google Scholar CrossRef Search ADS   Korneev V.A., Johnson L.R., 1996. Scattering of P and S waves by a spherically symmetric inclusion, Pure appl. Geophys.  147 675– 718. Google Scholar CrossRef Search ADS   Li L.-W., Kang X.-K., Leong M.-S., 2002. Spheroidal Wave Functions in Electromagnetic Theory  Wiley, 295 pp. Miles J.W., 1960. Scattering of elastic waves by small inhomogeneities, Geophysics  25 642– 648. Google Scholar CrossRef Search ADS   Mishchenko M.I., Hovenier J.W., Travis L.D., 2000. Light Scattering by Nonspherical Particles  Academic Press, 690 pp. Morse P.M., Feshbach H., 1953. Methods of Theoretical Physics  McGraw-Hill, 1978 pp. Sinha B.P., MacPhie R.H., Prasad T., 1973. Accurate computation of eigenvalues for prolate spheroidal wave functions, IEEE Trans. Antennas Propag.  21 406– 407. Google Scholar CrossRef Search ADS   Stratton J.A., Morse P.M., Chu L.J., Little J.D.S., Corbató F.J., 1956. Spheroidal Wave Functions  Wiley, 613 pp. Thompson W.J., 1997. Atlas for Computing Mathematical Functions  Wiley, 888 pp. van de Hulst H.C., 1957. Light Scattering by Small Particles  Wiley, 470 pp. Wu R., Aki K., 1985. Scattering characteristics of elastic waves by an elastic heterogeneity, Geophysics  50 582– 595. Google Scholar CrossRef Search ADS   Zhang S., Jin J., 1996. Computation of Special Functions  Wiley, 717 pp. APPENDIX A: DISPLACEMENTS AND TRACTIONS IN SPHEROIDAL COORDINATES Listed below are the components of the three displacement vectors L, M, N in spheroidal coordinates. The dependence on (η, ξ, ϕ, ω) as well as the superscript (j) in these vectors and the potential functions ψ, χ are omitted to improve clarity, and the standard notation of a subscript following a comma is used to denote a partial spatial derivative. Denoting the components of the displacement vectors in the form   $${\bf L} = L_{\eta } \hat{\bf a}_{\eta } + L_{\xi } \hat{\bf a}_{\xi } + L_{\phi } \hat{\bf a}_{\phi } ,$$ (A1)with similar equations for M and N, the individual components are listed below:   $$L_{\eta } = \frac{1}{c_{\alpha }} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}} \, \psi _{,\eta } ,$$ (A2a)  $$L_{\xi } = \frac{1}{c_{\alpha }} \sqrt{\frac{\xi ^2-1}{\xi ^2-\eta ^2}} \, \psi _{,\xi } ,$$ (A2b)  $$L_{\phi } = \frac{1}{c_{\alpha }}\frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}} \, \psi _{,\phi } ,$$ (A2c)  $$M_{\eta } = - \frac{\xi }{\sqrt{(1-\eta ^2)(\xi ^2-\eta ^2)}} \, \chi _{,\phi } ,$$ (A3a)  $$M_{\xi } = \frac{\eta }{\sqrt{(\xi ^2-1)(\xi ^2-\eta ^2)}} \, \chi _{,\phi } ,$$ (A3b)  $$M_{\phi } = \frac{\sqrt{(1-\eta ^2)(\xi ^2-1)}}{\xi ^2-\eta ^2} \, (\xi \chi _{,\eta }-\eta \chi _{,\xi }) ,$$ (A3c)  \begin{eqnarray} N_{\eta } & = & \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}} \Bigg [ \Big ( (\xi ^2-1) + 2 \xi ^2 \frac{1-\eta ^2}{\xi ^2-\eta ^2} \Big ) \chi _{,\eta } + 2 \eta \xi \frac{\xi ^2-1}{\xi ^2-\eta ^2} \chi _{,\xi } \nonumber \\ & &+\,\, \xi (\xi ^2-1)\chi _{,\eta \xi } - \eta \Big ( \Lambda _{mn}(c_\beta ) - m^2 \frac{1}{1-\eta ^2} - c_\beta ^2 \xi ^2 \Big ) \chi \Bigg ] , \end{eqnarray} (A4a)  \begin{eqnarray} N_{\xi } & = & \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \sqrt{\frac{\xi ^2-1}{\xi ^2-\eta ^2}} \Bigg [ -2 \eta \xi \frac{1-\eta ^2}{\xi ^2-\eta ^2} \chi _{,\eta } + \Big ( (1-\eta ^2) - 2\eta ^2\frac{\xi ^2-1)}{\xi ^2-\eta ^2} \Big ) \chi _{,\xi } \nonumber \\ & &+\,\, \eta (1-\eta ^2)\chi _{,\eta \xi } + \xi \Big ( \Lambda _{mn}(c_\beta ) + m^2 \frac{1}{\xi ^2-1} - c_\beta ^2 \eta ^2 \Big ) \chi \Bigg ] , \end{eqnarray} (A4b)  \begin{eqnarray} N_{\phi }= \frac{1}{c_\beta } \frac{1}{\xi ^2-\eta ^2} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}} \left[ \eta (1-\eta ^2) \chi _{,\phi \eta } + \xi (\xi ^2-1) \chi _{,\phi \xi} +\,\, (\xi ^2-\eta ^2) \chi _{,\phi }) \right]. \end{eqnarray} (A4c) Next consider the tractions on a surface of constant ξ associated with the vector L written in component form as   $${\bf T} \big \lbrace {\bf L} \big \rbrace = T_{\eta } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\eta } + T_{\xi } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\xi } + T_{\phi } \big \lbrace {\bf L} \big \rbrace \hat{\bf a}_{\phi } ,$$ (A5)with similar expressions for T{M} and T{N}. The individual components of the tractions are listed below:   $$T_{\eta } \big \lbrace {\bf L} \big \rbrace = \frac{2\mu }{c_{\alpha }} \frac{2}{h} \frac{\sqrt{(1-\eta ^2)(\xi ^2-1)}}{\xi ^2-\eta ^2} \bigg [ -\frac{\xi }{\xi ^2-\eta ^2}\psi _{,\eta } +\frac{\eta }{\xi ^2-\eta ^2}\psi _{,\xi } + \psi _{,\eta \xi } \bigg ] ,$$ (A6a)  \begin{eqnarray} T_{\xi } \big \lbrace {\bf L} \big \rbrace & = & \frac{\kappa }{c_\alpha } \frac{2}{h} \Bigg \lbrace \left(\frac{2\mu }{3\kappa } -1\right) c_\alpha ^2 \psi + \frac{2\mu }{\kappa } \frac{1}{\xi ^2-\eta ^2} \bigg [ -\eta \frac{1-\eta ^2}{\xi ^2-\eta ^2}\psi _{,\eta } \nonumber \\ & & -\,\,\xi \frac{2(\xi ^2-1) + (1-\eta ^2)}{\xi ^2-\eta ^2}\psi _{,\xi } + \left( \Lambda _{mn}(c_\alpha ) + \frac{m^2}{\xi ^2-1} - c_\alpha ^2 \xi ^2 \right) \psi \bigg ] \Bigg \rbrace , \end{eqnarray} (A6b)  $$T_{\phi } \big \lbrace {\bf L} \big \rbrace = \frac{2\mu }{c_{\alpha }} \frac{2}{h} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-\eta ^2)}} \bigg [ \psi _{,\phi \xi } - \frac{\xi }{\xi ^2-1} \psi _{,\phi } \bigg ] ,$$ (A6c)  \begin{eqnarray} T_{\eta } \left\lbrace {\bf M} \right\rbrace & = & \mu \frac{2}{h} \frac{1}{\sqrt{(1-\eta ^2)(\xi ^2-1)}}\frac{1}{\xi ^2-\eta ^2} \left[ \eta (1-\eta ^2) \chi _{,\eta \phi } - \xi (\xi ^2-1) \chi _{,\xi \phi }+\,\, \left( (\xi ^2+1) - (1-\eta ^2) \right) \chi _{,\phi } \right], \end{eqnarray} (A7a)  $$T_{\xi } \big \lbrace {\bf M} \big \rbrace = 2 \mu \frac{2}{h} \frac{\eta }{\xi ^2-\eta ^2} \bigg [ \chi _{,\xi \phi } - \frac{\xi }{\xi ^2-1} \chi _{,\phi } \bigg ] ,$$ (A7b)  \begin{eqnarray} T_{\phi } \left\lbrace {\bf M} \right\rbrace & = & \mu \frac{2}{h} \sqrt{\frac{1-\eta ^2}{\xi ^2-\eta ^2}}\frac{1}{\xi ^2-\eta ^2} \bigg [ - \frac{(\xi ^2-1)[(\xi ^2+1)-(1-\eta ^2)]}{\xi ^2-\eta ^2} \chi _{,\eta } \nonumber \\ & & +\,\, 2 \eta \xi \frac{2(\xi ^2-1)+(1-\eta ^2)}{\xi ^2-\eta ^2} \chi _{,\xi } + \xi (\xi ^2-1) \chi _{,\eta \xi } \nonumber \\ & & -\,\, \eta \left( \Lambda _{mn}(c_\beta ) + m^2 [ \frac{2}{\xi ^2-1} + \frac{1}{1-\eta ^2} ] - c_\beta ^2 \xi ^2 \right) \chi \bigg ] , \end{eqnarray} (A7c)  \begin{eqnarray} T_{\eta } \left\lbrace {\bf N} \right\rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{\sqrt{(\xi ^2-1)(1-\eta ^2)}}{(\xi ^2-\eta ^2)^2} \bigg [ \xi \bigg ( \Lambda _{mn}(c_\beta ) -1 + \frac{m^2}{\xi ^2-1} - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \nonumber \\ & & +\,\, 9 \frac{1-\eta ^2}{\xi ^2-\eta ^2} - 12 \xi ^2 \frac{1-\eta ^2}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\eta } \nonumber \\ & & -\,\, \eta \bigg ( \Lambda _{mn}(c_\beta ) -1 - \frac{m^2}{1-\eta ^2} -\frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \nonumber \\ & & -\,\, 3 \frac{\xi ^2-1}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\xi } \nonumber \\ & & +\,\, \bigg ( 5\xi ^2 - 2 - (1-\eta ^2) - 6\xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \bigg ) \chi _{,\eta \xi } \nonumber \\ & & +\,\, \frac{6 \eta \xi }{\xi ^2-\eta ^2} \left( \Lambda _{mn}(c_\beta ) + \frac{m^2}{2} \left[ \frac{1}{\xi ^2-1} - \frac{1}{1-\eta ^2} \right] - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] \right) \chi \bigg ] , \end{eqnarray} (A8a)  \begin{eqnarray} T_{\xi } \big \lbrace {\bf N} \big \rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{1}{(\xi ^2-\eta ^2)^2} \bigg [ \eta (1-\eta ^2) \bigg ( \Lambda _{mn}(c_\beta ) + \frac{m^2}{\xi ^2-1} - c_\beta ^2 \xi ^2 \nonumber \\ & & -\,\, \frac{7\xi ^2-3}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\eta } \nonumber \\ & & +\,\, \xi (\xi ^2-1) \bigg ( \Lambda _{mn}(c_\beta ) + 1 + \frac{m^2}{\xi ^2-1} - c_\beta ^2 + (1-\eta ^2)c_\beta ^2 \nonumber \\ & & -\,\, \frac{1-\eta ^2}{\xi ^2-1} - \frac{13\xi ^2-9}{\xi ^2-\eta ^2} + 12 \xi ^2 \frac{\xi ^2-1}{(\xi ^2-\eta ^2)^2} \bigg ) \chi _{,\xi } \nonumber \\ & & -\,\, \xi \eta (1-\eta ^2) \left( 1 + 6 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right) \chi _{,\eta \xi } \nonumber \\ & & +\,\, \bigg ( \Lambda _{mn}(c_\beta ) \left[ 5\xi ^2 - 3 - 6 \xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right] + m^2 \left[ 3 - \frac{1}{\xi ^2-1} + \frac{1-\eta ^2}{\xi ^2-1} - 6 \frac{\xi ^2}{\xi ^2-\eta ^2} \right] \nonumber \\ & & -\,\, c_\beta ^2 \left[ 3\xi ^2 (2\xi ^2-1) - 1 - (2\xi ^2-1)(1-\eta ^2) - 6 \xi ^4 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right] \bigg ) \chi \bigg ] , \end{eqnarray} (A8b)  \begin{eqnarray} T_{\phi } \big \lbrace {\bf N} \big \rbrace & = & \frac{2\mu }{c_\beta } \frac{2}{h} \frac{1}{\sqrt{1-\eta ^2}} \frac{1}{(\xi ^2-\eta ^2)^{3/2}} \bigg [ - \xi \eta (1-\eta ^2) \bigg ( \frac{1}{\xi ^2-1}+\frac{2}{\xi ^2-\eta ^2} \bigg ) \chi _{,\eta \phi } \nonumber \\ & & + \,\,\left( \xi ^2 - 2 + (1-\eta ^2) -2 \xi ^2 \frac{\xi ^2-1}{\xi ^2-\eta ^2} \right) \chi _{,\xi \phi } + \eta (1-\eta ^2) \chi _{,\eta \xi \phi } \nonumber \\ & & + \,\,\xi \left( \Lambda _{mn}(c_\beta ) -1 + \frac{m^2}{\xi ^2-1} - \frac{c_\beta ^2}{2} [(\xi ^2+1)-(1-\eta ^2)] - \frac{1-\eta ^2}{\xi ^2-1} \right) \chi _{,\phi } \bigg] . \end{eqnarray} (A8c) APPENDIX B: REMOVAL OF ϕ AND η DEPENDENCE FROM BOUNDARY EQUATIONS When written out in component form the boundary equations of eqs (22) and (23) consist of six scalar equations that contain the various components of displacement and traction listed in Appendix A that in turn contain the potential functions of eq. (8). These six equations are functions of the spheroidal coordinates (η, ξ, ϕ), and they are to be evaluated on the boundary where ξ = ξB. Furthermore, each term of these boundary equations involves infinite sums over the two indices m and n. The task of this appendix is to show how the dependence on η and ϕ in these equations can be removed. The removal of the ϕ dependence is straightforward and follows the same approach used in spherical problems. The only ϕ dependence in the boundary equations is the factor eimϕ that occurs in every term. Multiplying all of the boundary equations by $$e^{-im^{\prime } \phi } / 2 \pi$$ and then integrating over ϕ from 0 to 2π removes all of the ϕ dependence from the equations and leaves behind in each term the Kronecker delta function $$\delta _{m m^{\prime }}$$. Thus the infinite sums over m reduce to a single term where m = m΄. This means that the boundary equations separate with respect to the index m so that a separate set of equations exists for every value of m. The removal of the η dependence is more complicated and the general approach of Asano & Yamamoto (1975) is followed. The first step is to multiply the displacement boundary equations by the factors   $$u_\eta \ :\ \frac{(\xi ^2-\eta ^2)^{5/2}}{\xi (\xi ^2-1)^2} ,$$ (B1a)  $$u_\xi \ :\ \frac{(\xi ^2-\eta ^2)^{5/2}}{\xi ^2 (\xi ^2-1)^{3/2}} ,$$ (B1b)  $$u_\phi \ :\ \frac{(\xi ^2-\eta ^2)}{\xi (\xi ^2-1)^{1/2}} ,$$ (B1c)and to multiply the traction boundary equations by the factors   $$T_\eta \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^4}{(\xi ^2-1)^{7/2}} ,$$ (B1d)  $$T_\xi \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^4}{\xi (\xi ^2-1)^3} ,$$ (B1e)  $$T_\phi \ :\ \frac{h}{2} \frac{(\xi ^2-\eta ^2)^{5/2}}{(\xi ^2-1)^2} .$$ (B1f) After this step all of the η dependence of the boundary equations is isolated into four types of factors   $$Q_{mn}^{(k)}(c,\eta ) = (1-\eta ^2)^{k/2} S_{mn}(c,\eta ) , \ \ -1 \le k \le 8 ,$$ (B2a)  $$\hat{Q}_{mn}^{(k)}(c,\eta ) = \eta Q_{mn}^{(k)}(c,\eta ) , \ \ -1 \le k \le 7 ,$$ (B2b)  $$Q_{mn}^{^{\prime }(k)}(c,\eta ) = (1-\eta ^2)^{k/2} \frac{d}{d\eta } S_{mn}(c,\eta ) , \ \ 1 \le k \le 7 ,$$ (B2c)  $$\hat{Q}_{mn}^{^{\prime }(k)}(c,\eta ) = \eta Q_{mn}^{^{\prime }(k)}(c,\eta ) , \ \ 1 \le k \le 7 ,$$ (B2d) where k is an integer in the specified range. A fundamental difficulty is now apparent. Although the wave functions Smn(c, η) are orthogonal with respect to the variable η (see eq. 13), the four types of factors all contain additional dependence upon η that prevents the use of this orthogonality to obtain a separate set of boundary equations for each value of n. While it is not possible to decouple the boundary equations with respect to the separation index n, it is possible to remove the η dependence from these equations. The basic approach is to substitute for the spheroidal angle functions Smn(c, η) a series expansion in terms of associated Legendre functions (eq. 12), and then use the completeness, orthogonality, and other properties of the associated Legendre functions to eliminate the η dependence. The approach varies slightly, depending upon the values of k and m. In four of the six boundary equations, the η and ϕ components of displacement and traction, the factors of eq. (B2) only appear with k as an odd integer. So first consider the case of k an odd integer and m not equal to zero. Let $$P_{m-1+l}^{m-1}(\eta )$$ be an appropriate selection function with normalization   $$N_{ml} = \int _{-1}^{+1} [ P_{m-1+l}^{m-1}(\eta ) ]^2 \, d \eta = \frac{2(2m-2+l)!}{(2m-1+2l)(l)!} ,$$ (B3)and note that l is an unspecified integer index at this point. Then the η dependence can be removed from the factors of eq. (B2) by performing the integrations   $$q_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} Q_{mn}^{(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4a)  $$\hat{q}_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} \hat{Q}_{mn}^{(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4b)  $$q_{mnl}^{^{\prime }(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} Q_{mn}^{^{\prime }(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta ,$$ (B4c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \frac{1}{N_{ml}} \int _{-1}^{+1} \hat{Q}_{mn}^{^{\prime }(k)}(c,\eta ) P_{m-1+l}^{m-l}(\eta ) \, \text{d} \eta .$$ (B4d) The remaining task is to evaluate the integrals in these expressions. Let the spheroidal angular function be expanded in associated Legendre functions (eq. 12) and then eq. (B4a) becomes   $$q_{mnl}^{(k)}(c) = \frac{1}{N_{ml}} \sum _{r=0,1}^{\infty \, ^{\prime \prime }} \, d_r^{mn}(c) \int _{-1}^{+1} (1-\eta ^2)^{k/2} P_{m+r}^m(\eta ) P_{m-1+l}^{m-1}(\eta ) \, \text{d} \eta ,$$ (B5)with similar expressions for eqs (B4b)–(B4d). The integrals in these expressions can be evaluated by using properties of the associated Legendre functions (Erdélyi et al.1953). For the lowest required values of k the results are   $$q_{mnl}^{(-1)} (c) = - (2m-1+2l) \sum _{r=l}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{even} ,$$ (B6a)  $$\hat{q}_{mnl}^{(-1)} (c) = -l d_{l-1}^{mn}(c) - (2m-1+2l) \sum _{r=l+1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{odd} ,$$ (B6b)  $$q_{mnl}^{^{\prime }(1)}(c) = (m-1+l) l d_{l-1}^{mn}(c) - m(2m-1+2l) \sum _{r=l+1}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{odd} ,$$ (B6c)  \begin{eqnarray} \hat{q}_{mnl}^{^{\prime }(1)}(c) & = & \frac{(m-2+l)(l-1)l}{2m-3+2l} d_{l-2}^{mn}(c) + \frac{1}{2} \left[ (l-1)l + \frac{(2m-1+l)(2m+l)}{2m+1+2l} \right] d_l^{mn}(c) \nonumber \\ & & -\,\, m(2m-1+2l) \sum _{r=l+2}^{\infty \, ^{\prime \prime }} d_r^{mn}(c) , \ (n-m+l) \, \text{even} . \end{eqnarray} (B6d) Regarded as a function of the index l, half of the integrals will be zero because they involve the product of even and odd functions of η, so in the above expressions the condition on whether (n − m + l) is even or odd is given to denote the non-zero values. This practice of listing only the non-zero values will apply in what follows. The approach outlined in the preceding paragraphs, which will be termed the direct approach, must be applied for the lowest values of k and can be applied for all values of k. However, for large values of k it becomes rather cumbersome and tedious. A more efficient and more systematic approach is to use the direct approach for the smallest values of k and then use recurrence relationships for the larger values of k. The recurrence process makes use of the following five functions:   $$R^{(-2)}(m,l) = - \frac{(2m-1+l)(2m+l)}{(2m-1+2l)(2m+1+2l)} ,$$ (B7a)  $$R^{(-1)}(m,l) = \frac{2m+l}{2m+1+2l} ,$$ (B7b)  $$R^{(0)}(m,l) = - \frac{1}{2m+3+2l} \left[ 1 + 2 \frac{(2m+l)l)}{2m-1+2l} \right] ,$$ (B7c)  $$R^{(1)}(m,l) = \frac{l+1}{2m+1+2l} ,$$ (B7d)  $$R^{(2)}(m,l) = - \frac{(l+1)(l+2)}{(2m+1+2l)(2m+3+2l)} .$$ (B7e) For the present case of k odd and m not equal to zero, the desired results can be written   $$q_{mnl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{even} ,$$ (B8a)  $$\hat{q}_{mnl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B8b)  $$q_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k-1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B8c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l-1-j) d_{l-1-j}^{mn} (c) , \ (n-m+l) \, \text{even} .$$ (B8d) Assuming the expansion coefficients $$d_r^{mn}(c)$$ are already known, what remains is to determine the expansion coefficients $$p_j^{(k)}(m,r)$$. They are determined recursively by   $$p_j^{(k)}(m,l) = p_j^{(k-2)}(m,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{(k-2)}(m,l) ,$$ (B9a)  $$\hat{p}_j^{(k)}(m,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{(k)}(m,l) ,$$ (B9b)  $$p_j^{^{\prime }(k)}(m,l) = p_j^{^{\prime }(k-2)}(m,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{^{\prime }(k-2)}(m,l) ,$$ (B9c)  $$\hat{p}_j^{^{\prime }(k)}(m,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(m-1,l+1+i) p_i^{^{\prime }(k)}(m,l) .$$ (B9d) In using these recursive equations values of Ri(m, l) are taken to be zero when the index i is not in the range −2 ≤ i ≤ 2. The first terms that are needed to start the process are   $$p_{-1}^{(1)}(m,l) = -\frac{(2m-1+l)(2m+l)}{2m+1+2l} ,$$ (B10a)  $$p_{1}^{(1)}(m,l) = \frac{(l+1)(l+2)}{2m+1+2l} ,$$ (B10b)  $$\hat{p}_{-2}^{^{\prime }(3)}(m,l) = - \frac{(2m-2+l)(2m-1+l)(2m+l)(m+1+l)}{(2m-l+2l)(2m+1+2l)} ,$$ (B10c)  $$\hat{p}_{0}^{^{\prime }(3)}(m,l) = \frac{(2m+l)(l+1)[2(m+l)(m+1+l)-3m]}{(2m-l+2l)(2m+3+2l)} ,$$ (B10d)  $$\hat{p}_{2}^{^{\prime }(3)}(m,l) = - \frac{(m+l)(l+1)(l+2)(l+3)}{(2m+l+2l)(2m+3+2l)} .$$ (B10e) With these initial values and the recurrence expressions of eq. (B9) the coefficients $$p_j^{k}(m,l)$$ can be determined for all values of (k, m, l) that are needed, and then the $$q_{mnl}^{(k)}(c)$$ are determined from eq. (B8). For example, for k = 1 eq. (B8a) gives   $$q_{mnl}^{(1)}(c) = p_{-1}^{(1)}(m,l) d_{l}^{mn}(c) + p_{1}^{(1)}(m,l-2) d_{l-2}^{mn}(c) ,$$ (B11a)with the coefficients $$p_{-1}^{(1)}(m,l)$$ and $$p_{1}^{(1)}(m,l)$$ given in eqs (B10a) and (B10b). Then eq. (B8b) gives   $$\hat{q}_{mnl}^{(1)}(c) = \hat{p}_{-2}^{(1)}(m,l+1) d_{l+1}^{mn}(c) + \hat{p}_{0}^{(1)}(m,l-1) d_{l-1}^{mn}(c) + \hat{p}_{2}^{(1)}(m,l-3) d_{l-3}^{mn}(c) ,$$ (B11b) with eq. (B9b) combined with eqs (B10a) and (B10b) to obtain   $$\hat{p}_{-2}^{(1)}(m,l) = R^{(-1)}(m-1,l) p_{-1}^{(1)}(m,l) ,$$ (B12a)  $$\hat{p}_{0}^{(1)}(m,l) = R^{(1)}(m-1,l) p_{-1}^{(1)}(m,l) + R^{(-1)}(m-1,l+2) p_{1}^{(1)}(m,l) ,$$ (B12b)  $$\hat{p}_{2}^{(1)}(m,l) = R^{(1)}(m-1,l+2) p_{1}^{(1)}(m,l) .$$ (B12c) Following this general process, all of the $$q_{mnl}^{\cdot (k)}(c)$$ functions listed in eq. (B4) can be evaluated for larger values of k. A second case that is only slightly different is that of m equal to zero and k an odd integer. The selection function is now taken to be $$P_{l+1}^{1}(\eta )$$ with the normalization factor   $$N_{0l} = \int _{-1}^{+1} [ P_{l+1}^{1}(\eta ) ]^2 \, \text{d} \eta = \frac{2(l+2)!}{(2l+3)(l)!} .$$ (B13)The direct approach produces   $$q_{0nl}^{(-1)} (c) = \frac{2l+3}{(l+1)(l+2)} \sum _{r=l+2}^{\infty \, ^{\prime \prime }} d_r^{0n} , \ (n+l) \, \text{even} ,$$ (B14a)  $$\hat{q}_{0nl}^{(-1)} (c) = \frac{1}{l+1} d_{l+1}^{0n} + \frac{2l+3}{(l+1)(l+2)} \sum _{r=l+3}^{\infty \, ^{\prime \prime }} d_r^{0n} , \ (n+l) \, \text{odd} ,$$ (B14b)  $$q_{0nl}^{^{\prime }(1)} (c) = - d_{l+1}^{0n}(c) , \ (n+l) \, \text{odd} .$$ (B14c)  $$\hat{q}_{0nl}^{^{\prime }(1)} (c) = - \frac{1}{2l+1} d_{l}^{0n}(c) - \frac{l+3}{2l+5} d_{l+2}^{0n}(c) , \ (n+l) \, \text{even} .$$ (B14d) With the recursion approach for larger values of k the results are   $$q_{0nl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15a)  $$\hat{q}_{0nl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15b)  $$q_{0nl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k-1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) ,$$ (B15c)  $$\hat{q}_{0nl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l+1-j) d_{l+1-j}^{mn} (c) .$$ (B15d) The expansion coefficients $$p_j^{(k)}(0,l)$$ now satisfy   $$p_j^{(k)}(0,l) = p_j^{(k-2)}(0,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{(k-2)}(0,l) ,$$ (B16a)  $$\hat{p}_j^{(k)}(0,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{(k)}(0,l) ,$$ (B16b)  $$p_j^{^{\prime }(k)}(0,l) = p_j^{^{\prime }(k-2)}(0,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{^{\prime }(k-2)}(0,l) ,$$ (B16c)  $$\hat{p}_j^{^{\prime }(k)}(0,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(1,l-1+i) p_i^{^{\prime }(k)}(0,l) ,$$ (B16d) with starting values   $$p_{-1}^{(1)}(0,l) = \frac{1}{2l+1} ,$$ (B17a)  $$p_{1}^{(1)}(0,l) = - \frac{1}{2l+1} ,$$ (B17b)  $$p_{-2}^{^{\prime }(3)}(0,l) = \frac{l(l+l)}{(2l-1)(2l+1)} ,$$ (B17c)  $$p_{0}^{^{\prime }(3)}(0,l) = - 2 \frac{l(l+1)}{(2l-1)(2l+3)} ,$$ (B17d)  $$p_{2}^{^{\prime }(3)}(0,l) = \frac{l(l+1)}{(2l+1)(2l+3)} .$$ (B17e) In two of the six boundary equations, the ξ component of displacement and traction, the factors of eq. (B2) only appear with k as an even integer. The selection function is now $$P_{m+l}^{m}(\eta )$$ with the normalization factor   $$N_{ml} = \int _{-1}^{+1} [ P_{m+l}^{m}(\eta ) ]^2 \, \text{d} \eta = \frac{2(2m+l)!}{(2m+1+2l)(l)!} .$$ (B18)In this case the results for all values of k can be written   $$q_{mnl}^{(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} p_j^{(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{even} ,$$ (B19a)  $$\hat{q}_{mnl}^{(k)}(c) = \sum _{j=-(k+1)}^{(k+1) \, ^{\prime \prime }} \hat{p}_j^{(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B19b)  $$q_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-(k-1)}^{(k+1) \, ^{\prime \prime }} p_j^{^{\prime }(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{odd} ,$$ (B19c)  $$\hat{q}_{mnl}^{^{\prime }(k)}(c) = \sum _{j=-k}^{k \, ^{\prime \prime }} \hat{p}_j^{^{\prime }(k)}(m,l-j) d_{l-j}^{mn} (c) , \ (n-m+l) \, \text{even} .$$ (B19d) The expansion coefficients $$p_j^{(k)}(m,l)$$ satisfy   $$p_j^{(k)}(m,l) = p_j^{(k-2)}(m,l) + \sum _{i=-(k-2)}^{(k-2) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{(k-2)}(m,l) ,$$ (B20a)  $$\hat{p}_j^{(k)}(m,l) = \sum _{i=-k}^{k \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{(k)}(m,l) ,$$ (B20b)  $$p_j^{^{\prime }(k)}(m,l) = p_j^{^{\prime }(k-2)}(m,l) + \sum _{i=-(k-3)}^{(k-3) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{^{\prime }(k-2)}(m,l) ,$$ (B20c)  $$\hat{p}_j^{^{\prime }(k)}(m,l) = \sum _{i=-(k-1)}^{(k-1) \, ^{\prime \prime }} R^{(j-i)}(m,l+i) p_i^{^{\prime }(k)}(m,l) .$$ (B20d) The process begins with   $$p_0^{(0)}(m,l) = 1 ,$$ (B21a)  $$p_{-1}^{^{\prime }(2)}(m,l) = \frac{(2m+l)(m+1+l)}{(2m+1+2l)} ,$$ (B21b)  $$p_{1}^{^{\prime }(2)}(m,l) = - \frac{(m+l)(l+1)}{(2m+1+2l)} .$$ (B21c) APPENDIX C: COEFFICIENTS OF THE BOUNDARY EQUATIONS The boundary equations such as eqs (22) and (23) are functions of the vectors L, M and N (eq. 7) and each of these is a function of the spheroidal coordinates η, ξ and ϕ, as listed in Appendix A in component form. The process described in Appendix B is used to remove the η and ϕ dependence from the vector components and this leads to a transformation of the form   $$i^n L_{\eta mn}^{(j)} (\eta ,\xi ,\phi ,c_v) \rightarrow E_{\eta mnl}^{(Lj)}(v) ,$$ (C1a)  $$i^n T_{\eta } \big \lbrace {\bf L}_{mn}^{(j)}(\eta ,\xi ,\phi ,c_v) \big \rbrace \rightarrow F_{\eta mnl}^{(Lj)}(v) ,$$ (C1b) with similar transformations for the ξ and ϕ components and for the M and N vectors. With this transformation the boundary equations (22) and (23) become the boundary equations (24) and (25), with similar results for the ξ and ϕ components. The 18 coefficients of the transformed boundary equations are listed in this appendix. Note that in these coefficients the velocity v can take on any of the values α, β, α΄, or β΄, and the coordinate ξ is set equal to its boundary value ξ = ξB.   $$E_{\eta mnl}^{(Lj)}(v) = i^n \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } ,$$ (C2a)  \begin{eqnarray} E_{\eta mnl}^{(Nj)}(v) & = & i^n \Bigg [ \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{(3\xi ^2-1)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 \Big ) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} \Big ( 2 \hat{q}_{mnl}^{(1)}(c_v) + q_{mnl}^{^{\prime }(3)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C2b)  $$E_{\eta mnl}^{(Mj)}(v) = - i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) ,$$ (C2c)  $$E_{\xi mnl}^{(Lj)}(v) = i^n \frac{(\xi ^2-1)}{\xi ^2} \bigg \lbrace q_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} ,$$ (C3a)  \begin{eqnarray} E_{\xi mnl}^{(Nj)}(v) & = & i^n \Bigg [ \bigg \lbrace \Big [ \Lambda _{mn}(c_v) + \frac{m^2}{\xi ^2-1} \Big ] \Big [ q_{mnl}^{(0)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \Big ] \nonumber \\ & & -\,\, c_v^2 \Big [ q_{mnl}^{(0)}(c_v) - \frac{(\xi ^2-2)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) - \frac{1}{(\xi ^2-1)} q_{mnl}^{(4)}(c_c) \Big ] \nonumber \\ & & -\,\, \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(2)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v \xi } \nonumber \\ & & +\,\, \frac{1}{\xi ^2} \bigg \lbrace -2 q_{mnl}^{(0)}(c_v) + 3 q_{mnl}^{(2)}(c_v) + \hat{q}_{mnl}^{^{\prime }(2)}(c_v) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)} \Big [ q_{mnl}^{(4)}(c_v) + \hat{q}_{mnl}^{^{\prime }(4)}(c_v) \Big ] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C3b)  $$E_{\xi mnl}^{(Mj)}(v) = i^{n+1} \frac{m}{\xi ^2} \bigg \lbrace \hat{q}_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(4)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) ,$$ (C3c)  $$E_{\phi mnl}^{(Lj)}(v) = i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } ,$$ (C4a)  \begin{eqnarray} E_{\phi mnl}^{(Nj)}(v) & = & i^{n+1} m \left[ \left\lbrace q_{mn}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} \left[ q_{mn}^{(1)}(c_v) + \hat{q}_{mn}^{^{\prime }(1)}(c_v) \right] \right\rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } + q_{mnl}^{(-1)}(c_v) \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \right] , \end{eqnarray} (C4b)  $$E_{\phi mnl}^{(Mj)}(v) = i^n \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) R_{mn}^{(j)}(c_v,\xi ) - \hat{q}_{mnl}^{(1)}(c_v) \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{\xi } \bigg \rbrace ,$$ (C4c)  \begin{eqnarray} F_{\eta mnl}^{(Lj)}(v) & = & 2 \mu i^n \Bigg [ - \frac{\xi ^2}{(\xi ^2-1)}\bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)}[\hat{q}_{mnl}^{(1)}(c_v)+3 q_{mnl}^{^{\prime }(3)}(c_v)] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [2 \hat{q}_{mnl}^{(3)}(c_v) + 3 q_{mnl}^{^{\prime }(5)}(c_v)] + \frac{1}{(\xi ^2-1)^3} [\hat{q}_{mnl}^{(5)}(c_v) + q_{mnl}^{^{\prime }(7)}(c_v) ] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5a)  \begin{eqnarray} F_{\eta mnl}^{(Nj)}(v) & = & \frac{2 \mu i^n}{(\xi ^2-1)} \Bigg [ \xi ^2 \bigg \lbrace \Lambda _{mn}(c_v) \Big ( q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} [3 \hat{q}_{mnl}^{(1)}(c_v) + q_{mnl}^{^{\prime }(3)}(c_v)] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [6 \hat{q}_{mnl}^{(3)}(c_v) + q_{mnl}^{^{\prime }(5)}(c_v)] \Big ) \nonumber \\ & & -\,\, \frac{c_v^2}{2} \Big ( (\xi ^2+1) q_{mnl}^{^{\prime }(1)}(c_v) +6 \frac{\xi ^2+1}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{(\xi ^2+3)}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{(\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) + \frac{1}{(\xi ^2-1)^2} [12 \hat{q}_{mnl}^{(3)}(c_v)-6\hat{q}_{mnl}^{(5)}(c_v)-q_{mnl}^{^{\prime }(7)}(c_v)] \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( -3 \hat{q}_{mnl}^{(-1)}(c_v) + q_{mnl}^{^{\prime }(1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & + \,\,\frac{1}{(\xi ^2-1)^2} [ 3 \hat{q}_{mnl}^{(3)}(c_v) + q_{mnl}^{^{\prime }(5)}(c_v) ] \Big ) \nonumber \\ & & -\,\, q_{mnl}^{^{\prime }(1)}(c_v) - \frac{(5\xi ^2+7)}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(3)}(c_v) + \frac{8}{(\xi ^2-1)^2} q_{mnl}^{^{\prime }(5)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace - \Lambda _{mn}(c_v) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(5)}(c_v) \Big ) \nonumber \\ & & +\,\,\frac{c_v^2}{2} \Big ( (\xi ^2+1) \hat{q}_{mnl}^{(1)}(c_v) + \frac{(\xi ^2+3)}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \nonumber \\ & & -\,\, \frac{(\xi ^2-3)}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(5)}(c_v) - \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(7)}(c_v) \Big ) \nonumber \\ & & +\,\, m^2 \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & -\,\, (\xi ^2+2)q_{mnl}^{^{\prime }(1)}(c_v) + 3 q_{mnl}^{^{\prime }(3)}(c_v) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} [ 4 (2\xi ^2+1) \hat{q}_{mnl}^{(1)}(c_v) - 5 \hat{q}_{mnl}^{(3)}(c_v) ] \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^2} [ 3 \xi ^2 q_{mnl}^{^{\prime }(5)}(c_v) + \hat{q}_{mnl}^{(5)}(c_v) - q_{mnl}^{^{\prime }(7)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5b)  \begin{eqnarray} F_{\eta mnl}^{(Mj)}(v) & = & \mu i^{n+1} m \Bigg [ \frac{1}{(\xi ^2-1)} \bigg \lbrace (\xi ^2+1) q_{mnl}^{(-1)}(c_v) + 2 \frac{(\xi ^2+2)}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \hat{q}_{mnl}^{^{\prime }(1)}(c_v) \nonumber \\ & & +\,\, \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(3)}(c_v) + \frac{3}{(\xi ^2-1)^2} [ 2q_{mnl}^{(3)}(c_v) + \hat{q}_{mnl}^{^{\prime }(5)}(c_v) ] \nonumber \\ & & -\,\, 2 \frac{(\xi ^2-2)}{(\xi ^2-1)^3} q_{mnl}^{(5)}(c_v) - \frac{1}{(\xi ^2-1)^3} [ q_{mnl}^{(7)}(c_v) - \hat{q}_{mnl}^{^{\prime }(7)}(c_v) ] \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & -\,\, c_v \xi \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{3}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v)+ \frac{3}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(5)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C5c)  \begin{eqnarray} F_{\xi mnl}^{(Lj)}(v) & = & (\frac{2}{3} \mu - \kappa ) i^n c_v \xi \frac{(\xi ^2-1)}{\xi ^2} \Big ( q_{mnl}^{(0)}(c_v) + \frac{4}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{6}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{4}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) + \frac{1}{(\xi ^2-1)^4} q_{mnl}^{(8)}(c_v) \Big ) R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & +\,\, 2\mu i^n \Bigg [ \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 + \frac{m^2}{\xi ^2-1} \Big ) \Big ( q_{mnl}^{(0)}(c_v) + \frac{3}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{3}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)} \Big ( \hat{q}_{mnl}^{^{\prime }(2)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(6)}(c_v) \Big ) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v \xi } \nonumber \\ & & -\,\, \bigg \lbrace 2 q_{mnl}^{(0)}(c_v) + \frac{5}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\, \frac{4}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{1}{(\xi ^2-1)^3} q_{mnl}^{(6)}(c_v) \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] , \end{eqnarray} (C6a)  \begin{eqnarray} F_{\xi mnl}^{(Nj)}(v) & = & 2 \mu i^n \Bigg [ \frac{1}{(\xi ^2-1)} \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 + \frac{m^2}{\xi ^2-1} \Big ) \cdot \nonumber \\ & & \cdot \Big ( \hat{q}_{mnl}^{^{\prime }(2)}(c_v) + \frac{2}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) + \frac{1}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \Lambda _{mn}(c_v) \Big (- (\xi ^2+3) q_{mnl}^{(0)}(c_v) + 2 \frac{(2\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{(5\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \Big ) \nonumber \\ & & + c_v^2 \Big ( (3\xi ^2+1) q_{mnl}^{(0)}(c_v) - \frac{(4\xi ^4-3\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & -\,\, \frac{(2\xi ^4+3\xi ^2-3)}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) + \frac{(2\xi ^2-1)}{(\xi ^2-1)^2} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{m^2}{(\xi ^2-1)} \Big ( -(3\xi ^2+4) q_{mnl}^{(0)}(c_v) + \frac{(\xi ^2-9)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & +\,\,\frac{5\xi ^2-6}{(\xi ^2-1)^2}q_{mnl}^{(4)}(c_v) +\frac{1}{(\xi ^2-1)^2} q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & + \,\,\frac{(5\xi ^2+3)}{(\xi ^2-1)} \hat{q}_{mnl}^{^{\prime }(2)}(c_v) - \frac{(7\xi ^2-3)}{(\xi ^2-1)^2} \hat{q}_{mnl}^{^{\prime }(4)}(c_v) \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & + \,\,\bigg \lbrace \Big ( \Lambda _{mn}(c_v) + \frac{m^2}{(\xi ^2-1)} \Big ) \Big ( q_{mnl}^{(0)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(4)}(c_v) \Big ) \nonumber \\ & & -\,\, c_v^2 \Big ( q_{mnl}^{(0)}(c_v) - \frac{(\xi ^2-3)}{(\xi ^2-1}) q_{mnl}^{(2)}(c_v) - \frac{(2\xi ^2-3)}{(\xi ^2-1)^2}q_{mnl}^{(4)}(c_v)-\frac{1}{(\xi ^2-1)^2}q_{mnl}^{(6)}(c_v) \Big ) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)} \Big ( 8 q_{mnl}^{(0)}(c_v) - 4 \frac{(3\xi ^2-2)}{(\xi ^2-1)} q_{mnl}^{(2)}(c_v) \nonumber \\ & & -\,\, 7 \hat{q}_{mnl}^{^{\prime }(2)}(c_v) -\frac{1}{(\xi ^2-1)} [ q_{mnl}^{(4)}(c_v) + 8 \hat{q}_{mnl}^{^{\prime }(4)}(c_v) ] \Big ) \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)^2} [q_{mnl}^{(6)}(c_v) + \hat{q}_{mnl}^{^{\prime }(6)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Big ] , \end{eqnarray} (C6b)  \begin{eqnarray} F_{\xi mnl}^{(Mj)}(v) & = & 2 \mu i^{n+1} m \Bigg [ \hat{q}_{mnl}^{(0)}(c_v) + \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{(2)}(c_v) + \frac{3}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(4)}(c_v) \nonumber \\ & & +\,\, \frac{1}{(\xi ^2-1)^3} \hat{q}_{mnl}^{(6)}(c_v) \Bigg ] \Bigg [ \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{\xi } - \frac{1}{(\xi ^2-1)} R_{mn}^{(j)}(c_v,\xi ) \Bigg ] , \end{eqnarray} (C6c)  \begin{eqnarray} F_{\phi mnl}^{(Lj)}(v) & = & 2\mu i^{n+1} m \bigg \lbrace q_{mnl}^{(-1)}(c_v) + \frac{2}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)^2} q_{mnl}^{(3)}(c_v) \bigg \rbrace \nonumber \\ & & \bigg \lbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} - \frac{\xi ^2}{(\xi ^2-1)} \frac{ R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \bigg \rbrace , \end{eqnarray} (C7a)  \begin{eqnarray} F_{\phi mnl}^{(Nj)}(v) & = & \frac{2\mu i^{n+1} m}{(\xi ^2-1)} \Bigg [ \xi ^2 \bigg \lbrace \Big ( \Lambda _{mn}(c_v) + \frac{m^2}{(\xi ^2-1)} \Big ) \Big ( q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} q_{mnl}^{(1)}(c_v) \Big ) \nonumber \\ & & -\,\, \frac{c_v^2}{2} \Big ( (\xi ^2+1) q_{mnl}^{(-1)}(c_v) + \frac{1}{(\xi ^2-1)} [ 2 q_{mnl}^{(1)}(c_v) - q_{mnl}^{(3)}(c_v) ] \Big ) \nonumber \\ & & -\,\, q_{mnl}^{(-1)}(c_v) - \frac{1}{(\xi ^2-1)} [ 2q_{mnl}^{(1)}(c_v)+3\hat{q}_{mnl}^{^{\prime }(1)}(c_v)] \nonumber \\ & & -\,\, \frac{1}{(\xi ^2-1)^2} [q_{mnl}^{(3)}(c_v) + \hat{q}_{mnl}^{^{\prime }(3)}(c_v)] \bigg \rbrace \frac{R_{mn}^{(j)}(c_v,\xi )}{c_v\xi } \nonumber \\ & & +\,\, \bigg \lbrace -(\xi ^2+2) q_{mnl}^{(-1)}(c_v) + \frac{(2\xi ^2-3)}{(\xi ^2-1)} q_{mnl}^{(1)}{c_v} \nonumber \\ & & +\,\, \hat{q}_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} [q_{mnl}^{(3)}{c_v} + \hat{q}_{mnl}^{^{\prime }(3)}(c_v)] \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ] . \end{eqnarray} (C7b)  \begin{eqnarray} F_{\phi mnl}^{(Mj)}(v) & = & \mu i^n \Bigg [ - \frac{1}{(\xi ^2-1)} \bigg \lbrace \Big ( \Lambda _{mn}(c_v) - c_v^2 \xi ^2 \Big ) \Big ( \hat{q}_{mnl}^{(1)}(c_v) + \frac{1}{(\xi ^2-1)} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, m^2 \Big ( \hat{q}_{mnl}^{(-1)}(c_v) + \frac{3}{(\xi ^2-1)} \hat{q}_{mnl}^{(1)}(c_v) + \frac{2}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}(c_v) \Big ) \nonumber \\ & & +\,\, (\xi ^2+1) q_{mnl}^{^{\prime }(1)}(c_v) - q_{mnl}^{^{\prime }(3)}(c_v) \bigg \rbrace R_{mn}^{(j)}(c_v,\xi ) \nonumber \\ & & +\,\, c_v \xi \bigg \lbrace q_{mnl}^{^{\prime }(1)}(c_v) + \frac{1}{(\xi ^2-1)} [ 4\hat{q}_{mnl}^{(1)}{c_v}+q_{mnl}^{^{\prime }(3)}(c_v) ] \nonumber \\ & & +\,\, \frac{2}{(\xi ^2-1)^2} \hat{q}_{mnl}^{(3)}{c_v} \bigg \rbrace \frac{R_{mn,\xi }^{(j)}(c_v,\xi )}{c_v} \Bigg ]. \end{eqnarray} (C7c) Similar to the $$q_{mnl}^{(k)}$$ functions of Appendix B, many of the $$E_{.mnl}^{(.j)}$$ and $$F_{.mnl}^{(.j)}$$ functions are zero. The non-zero values are for (n − m + l) even   $$E_{\eta mnl}^{(Mj)}, \, E_{\xi mnl}^{(Lj)}, \, E_{\xi mnl}^{(Nj)}, \, E_{\phi mnl}^{(Lj)}, \, E_{\phi mnl}^{(Nj)}, \, F_{\eta mnl}^{(Mj)}, \, F_{\xi mnl}^{(Lj)}, \, F_{\xi mnl}^{(Nj)}, \, F_{\phi mnl}^{(Lj)}, \, F_{\phi mnl}^{(Nj)} ,$$ (C8a)and for (n − m + l) odd   $$E_{\eta mnl}^{(Lj)}, \, E_{\eta mnl}^{(Nj)}, \, E_{\xi mnl}^{(Mj)}, \, E_{\phi mnl}^{(Mj)}, \, F_{\eta mnl}^{(Lj)}, \, F_{\eta mnl}^{(Nj)}, \, F_{\xi mnl}^{(Mj)}, \, F_{\phi mnl}^{(Mj)} .$$ (C8b) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### DeepDyve Freelancer ### DeepDyve Pro Price FREE$49/month

\$360/year
Save searches from Google Scholar, PubMed
Create lists to organize your research
Export lists, citations