When writing code, it is important to test it.  For moment tensor inversion this can be accomplished by creating a set of synthetic seismograms for a given moment tensor and then determining if the inversion codes yield the known moment tensor. 

For routine moment tensor inversion of earthquake, I prefer using wvfgrd96 (which is equivalent to wvfmtgrd96 -DC) because both have an efficient  algorithm to handle the time shifts necessary to align the observed and predicted waveforms  because of location or model error. The deviatoric and full moment tensor inversion codes wvfmt96 and wvfmtd96 are fast but do not have a good way to handle time shifts. On the other hand, the wvfmtgrd, wvfmtgrd -DC and wvfmtgrd -DEV  handle the time shifts well by using the same algorithm used by wvfgrd96. The wvfmtgrd and wvfmtgrd -DEV run more slowly because the grid search is over 5 and four parameters, respectively. For  the noise free case, the examples below show that  the wvfmt96 and wvfmtgrd96 give the same results as do the wvfmtd96 amd wvfmtgrd96 -DEV. Until someone rewrites the code for wvfmt96 and wvfmtd96 so that time shifts are handled well, I recommend the use of wvfmtgrd -DC and wvfmtgrd -DEV.

This tutorial presents the results of processing moment tensor inversion codes for different synthetic data sets. For each data set, an inversion is performed forthe


To repeat the computations download the file MTIVTEST.tgz. Unpack this with the command

gunzip -c MTIVTEST.tgz | tar xf -

This will create the following directory structure:

├── 0XXXREG          # these are the prototype directories with scripts are used for all tests         
├── 20210813115735   # This is a real data set. The dsta in distribution is used in the synthetic test
│   ├── DAT.REG
│   │   └── NOUSE
│   ├── GRD.REG
│   ├── HTML.REG
│   ├── MAP.REG
│   ├── MLG.REG
│   ├── ML.REG
│   ├── MTD.REG
│   ├── MTGRD.REG
│   ├── MTGRD.REG.DC
│   ├── MT.OTHER
│   ├── MT.REG
│   ├── NEW2.REG
│   └── SYN.REG
├── 20210813115735_0.06  # this has the same structure as above, but the lower frequencies are used for the inversion  
│   ├── ...
├── bin  # executables that are specific to moment tensor inversion
|── MTDC # organization for testing the moment tensor inversion of a double couple source │   ├── DAT.REG # waveforms for inversion │   ├── DAT.SYN # create synthetic data set and place into DAT.REG │   ├── GRD.REG │   ├── HTML.REG │   ├── MLG.REG │   ├── ML.REG │   ├── MTDC │   ├── MTD.REG │   ├── MTGRD.REG │   ├── MTGRD.REG.DC │   ├── MTGRD.REG.DEV │   ├── MT.OTHER │   └── MT.REG ├── MTDEV # organization for testing the moment tensor inversion of a deviatoric source │   ├── ... (see above) ├── MTFULL # organization for testing the moment tensor inversion of a full moment tensor source
├── ... (see above)
├── MTFULL.NOISE # organization for testing the moment tensor inversion of a full moment tensor source with noise
├── ... (see above) ├── src # source code for executables to be placed in the ../bin directory └── TGREEN # precomputed Green's functions at only those distances for the test. └── CUS.REG ├── 0010 ├── 0020 ├── ... ├── 0270 ├── 0280 └── 0290      

After downloading, set the environment to point to the CPS bin directory and to the CPS source codes. Here are some examples of how to do this.

LINUX example 1
In ~/.profile
# set PATH so it includes user's PROGRAMS.330/bin if it exists
if [ -d "$HOME/PROGRAMS.330/bin" ] ; then
In ~/.bashrc
export CPS=${HOME}/PROGRAMS.330
LINUX example 2
In ~/.profile
In ~/.bashrc
export CPS='/home/rbh/PROGRAMS.310t/PROGRAMS.330/'
OSX example
In ~/.profile

It is also necessary to point to the Green's functions before trying to duplicate this tutorial. For this tutorial, the required subset of  Green's functions is provided.

export GREENDIR=`pwd`

Finally compile the executables:

rm -fr bin
mkdir bin
cd src
make all
cd ..


To test the codes and to compare results the analyses were computed in each of the following directories.The link at the end of each link points to the documentation for each data set.

Synthetic data sets

The first thing that is required is a station distribution, which is defined here in the file list1 which has the following columns entries: Station_name, Network_name, epicentral distance (km), azimuth from epicenter to the station (degrees), station latitude and station longitude.

BLO   NM 392.847  339.053 39.1719 -86.5222
CASEE CO 203.703  118.211 34.993  -82.9317
CPCT  ET  58.4465 144.044 35.45   -84.52
GOGA  US 303.381  153.962 33.4112 -83.4666

The next thing required is a precomputed set of Green's functions. This is described in the tutorial There is a directory that contains the Green's functions and an environment parameter GREENDIR that points tot he directory. If one does a ls ${GREENDIR}, one may see

Models/      CUS.REG
The Models directory will have the velocity model in the model96  format, e.g.,  CUS.mod. The CUS.REG directory has the Green's functions computed for various source depths and distances. The directory names indicate the source depth. These names are of the form DDDd which represents a depth of DDD.d km.

0005/  0030/  0060/  0090/  0120/  0150/  0180/ 
0010/  0040/  0070/  0100/  0130/  0160/  0190/ 
0020/  0050/  0080/  0110/  0140/  0170/  0200/ 

The subdirectories contain the Green's functions for a given source depth. Here 0150 represents a source depth of 15.0 km.

Finally each depth directory has Green's functions.

009800100.RDD  009800100.RDS  009800100.REX  009800100.RSS 
009800100.TDS  009800100.TSS  009800100.ZDD  009800100.ZDS     
009800100.ZEX  009800100.ZSS  W.CTL
  gives synthetics at a distance of 0098.0 km for a source depth of 010.0 km. The W.CTL file is uses as an index. A few lines of which are

 94 0.25 512 6.75  0 0100 009400100
 96 0.25 512 7     0 0100 009600100
 98 0.25 512 7.25  0 0100 009800100
100 0.25 512 7.5   0 0100 010000100
105 0.25 512 8.125 0 0100 010500100

which gives the distance, sample interval, number of data points the t0 and vred to compute the first time sample, the depth directory and the corresponding Green's function prototype. The shell script DOSTA  below will search for the Green's function distance closed to the observed data. Thus if the epicentral distance is 97.5 km, the script will select 009800100 and thus the Green's functions listed above.

There are two steps: make the synthetics for each station in list1 and then, optionally, add noise. Within each source inversion directory, there is a DAT.SYN directory with a script that reads the file list1. Although the scripts start with the initial full moment tensor, some of the scripts will create synthetics for the full moment tensor, for the deviatoric component of the moment tensor and for the major double couple component of the moment tensor.


This shell script to make synthetics for the full moment tensor is as follows:


# make synthetics for a given moment tensor and mw

#    define the origin time # this is so that the observed and synthetics can be plotted on the same absolute time scale 


HS=0010                     # This is important since it defines the source depth  of 1.0km for the synthetics. This naming matches the
# directory structure in ${GREENDIR}/CUS.REG/
mtinfo -XX $MXX -XY $MXY -XZ $MXZ -YY $MYY -YZ $MYZ -ZZ $MZZ -a > mtinfo.txt GREEN=${GREENDIR}/CUS.REG # Change this to use Greens functions if another model is used. MOMENT=`echo $MXX $MXY $MXZ $MYY $MYZ $MZZ| awk '{print sqrt(0.5*($1*$1 + 2*$2*$2 + 2*$3*$3 + $4*$4 + 2*$5*$5 + $6*$6))}' `
echo $MOMENT
fmplot -FMPLMN -P -XX $MXX -XY $MXY -XZ $MXZ -YY $MYY -YZ $MYZ -ZZ $MZZ ##### # begin the computation of synthetics ##### while read STA NET DIST AZ STLA STLO do echo $STA $NET $DIST $AZ ##### # search over source depth These depths are the subdirectory # names in the Green's Function Directory ##### cat > awkprog << FOE # This works under gawk - on Solaris try nawk BEGIN { MDIF = 10000.0 } {DIF = $DIST - \$1 ; if( DIF < 0 ) DIF = - DIF ; if(DIF < MDIF) { MDIF = DIF ; Dfile = \$7 ; Rate = \$2 ; Dist = \$1 } } END { print Dfile , Rate, Dist } FOE cat ${GREEN}/${HS}/W.CTL | \ awk -f awkprog > j DFILE=`awk '{print $1}' < j ` PROTO=${GREEN}/${HS}/${DFILE} echo $STA $NET $DIST $AZ $PROTO gsac << EOF mt to ZRT MXX $MXX MXY $MXY MXZ $MXZ MYY $MYY MYZ $MYZ MZZ $MZZ AZ $AZ FILE ${PROTO} w rh T.? ch NZYEAR $NZYEAR NZJDAY $NZJDAY NZHOUR $NZHOUR NZMIN $NZMIN NZSEC $NZSEC NZMSEC $NZMSEC ch ocal $NZYEAR $NZMON $NZDAY $NZHOUR $NZMIN $NZSEC $NZMSEC ch KSTNM $STA KNETWK $NET STLA $STLA STLO $STLO ch lcalda false ch evla ${EVLA} evlo ${EVLO} wh quit EOF mv T.Z ../DAT.REG/${STA}${NET}BHZ mv T.R ../DAT.REG/${STA}${NET}BHR mv T.T ../DAT.REG/${STA}${NET}BHT done < list1 rm -f j awkprog ##### # create the true solution for the comparison panel ##### mtinfo -xx $MXX -yy $MYY -zz $MZZ -xy $MXY -xz $MXZ -yz $MYZ > mtinfo.out cp mt.msg ../MT.OTHER/true

The script is commented. The upper part defines the origin time and location of the event as well as its moment tensor values. The source depth is given as 0010 to agree with the organization of the Green's function directory. If the Green's functions were ground velocity in cm/s (default when using KM and GM/CM^3 for model units), then the output here will be ground velocity in m/s.The gsac mt command make the synthetics for the given moment tensor and azimuth using the Green's functions and also ensures the final units.s

The resulting synthetics can be compared to observed data since they have the same start and end times, distances and azimuths.

The script for the deviatoric moment tensor synthetics has these lines after the moment tensor is defined:

# make deviatoric
MI=`echo $MXX $MYY $MZZ | awk '{printf "%11.3e", ($1+$2+$3)/3.}' `
MXX=`echo $MXX $MI | awk '{printf "%11.3e", $1 - $2}' `
MYY=`echo $MYY $MI | awk '{printf "%11.3e", $1 - $2}' `
MZZ=`echo $MZZ $MI | awk '{printf "%11.3e", $1 - $2}' `

The deviatoric moment tensor is formed by subtracting the isotropic component from the moment tensor.  Thus if the moment tensor is changed, then this part of the code will give the correct deviatoric moment tensor.

To make synthetics for the double couple source, the strike, dip and rake angles can be specified with the moment magnitude. However, I chose to manually run mtinfo to determine the moment tensor corresponding to the major double couple. This was then manually inserted into the DOSYN script in MTDC/DAT.SYN. Thus that script has

#These are the from the major double couple of the original
#moment tensor which is obtained from mtinfo



This script DOSYN script in MTFULL.NOISE/DAT.SYN is given here. This script makes use fo the program sacnoise which was part of the tutorial Moment Tensor Sensitivity to Noise. This code is in the disctributed src directory.


export PATH=:../../bin:$PATH

# add noise to existing file

PVAL=0.4                # PVAL=0 is NLNM noise while PVAL=1.0 is NHNM noise

#    add noise to each file in ../DAT.REG
echo ========= FOR LOOP ===========
for TRACE in TMP/*[ZRT]
        cp $TRACE T.sac
        gsac << EOF
rh T.sac
synchronize B
        #  different random seed for each waveform
        echo ======================= $i =======================
        #    get information about the file
        NPTS=`saclhdr -NPTS $i`
        DELTA=`saclhdr -DELTA $i`

        NPTS=`saclhdr -NPTS ${i}`
        B=`saclhdr -B ${i}`
        E=`saclhdr -E ${i}`
        A=`saclhdr -A ${i}`
        O=`saclhdr -O ${i}`
        DELTA=`saclhdr -DELTA ${i}`
        DIST=`saclhdr -DIST ${i}`
        AZ=`saclhdr -AZ ${i}`
        BAZ=`saclhdr -BAZ ${i}`
        EVLA=`saclhdr -EVLA ${i}`
        EVLO=`saclhdr -EVLO ${i}`
        EVDP=`saclhdr -EVDP ${i}`
        STLA=`saclhdr -STLA ${i}`
        STLO=`saclhdr -STLO ${i}`
        KCMPNM=`saclhdr -KCMPNM ${i}`
        KSTNM=`saclhdr -KSTNM $i`
        KNETWK=`saclhdr -KNETWK $i`
        NZYEAR=`saclhdr -NZYEAR ${i}`
        NZMON=`saclhdr -NZMON ${i}`
        NZDAY=`saclhdr -NZDAY ${i}`
        NZJDAY=`saclhdr -NZJDAY ${i}`
        NZHOUR=`saclhdr -NZHOUR ${i}`
        NZMIN=`saclhdr -NZMIN ${i}`
        NZSEC=`saclhdr -NZSEC ${i}`
        NZMSEC=`saclhdr -NZMSEC ${i}`
        sacnoise -dt ${DELTA} -s ${RVAL} -p ${PVAL} -npts ${NPTS}
        #  sacnoise creates O.sac
        #  To get noise before the synthetic
        #    for the synthetic. The synthetic starts at zero
        #    use shift to set the correct B time

        #    then set the time stamp
        #    then synchronize O
        #    then set the A time for the P first arrival

        gsac << EOF

echo $B
r O.sac
ch lcalda false
shift Fixed ${B}
cut ${B} ${E}
cut off
        ch a $A o $O
        ch KSTNM ${KSTNM}
        ch KNETWK ${KNETWK}
        ch KCMPNM ${KCMPNM}
        ch DIST $DIST AZ $AZ BAZ $BAZ
        # the output of sacnoise is an acceleration history in M/S^2
        transfer from none to none freqlimits 0.005 0.01 10 20
        #convert to velocity
        w noise
        rh noise
        synchronize o
#    now add noise tot he original file
        r $TRACE noise
        addf master 1
        w 1 2
B=`basename ${TRACE}`
mv 1 ../DAT.REG/${B}
rm 2

The PFAC=0.4 sets the noise level. 0.0 is for the NLNM and 1.0 is for the NHNM. most of the script is dedicated to making the noise start and end at the same time as the noise-free time series, so that they can be summed with the result having the correct headers.

The result of this simulation is shown in the next figure:

Time series: Noise free (red) and with noise (blue)

Fourier velocity spectra: Noise free (red) and with noise (blue)

Last Changed May 30, 2022