Copyright © 1996 by C. J. Thomson, Queen's University
Please Note: This document on the program Rmatrix is accompanied by some notes on the basic theory of waves in layered media. Familiarity with the equations will assist the user with notation, which is rather similar in the FORTRAN code. Realistically, the theory needs to be understood in order to use the program correctly. The code is written with the same cartesian coordinates (x1,x2,z) or (x,y,z) in mind, with z positive downwards. SI units are used.
These notes are not exhaustive and they will grow and become more useful as users' questions are answered.
Rmatrix is a subroutine which returns the reflection and transmission matrices RU, RD, TU and TD for a stack of arbitrarily anisotropic homogeneous plane layers (see theory notes). These matrices embody everything about waves in the stack and may be used to model all kinds of propagation problems. For example, if many thin layers are used to simulate a gradient one may study back-scattered waves from the gradient or quasi-shear wave forward coupling. One can model plane waves incident from below a crustal/upper-mantle structure (e.g. the SKS and receiver function techniques). Surface waves in anisotropic media can be explored.
The basic call to rmatrix is
call rmatrix(px,py,omega,nlay,ifree,tus,rus,tds,rds,iter)
and for the given slownesses and frequency it returns the 3x3 reflection/transmission matrices TU, RU, TD and RD for a stack of anisotropic layers between two half spaces. The variables ifree and iter are the least obvious. Normally ifree is zero and then the values of TU and RD are those just above the first interface in the model (i.e. at the bottom of the first layer, so the thickness of that layer is irrelevant). When ifree is nonzero the phase shift across the top layer is included just before rmatrix is exited (useful when free-surface effects are to be included at the top of the first layer). The bottom layer of the model is taken as a lower half space, so its thickness is always irrelevant. When iter>1, the frequency-independent R/T coffecients at each single interface are not recalculated, but are assumed to have been stored for the given (px,py) in a previous call to rmatrix with iter=1.
Once the R/T matrices are found one may use them to solve various initial/boundary-value problems. For example, the propagator matrix for displacement/stress vectors can be constructed from equations (4.2) in the theory notes. Subroutines to add, multiply and invert 3x3 matrices are part of the program and these can be exploited.
Important point: Care should be taken with the ordering of the eigenvalues and eigenvectors. In the isotropic case the six columns are easily ordered into P-down, SV-down, SH-down, P-up, SV-up and SH-up (see subroutine isotroc). In anisotropic media the polarisations can sometimes behave in a complicated way, so little effort has been expended yet to provide a general sorting routine. The ordering currently depends on which (numerical) method is used to find the vertical slownesses (i.e. eigenvalues -- see comments in subroutine evset). If the most general routine is used (aniso6cg) or Georg Ruempker's routines (called from aniso3c) are used, the vertical slownesses are ordered in decreasing size (of real or imaginary part as appropriate), so the columns of the eigenvector matrix are possibly like P-down, S1-down, S2-down, S2-up, S1-up and P-up. If Sean Guest's routines are used (via aniso6c1) the columns are ordered more like that in the isotropic case. The user can order the columns as desired, so long as the first three correspond to downward propagating or decaying waves and the last three upward propagating or decaying waves. Awareness of the eigenvector orderings is important for the correct interpretation of the R/T coefficents, ray expansions, etc., and for anisotropic media the investigator will no doubt be looking closely at these quantities and the polarisations.
The program is divided up into 8 files. The first file (rmatrix1.f) has a main program, which calls rmatrix, two input subroutines and rmatrix itself. The second file (rmatrix2.f) contains the core routines of the reflection matrix calculation and its recursion scheme. The next 4 files (isotroc.f, aniso6cg.f, aniso6c1.f & aniso3c.f) contain subroutines to find the eigenvalues and vectors for a given layer. The user may supply personal versions of these if desired (e.g. efficient ones specific to, say, orthorhombic media). Finally there are two files (cg.f & deigv.f) containing public-domain matrix eigen-problem software, which can also be replaced at will by the user.
Most of the variables in the input file (rmatrix.inp) are obvious. The first two lines are like
5 1 [nlay, ifree]
(1.d0,0.d0) (0.5D-5,0.d0) (0.5d-5,0.d0)
[omega, px, py (complex*16 values)]
For an isotropic layer there are two lines of the input file such as
2 4.0d3 1000.0000000000
1.5d3 0.0d3
Many of the common blocks are self-explanatory. For example
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),rho(nlaymx)
contains the Earth model and
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)
contains the eigenvalues and vectors in the layers for the current value of slowness.
For each buried interface in the model the 3x3 blocks of the local wave matrix Q and the local R/T matrices TD, ..., are stored in
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
etc
and the cummulative coefficients from the bottom of the stack up to this interface are stored in
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
etc.
Some quantities used in symmetry testing are stored in
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
Certain variables used for group velocity calculations are found for forward and reversed modes and stored in
common /groupvel/vece(6,6,nlaymx,2),vale(6,nlaymx,2)
common /grouprt/tdslay(3,3,nlaymx,2),rdslay(3,3,nlaymx,2),cdsto(3,2)
The source and receiver depths are stored in
common /srcrec/zsrc,zrec,ifmode
The programs can be compiled using the standard f77 command under Unix. The code is self contained and no external libraries are needed. It has been developed on a Sun system, where the command
f77 -o rmatrix rmatrix1.f rmatrix2.f isotroc.f aniso6cg.f
aniso6c1.f aniso3c.f deigv.f cg.f -C
has been used.
Only one output file is produced (rmatrix.out) and by changing various ``iwrite'' variables throughout the program greater or lesser amounts of information can be dumped (e.g. eigenvalues and vectors in the layers).
There are three types of symmetry which can be tested. This is best discussed once the theory is understood.
The program can be obtained via anonymous ftp and information on how to get it can be obtained by sending email to thomson@geol.queensu.ca.
The paper is available in PostScript (77 kB) and GZIPped PostScript (21 kB).
SW3D
- main page of consortium Seismic Waves in Complex 3-D Structures .