Implementation

The problem is to find the solutions of the equation $l_{21}(V_L)=0$ . The computation of the element $l_{21}(z_0)$ requires the multiplication of $n$ 2x2 matrices, whereas three of the four elements of the final matrix $L(z_0)$ are useless. From the formulae of the product of two 2x2 matrices, only $l_{21}(z_1)$ and $l_{22}(z_1)$ are necessary for computing $l_{21}(z_0)$ which, in turn, requires $l_{21}(z_2)$ and $l_{22}(z_2)$ (equation 3.21).

$$\underbrace{\left(\begin{array}{cc} - & - \ l_{21}(z_0) & - \end{array}\right)}_{L(z_0)} = T_n^{-1} G_n \ldots{} G_1$$ (3.21)

$$= \underbrace{\left(\begin{array}{cc} - & - \ l_{21}(z_n) & l_{22}(z_n)\end{array}\right)}_{L(z_n)=T_n^{-1}} G_n \ldots{} G_1$$

$$\ldots{}$$

$$= \underbrace{\left(\begin{array}{cc} - & - \ l_{21}(z_2) & l_{22}(z_2)\end{array}\right)}_{L(z_2)}G_2 G_1$$

$$= \underbrace{\left(\begin{array}{cc} - & - \ l_{21}(z_1) & l_{22}(z_1)\end{array}\right)}_{L(z_1)}G_1$$

To calculate $L(z_n)=T_n^{-1}$ , the factors $e^{-\nu_n z}$ are dropped from elements $l_{21}(z_n)$ and $l_{22}(z_n)$ (equation (3.16) because we are looking for the roots of $l_{21}(V_L)$ .

From equation (3.17), $\nu_n$ is imaginary if $k$ is less than the wavenumber of S-waves $k_{sn}=\frac{\omega^2}{V_{sn}^2}$ . To avoid using complex number libraries, the $sinh$ and $cosh$ functions of equation (3.19) are replaced by the corresponding trigonometric functions $sin$ and $cos$ . For real values of $\nu_n$ , the hyperbolic functions do not tolerate high arguments. They are preferably computed from their analytical formulae:

$$sinh(ix)=\frac{e^{ix} - e^{-ix}}{2}= e^x \frac{1- e^{-2x}}{2}$$ (3.22)

$$cosh(ix)=\frac{e^{ix} + e^{-ix}}{2}= e^x \frac{1+e^{-2x}}{2}$$

Hence, an exponential factor can be dropped from the expression of $G$ in equation (3.19) because we are seeking for roots. When factor $e^x$ is dropped, the computation of both hyperbolic functions require the calculation of only one exponential function. As we are working in double precision (floating points of 64 bits), all exponential values $e^{-2x}$ less than $10^{-19}$ are equivalent to zeros in equations (3.22). Thus, in equation 3.19, the hyperbolic functions reduce to $\frac{1}{2}$ if

$$\nu_{n-1} (z_n-z_{n-1})>\frac{19}{2}*ln(10)\approx21.2$$ (3.23)

In each layer, values of $l_{21}(z_i)$ are scaled to fit in the range between $-10^5$ to $10^5$ , to avoid overflow when propagating across a stack with many layers.

In figure 3.2, the values taken by $l_{21}(z_0)$ for all couples $(\omega,k)$ are shown in the case of a two-layer model: 200 m/s for $V_s$ in the first 25 m thick layer, and 1000 m/s in the half-space. The first 12 modes located at the root of the function are highlighted by black lines. The negative values of the function are not represented (white areas). The normal modes of Love are observed between the minimum and the maximum $V_s$ of the model. The fundamental mode is present over the whole frequency range, whereas each higher mode has its own threshold frequency under which it does not exist.

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

Finding the roots of $l_{21}(z_0)$ is not straightforward. This issue is treated in section 3.1.5 together with the overall performances of this algorithm.