Solution (source code)

= Solution

In the lower uniform <half-space> the <causal vertical wavenumber branch> has $\operatorname{Im}k_{\beta,b}>0$. The downward factor $e^{ik_{\beta,b}(z-z_b)}$ decays as $z\to+\infty$, whereas the upward factor grows there and represents an incoming wave from the lower infinity. The outgoing <radiation condition> therefore retains only the downward factor and gives
$$
\boxed{Z(z_b)=q_b=\mu_bp_{\beta,b}.}
$$
For each uniform layer $j$ of thickness $h_j$, work upwards from its known bottom load $Z_{j,b}$. With $q_j=\mu_jp_{\beta,j}$ and $\delta_j=\omega p_{\beta,j}h_j$, part (iv) gives
$$
\boxed{Z_{j,t}=\frac{Z_{j,b}\cos\delta_j-iq_j\sin\delta_j}
{\cos\delta_j-i(Z_{j,b}/q_j)\sin\delta_j}.}
$$
<velocity> and <traction> are continuous across a bonded interface, so this top impedance becomes the bottom load of the next shallower layer. Repeat up to $z_a$ to obtain $Z_{\rm in}=Z(z_a)$. Poles can be handled by carrying $(V,T)$ or the reciprocal ratio. Since the lower load has positive real part, part (iii) actually ensures that the outgoing whole-stack impedance stays finite with positive real part at all shallower depths in the upper <frequency> half-plane.

The upper <half-space> has directional impedance $q_a=\mu_ap_{\beta,a}$. At $z_a$, an incident <velocity> $V_0$ and reflected <velocity> $RV_0$ give
$$
V(z_a)=V_0(1+R),\qquad T(z_a)=-q_aV_0(1-R).
$$
Taking their ratio and solving establishes <reflection from a layered SH impedance>:
$$
Z_{\rm in}=q_a\frac{1-R}{1+R},\qquad
\boxed{R=\frac{q_a-Z_{\rm in}}{q_a+Z_{\rm in}}.}
$$
Thus no individual deeper reflection coefficients are needed after $Z_{\rm in}$ is known: it contains all underlying interfaces and their coherent multiples. A matched load gives $R=0$; the real-frequency limits of zero <traction> and zero <velocity> give $R=+1$ and $R=-1$. For real <frequency> and propagating upper and lower <half-space> waves, the usual energy-flux balance provides a further check on the computed reflection and transmission. A stable implementation for thick layers uses the ratio update or its decaying-exponential version instead of multiplying growing transfer <matrices>.