Rongjiang Wang
A simple orthonormalization method for stable and efficient computation of Green's function
Bulletin of the Seismological Society of America (June 1999), 89(3):733-741

Abstract:
Some major matrix methods for computation of Green's functions of a layered half-space model are compared. It is known that the original Thomson-Haskell propagator algorithm has the loss-of-precision problem when waves become evanescent. The minor matrix extension solves to some extent the numerical problem but at the expense of the computational efficiency. At present, the recursive technique (the reflectivity method) is used by most seismologists. An alternative algorithm is presented that uses the same strategy as the original propagator algorithm but avoids the numerical problem by inserting an additional numerical procedure into the matrix propagation loop. The new technique is simple and efficient, and in particular, it is not only applicable to the matrix method but also to any numerical integration method for general use. High-frequency synthetic seismograms are computed for modeling significant site effects such as observed in the aftershocks of the Latur earthquake, 1993, India.