Implementation

The detailed expressions of the determinants of $R(z_i)$ are given in appendix A. The six-component vector $R(z_{n-2})$ is obtained by combining the matrix $G_{n-1}$ and the vector $R(z_{n-1})$ in the same way as equations (A.12) for $G_n$ and $T_n^{-1}$ . The computation is done for each layer for bottom to top up to the first one at the free surface. As $T_{1214}$ and $T_{1223}$ (appendix A) are equal (from equation (3.31)), it is obvious that $R_{1223}(z_{n-1}) = R_{1214}(z_{n-1})$ . From bottom to top, the two components are always equal (equation (3.30)) and we can reduce the number of components to five rather than six. Also, $R_{1214}(z_{n-1})$ like $T_{1214}$ is the only imaginary component and this feature is preserved across the layered medium. Thus, the 6 components of $t_n^{-1}$ reduce to 5 and the matrix $g_n$ to 5x5 components.

To speed up the computation, we slightly modified the Dunkin's original formulation to reduce the total number of operations, preferring subtractions, additions and multiplications to divisions. The $sinh$ and $cosh$ functions are calculated in the same way as for Love case, including the real and imaginary cases (equation (3.22), section 3.1.3). A frequency factor of $\omega^2$ (equation (A.12), appendix A) has been introduced in $R_{1212}$ to avoid unscaled vector at low frequencies. For each layer, values of sub-determinants are scaled to fit in the range between $-10^5$ to $10^5$ to avoid overflow when propagating across a stack with many layers. Compared to Herrmann's formulation (1994) in the same conditions (not in its original Fortran code but already translated in C++), this implementation reduces by 25% the time consumption.

Figure 3.3: Values taken by $R_{1212}(z_0)$ at different couples (frequency, velocity) for a two-layer model.

In figure 3.3, the values taken by $R_{1212}(z_0)$ for all couples $(\omega,k)$ are shown in the case of a two-layer model: 1350 and 250 m/s for $V_p$ and $V_s$ , respectively in the first 25 m thick layer, and 2000 and 1000 m/s in the half-space.