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 (*x _{1},x_{2},z*)
or (

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
**R*** _{U}*,

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 **T*** _{U}*,

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,`

*p*(complex*16 values)]_{x}, p_{y}

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 **T*** _{D}*, ...,
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).

In: Seismic Waves in Complex 3-D Structures, Report 7, pp. 163-167, Dep. Geophys., Charles Univ., Prague, 1998.

SW3D - main page of consortium