Waveform Modeling and Moment Tensor Inversion Tools

Nevada Seismological Laboratory, University of Nevada

Feburary 6, 2000

Updated Feburary 27, 2000

Introduction

Waveform Modeling Tools are a set of tools which compute Green's functions for a multilayered halfspace using f-K integration/reflectivity method, and then estimate earthquake parameters (fault strike, dip, rake, and seismic moment) by fitting synthetic ground displacement to observed seismograms. A grid search and moment tensor inversion code are also part of the package. These programs are coded in the C computer language and should be portable to all Unix platforms. Additional programs and libraries are needed to run this packages: Seismic Analysis Code, and Pgplot. These will need to be downloaded and installed before compiling the source code. Binary executables are also supplied for LinuxOS distribution. These tools run in the Unix environment and use C-shell scripts or command line interface. Postscript and GIF hardcopies can be obtained through the Pgplot libraries and are preprogrammed into the tools.

Sample output

Plot 1 Three component ground displacement (cm) overlayed with synthetics. A layered halfspace can be used to model long period ground motion very well.

Plot 2 Gridsearch least square residuals versus model parameters: fault strike, dip, rake, and moment magnitude. Also residuals for simulated annealing downhill simplex program.

Plot 3 A least squares moment tensor inversion of the same event.

Download Distribution Here...

Intel processor LinuxOS linux.tar.gz - Gzipped Tar file 6654687 bytes

Sun Microsystems SolarisOS sun.tar.gz - Gzipped Tar file 11306068 bytes

Installation and Running an Example

There is a sample directory with data and ready to run shell scripts in the distribution directory ./distribution/FishLakeValley/19971102085153. Here is another example with general instructions on how to download data and run all of the programs from scratch. Modeling tips have a blue background.

Install SAC.
Download the Waveform Tools from above.

uncompress and untar the file.

gzip -d linux.tar.gz 
tar xf linux.tar

or on Solaris or SunOS

gzip -d sun.tar.gz
tar xf sun.tar

Make the executables
```
cd ./distribution
make clean
make all
```
If the compilation fails due to the PGPLOT libraries, then you will need to download, compile and install the libraries from: http://phobos.caltech.edu/~tjp/pgplot/. Copy the new libraries libcpgplot.a and libpgplot.a into the ./lib directory. Or if they are installed in /usr/lib then change the src/Makefile file to make the FPGPLOTLIBDIR variable point to this directory.
Make a directory named Event and then make a directory say 01.
```
mkdir Event
mkdir Event/01
cd ./Event/01
```
Download or copy some data into the ./Event/01 directory.
Convert it to sac format. If it is in sac format make sure it is the correct byte order.
Big endian to little endian (SUN to Intel processor)
```
sacsun2linux *.SAC
```
Little endian to big endian (Intel processor to SUN processor)
```
saclinux2sun *.SAC
```
If the sac files are made from the same platform as where the processing is made then no conversion is needed. Both programs accept wild cards and old files will be written over.
Make sure the file names are in convention with RDSEED because the C-shell script PROCESS_MACRO only reads files with this naming convention. If the names are different then you may change the C-shell script PROCESS_MACRO or form to this convention by renaming the files. Here is the RDSEED naming convention
```
1997.306.08.40.57.0655.BK.KCC..BHE.SAC
```
where the fields are separated by periods
```
sprintf(sacfile, "%4d.%03d.%02d.%02d.%02d.%04d.%-2.2s.%-4.4s.%-3.3s.SAC", 
```
If you do not already have data then go to the IRIS Wilbur Web site and download some data. Pick a recent event in the western US (Mag > 4.5) with seismograms from at least one of the following stations ANMO, NEW, TUC, COR, DUG, CMB, and PAS. Download the data as a RDSEED file and then use the command
```
rdseed -f WLBR_SPYD_990801_1606.seed -d
```
to extract sac files. If you do not have RDSEED then use the sac binary option but make sure its in the correct byte order.

make a file called runall.csh using your favorite text editor.

vi runall.sh

Say we use the following event from the Wilbur web site.

1999/08/01 16:06:22.0 SPYDERŽ 5.6  37.39 -117.07    5.00 CALIFORNIA-NEVADA BORDER REGION

then type this in the file runall.sh.

#!/bin/csh
# An example of the Aug 1, 1999 Scottys Junction, Nevada Earthquake.

############################################################################################################
#             sta  lc   hc   nt  dt   vel elat   elon   z   yr   jday hr mn sec  database file   resp dir  #
############################################################################################################
PROCESS_MACRO tuc  0.01 0.05 256 1.75 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO dug  0.01 0.05 256 1.00 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO cmb  0.01 0.05 256 0.80 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO pas  0.01 0.05 256 1.50 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO anmo 0.01 0.05 256 1.75 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO new  0.01 0.05 256 2.00 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response

This C-shell script can be cloned and edited for different parameters and events and will remain as a usefull record of the parameters used. Use the # symbol to make comments. See the man page PROCESS_MACRO and MAKE_SYN for the input format.

The source location and origin time can be typed in from the NEIC catalog location from the Wilbur web site. Use the best location information you can find for your events. We incorporate a shift between the data and synthetics for the maximum waveform correlation before the L2 raw amplitude misfit is measured which will take into account any errors in source location and origin time. Others have shown that greens functions do not vary much as a function of distance so that this is a reasonable procedure.

If you have problems downloading or converting the files then you can ftp the file from my website. WLBR_SPYD_990801_1606.SAC.tar.gz (Big Endian byte order)

You can use any event anywhere in the world but you must include the station location in the file ../../stalist and supply instrument responses in the directory ../../Response.

Start easy with a large (M 5) earthquake recorded by many stations and filter the data around 50 to 100 seconds period. Model the long period waves using a simple velocity model and then find a good solution with an optimal depth. Make sure the solution agrees with some ground truth. Then decrease the period by filtering and find models which will fit well. By calibrating the velocity models at long periods, with a ground truth event, will help in trouble shooting problems. It can be real difficult to trouble shoot with smaller events with more noiser data in seperating source, instrument or structure effects?

Now run the C-shell script by making it executable chmod a+x runall.csh or typing
```
csh runall.csh
```

Your directory ./Event/01/ should look like this:

1999.213.16.06.04.0000.BK.CMB..BHE.SAC      anmo.data                                   new.gbin
1999.213.16.06.04.0000.BK.CMB..BHN.SAC      anmo.gbin                                   new.gbin.bak
1999.213.16.06.04.0000.BK.CMB..BHZ.SAC      anmo.gbin.bak                               new.par
1999.213.16.06.08.9940.II.PFO..BHE.SAC      anmo.par                                    new.r.dis
1999.213.16.06.08.9940.II.PFO..BHN.SAC      anmo.r.dis                                  new.t.dis
1999.213.16.06.08.9940.II.PFO..BHZ.SAC      anmo.t.dis                                  new.z.dis
1999.213.16.06.13.0000.TS.PAS..BHE.SAC      anmo.z.dis                                  pas.data
1999.213.16.06.13.0000.TS.PAS..BHN.SAC      cmb.data                                    pas.gbin
1999.213.16.06.13.0000.TS.PAS..BHZ.SAC      cmb.gbin                                    pas.gbin.bak
1999.213.16.06.26.9880.US.DUG..BHE.SAC      cmb.gbin.bak                                pas.par
1999.213.16.06.26.9880.US.DUG..BHN.SAC      cmb.par                                     pas.r.dis
1999.213.16.06.26.9880.US.DUG..BHZ.SAC      cmb.r.dis                                   pas.t.dis
1999.213.16.07.06.0106.IU.TUC.00.BHE.SAC    cmb.t.dis                                   pas.z.dis
1999.213.16.07.06.0106.IU.TUC.00.BHN.SAC    cmb.z.dis                                   runall.sh
1999.213.16.07.06.0106.IU.TUC.00.BHZ.SAC    dug.data                                    tuc.data
1999.213.16.07.28.9858.IU.ANMO.00.BH1.SAC   dug.gbin                                    tuc.gbin
1999.213.16.07.28.9858.IU.ANMO.00.BH2.SAC   dug.gbin.bak                                tuc.gbin.bak
1999.213.16.07.28.9858.IU.ANMO.00.BHZ.SAC   dug.par                                     tuc.par
1999.213.16.07.55.0070.US.NEW..BHE.SAC      dug.r.dis                                   tuc.r.dis
1999.213.16.07.55.0070.US.NEW..BHN.SAC      dug.t.dis                                   tuc.t.dis
1999.213.16.07.55.0070.US.NEW..BHZ.SAC      dug.z.dis                                   tuc.z.dis
WLBR_SPYD_990801_1606.SAC.tar               new.data

The files that end in *.data are the 3 component data in binary format. Files ending in *.gbin are the 10 Greens functions in binary format. The files ending *.par are text files which serve as headers for *.data and *.gbin binary files. Files ending in *.[r,t,z].dis are SAC files which have been instrument corrected into ground displacement (cm).

Make a file called event.in using a text editor.
```
vi event.in
```

Type this into the file event.in

Scottys Junction Earthquake 8/1/1999
5
1000 0.1 1.0E-5
0 0 0    5.5
10 10 10  0.1
anmo.gbin anmo.data anmo.par 10 1.0 0.0 0.0
dug.gbin  dug.data  dug.par  10 1.0 0.0 0.0
new.gbin  new.data  new.par  10 1.0 0.0 0.0
pas.gbin  pas.data  pas.par  10 1.0 0.0 0.0
tuc.gbin  tuc.data  tuc.par  10 1.0 0.0 0.0

See the man page GridSearch4 or sawn for the input format of the file event.in. The input format is basically the same for both programs.

Do a grid search. This will do about 3000 iterations over a course grid and then 9000 iterations over a finer grid given the lowest misfit from the first course grid search. Depending on the speed of you computer and number of statons in the search, this may take a few minutes to a half hour.
```
GridSearch4 par=event.in > gs.out
```
See the man page GridSearch4 for command line options.

See the results of misfit versus model parameters. For this case you should of obtained a solution of

nodal plane 1: azm=84 dip=74 rak=-51 Mw=5.45 resid=7.62855e-05
nodal plane 2:    193     42    -156

plotgrid if=gs.out range=0.5 logplot=1

The plot should look something like this grid.gif

A successful regional moment tensor inversion depends on well calibrated greens functions (velocity models). In regions of sparse covreage, velocity models will never be known well enough at higher frequencies and 1D models may not reflect realistic 3D earth. My philosphy is that a iterative approach is more suitable than a maxtrix inversion. A large number of regional observations are made and misfits are averaged over each observation.

See the man page plotgrid for the command line arguments.

To use the waveform optimization program, run sawm. Depending on the speed of you computer, this may take a few minutes. In this case it is only 1000 iterations.
```
sawm -iff par=event.in > sawm.out
```
See the man page sawm for the command line arguments.

To see the waveform data and synthetic fits, use plplotsyn and type the commands with the solutions obtained from GridSearch4 or sawm. Say you got these values from sawm:

plplotsyn par=event.in azm=184.332 dip=41.8908 rak=-158.277 Mw=5.5

And these values from GridSearch4:

plplotsyn par=event.in azm=84 dip=74 rak=-51 Mw=5.45

The plot should look something like this waveform.gif

Your results may vary and the best solution is obtained by rerunning sawm once or twice more to make sure you get back to this solution. If not then you may need to run the program GridSearch4. As you can see for this case plplotsyn shows that the synthetics fit the data pretty well and that the sawm solution is only 10 degrees off in strike but nailed the other parameters. The gridsearch results show that the best solution at this depth still has some misfit at station TUC and DUG. You can use the script FILTER_SAC to further filter the data and synthetics for a better fit. The band between 10 and 100 seconds will be more independent of the velocity structure.

Type the commands at the prompt.

FILTER_SAC tuc 0.01 0.02
FILTER_SAC dug 0.01 0.02

Then retry GridSearch4 or sawm.

See the man page plplotsyn for the command line arguments.

You can rerun PROCESS_MACRO using different depths and then GridSearch4 or sawm to find the optimal depth. Heres one way you can do it in a C-shell script.

#!/bin/csh

cat > depth.in << EOF
Scottys Junction Earthquake 8/1/1999
5
1000 0.1 1.0E-5
0 0 0    5.5
10 10 10  0.1
anmo.gbin anmo.data anmo.par 10 1.0 0.0 0.0
dug.gbin  dug.data  dug.par  10 1.0 0.0 0.0
new.gbin  new.data  new.par  10 1.0 0.0 0.0
pas.gbin  pas.data  pas.par  10 1.0 0.0 0.0
tuc.gbin  tuc.data  tuc.par  10 1.0 0.0 0.0
EOF

PROCESS_MACRO tuc  0.01 0.05 256 1.75 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO dug  0.01 0.05 256 1.00 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO cmb  0.01 0.05 256 0.80 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO pas  0.01 0.05 256 1.50 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO anmo 0.01 0.05 256 1.75 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response
PROCESS_MACRO new  0.01 0.05 256 2.00 wus 37.39 -117.07 5.0 1999 213 16 06 22 ../../stalist ../../Response

foreach depth ( 02 04 06 08 10 12 14 16 18 20 )

mpf sta=tuc  dt=1.75 nt=256 fmax=0.02 mod=wus ela=37.39 elo=-117.07 edp=$depth db=../../stalist
mpf sta=dug  dt=1.00 nt=256 fmax=0.02 mod=wus ela=37.39 elo=-117.07 edp=$depth db=../../stalist
mpf sta=pas  dt=1.50 nt=256 fmax=0.05 mod=wus ela=37.39 elo=-117.07 edp=$depth db=../../stalist
mpf sta=anmo dt=1.75 nt=256 fmax=0.05 mod=wus ela=37.39 elo=-117.07 edp=$depth db=../../stalist
mpf sta=new  dt=2.00 nt=256 fmax=0.05 mod=wus ela=37.39 elo=-117.07 edp=$depth db=../../stalist

MAKE_SYN tuc  0.01 0.02
MAKE_SYN dug  0.01 0.02
MAKE_SYN pas  0.01 0.05
MAKE_SYN anmo 0.01 0.05
MAKE_SYN new  0.01 0.05

GridSearch4 par=depth.in > depth${depth}.out
end

mtinvl5

I have also written a moment tensor inversion program but it still needs some work on the problem of finding the best correlation between synthetics and data. The previous programs are iterative so that at each step we use a trial solution with a shift corresponding to the maximum waveform correlation. In a matrix inversion, the shift is not known and becomes a nonlinear model parameter which has to be inverted along with the moment tensor and optimal source depth. Some have removed this shift value by simultaneously solving for the shift and moment tensor by inverting for a Centroid Moment Tensor (CMT-Harvard Group) where the location and origin time become part of the inversion. The reason why we cannot use this procedure is that we do not know the velocity model well enough on the regional scale to use the CMT technique. This makes the iterative solution a more suitable but cannot compete with the efficiency of the matrix inversion.

The moment tensor inversion code can be run by adding these lines to runall.sh. See example in directory ./LeeVining/19901024061519. Notice that rather than PROCESS_MACRO we run PROCESS_MACRO2 which runs makegreen6 rather than makegreen5. The difference is in the format of the Greens function output as ZSS,ZDS,ZDD,ZEP,RSS,....,TDS.

PROCESS_MACRO2 pas  0.02 0.05 256 0.9 br1 38.0523338 -119.126663 12 1990 297 06 15 19 ../../stalist ../../Response
PROCESS_MACRO2 gsc  0.02 0.05 256 0.9 wus 38.0523338 -119.126663 12 1990 297 06 15 19 ../../stalist ../../Response

The data is processed and the Greens functions should now be computed. Now to do the inversion make a file name something like mtinv.in.

1990 OCt 24 Lee Vining 
2
0 0  0 0 0 0 
gsc.gbin gsc.data gsc.par 1 1.0 0.0 0.0 
pas.gbin pas.data pas.par 1 1.0 0.0 0.0

First line is a comment line, second is the number of stations, third line is the 6 moment tensor elements Mxx,Myy,Mzz,Mxy,Mxz,Myz. The rest of the file is very similar in format to the sawm.in input files except that the lagtol field now specifies a manual shift in the data rather than an automatic shift for the alignment bewteen synthetics and data. Run the program by typing this at the command prompt:

mtinvl5 par=mtinv.in

Play around with the lagtol field to manually shift the seismograms to get the maximum cross correlation. Also play around with different depths. I get an optimal depth of 12 km. You should get an answer like...

Mxx,Myy,Mzz,Mxy,Mxz,Myz
-6.93651e+07 7.0955e+07 48310.1 -2.8499e+06 1.01042e+07 6.43817e+06 The deviatoric moment tensor due to a major double couple: T-axis trend= 91 plunge= 5 P-axis trend= 182 plunge= 8 NP1: stk= 226 dip= 81 rak= -2 NP2: stk= 316 dip= 88 rak= -171

Run mtfwd to visuallize the results. Edit the file mtinv.in by adding the solution to the third line.

1990 OCt 24 Lee Vining 
2
-6.93651e+07 7.0955e+07 48310.1 -2.8499e+06 1.01042e+07 6.43817e+06
gsc.gbin gsc.data gsc.par 1 1.0 0.0 0.0 
pas.gbin pas.data pas.par 1 1.0 0.0 0.0

mtfwd par=mtinv.in

The output should look something like this.

Manual Pages Index

MAKE_SYN -- C-shell script that makes and filters synthetics.
FILTER_SAC -- C-shell script that re-filters both data and synthetics.
GridSearch4 -- Perform a grid search for minimum misfit of data and synthetics over fault strike, dip, rake, and seismic moment.
PROCESS_MACRO -- C-shell script that filters and processes seismograms.
gbin2sac -- Converts binary data to SAC files.
makedata -- Converts SAC files to binary data.
makedisp -- Given src-rec coordinates, source depth, fault strike, dip, rake, and seismic moment, makedisp will compute 3 compoent (North,East,Vertical) displacement synthetics in SAC format.
makegreen5 -- Computes Green's functions.
makegreen6 -- Computes Green's functions and outputs in fundamental faulting orientations for use in moment tensor inversion.
mpf -- Makes text *.par files that carries header information for each station.
myjulday -- returns Julian day
plotgrid -- Plots mean amplitude misfit of data and synthetics versus fault strike, dip, rake, and seismic moment. Takes input from prgrams GridSearch4 and sawm.
plplotsyn -- Plots data and synthetics given a fault strike, dip, rake, and seismic moment.
sac2gbin -- Converts SAC files to binary data.
sacsun2linux -- Byte swapping.
sawm -- Perform a simulated annealing - downhill Simplex iterative waveform inversion to determine the fault strike, dip, rake, and seismic moment.
mtinvl5 -- Moment Tensor Inversion
mtfwd -- Display Moment Tensor Inversion Results.

MAKE_SYN

NAME: Make Synthetics (MAKE_SYN)
VERSION 1.0
AUTHOR:
LANGUAGE: C-shell Script
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO: makegreen5, gbin2sac, SAC, sac2gbin
PROGRAM EXAMPLE:
```
MAKE_SYN kcc 0.01 0.1 
```
MAKE_SYN will execute makegreen5 kcc.par and output kcc.gbin (Green's functions) It will then execute gbin2sac kcc.par to demultiplex the Green's functions into 10 components in SAC files kcc.g01, kcc.g02, ..., kcc.g10. MAKE_SYN will then execute SAC and bandpass the files from 0.01 to 0.1 Hz using the same parameters as with the data. sac2gbin kcc.par is then executed overwrite the kcc.gbin file. The original unfiltered Green's functions are saved in file kcc.gbin.bak so that FILTER_SYN can be executed to refilter the synthetics if needed in which case the Green's functions will not have to be recalculated.

PROCESS_MACRO

NAME: Data Processing Macro
VERSION 1.0
AUTHOR:
LANGUAGE: C-shell Script
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO: mpf, makedata, SAC
PROGRAM EXAMPLE:
```
PROCESS_MACRO kcc  0.01 0.10 256 0.30 wus 39.1234 -118.12345 1999 365 23 59 59 ../../stalist ../../Response
```
This example of PROCESS_MACRO processes SAC files from station kcc. The entered fields are used to execute program mpf which returns important values from a look up table of station locations in a formated file named ../../stalist and velocity models embedded in the program itself. Additional stations can be entered in file stalist using a text editor. SAC is then executed to read the SAC files assuming the RDSEED file nameing convention. The values 39.1234 -118.12345 is the earthquake location. The values 1999 365 23 59 59 are the origin time of the earthquake and get set in the SAC header. The zero time in all the files are set to this origin time. The code wus stands for the Western US velocity model determine by Ritsema and Lay (1995). There are several velocity models to choose from and the *.par files for each station can be edited or created using a text editor to input a model of your choice. Commonly used models can be input to the source code of program mpf.c. All of this information is stored in a header file kcc.par which can be view and edited with a text editor.
The instrument correction is applied from Poles and Zeros files and can be located in any directory like specified here as ../../Response. The files must have names like KCC_ZVBB.zp KCC_NVBB.zp KCC_EVBB.zp. The C-shell Script can be edited to adapt to different *.zp file names. The files are rotated, itergrated to displacement, converted to cm units of ground motion, interpolated so that it will have the same number of points npts=256 and length ( npts=256 ) x ( dt=0.30 seconds ) = 76.8 seconds as the corresponding synthetics and finally bandpass filter from 0.01 to 0.1 Hz. Original displacement files are saved as files kcc.z.dis kcc.r.dis kcc.t.dis so that the script FILTER_SAC can be used to refilter the files without having to reprocess the data. makedata kcc.par is finally executed to convert displacement files into one kcc.data file.

FILTER_SAC

NAME: Filter binary files through SAC (FILTER_SAC)
VERSION 1.0
AUTHOR:
LANGUAGE: C-shell Script
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO: makedata, gbin2sac, sac2gbin
PROGRAM EXAMPLE:
```
FILTER_SAC kcc 0.02 0.05
```
This example uses FILTER_SAC to read files kcc.gbin.bak created by MAKE_SYN and runs gbin2sac kcc.par to convert multiplex gbin file to demultiplex sac files and then runs SAC to bandpass the synthetics from 0.02 to 0.05 Hz. The files are then converted back to kcc.gbin format using the command sac2gbin kcc.par. The data is then filtered by SAC using the kcc.z.dis kcc.r.dis kcc.t.dis files saved by PROCESS_MACRO command. The command makedata kcc.par is finnally run to convert the displacement SAC files into kcc.data file.

GridSearch4

NAME: Grid Search
VERSION 4.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO: plotgrid

PROGRAM EXAMPLE:

Create the input file by using a text editor or by cloning and editing a previous version. The format of file gridsearch.in:

{event comment}
{number of stations}
{number of iterations}  {temperature}   {tolerance}
{initial strike}        {initial dip}   {initial rake}   {initial Mag}
{increment strike}      {increment dip} {increment rake} {increment Mag}
{bin file}              {data file}     {parameter file} {lag tolerance}  {% weight} {rise time} {pulse width}
2nd station
3rd station
.
.
.
{bin file}              {data file}     {parameter file} {lag tolerance}  {% weight} {rise time} {pulse width}

event comment - optional text string.
number of stations - for now you can gridsearch up to 11 stations.
number of iterations - Ignore, not used in this program
temperature - Ignore, not used in this program
tolerance - Ignore, not used in this program
initial strike - Ignore, not used in this program
initial dip - Ignore, not used in this program
initial rake - Ignore, not used in this program
initial Magnitude - set to minimum magnitude possible for this event
increment strike - Ignore, not used in this program
increment dip - Ignore, not used in this program
increment rake - Ignore, not used in this program
increment Magnitude - set search increment for the range of magnitudes
bin file - name of Green's function file e.g. kcc.gbin
data file - name of data file e.g. kcc.data
parameter file - name of parameter file e.g. kcc.par
lag tolerance - the maximum allowed number of points to shift the synthetic for maximum correlation with the synthetics. A value of zero means that the synthetics will not be shifted.
% weight - the percent factor to use for this station in determining the average residual over all of the stations.
rise time - in seconds. This is the rise time for a trapazoid source time function which is convolved with the synthetics.
pulse width - in seconds. The width of the top flat portion of the trapazoid source time function.

Example of file gridsearch.in:

02 NOV 1997 08:51 UTC Fish Lake Valley, Nevada 4 500 0 1.0E-12 70 70 10 4.9 15 1 15 0.1 kcc.gbin kcc.data kcc.par 10 1.0 0.0 0.0 tuc.gbin tuc.data tuc.par 10 1.0 0.0 0.0 cor.gbin cor.data cor.par 10 1.0 0.0 0.0 anmo.gbin anmo.data anmo.par 10 1.0 0.0 0.0

Execute GridSearch4 at the command prompt

GridSearch4 par=gridsearch.in > gridsearch.out

This example will grid search for the best strike dip rake and magnitude over four stations kcc, tuc, cor, and anmo. An initial rough search is made and then a second refined search is made over a smaller range equivalent to the increments. Strike dip and rake are usually resolveable to within 3 degrees. and moment magnitude to within 0.1 units. An inital magnitude of 4.9 is used and will search over 5 magnitudes 4.9, 5.0, 5.1, 5.2, 5.3 with increment of 0.1. Chose a wider magnitude range if needed. See plotgrid to see the results.

gbin2sac

NAME: Gbin file to SAC (gbin2sac)
VERSION 1.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION:
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE:
You can execute gbin2sac at the command prompt by
```
gbin2sac kcc.par
```
but this program usually executes in a C-shell scripts MAKE_SYN and FILTER_SAC.
Reads the binary file kcc.gbin with the 10 Green's functions and writes ten SAC files named kcc.g01 kcc.g02 kcc.g03 kcc.g04 kcc.g05 kcc.g06 kcc.g07 kcc.g08 kcc.g09 kcc.g10.

makedata

NAME: Make Data files (makedata)
VERSION
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE: You can execute makedata at the command prompt by
```
makedata kcc.par
```
but this program usually executes in the C-shell script PROCESS_MACRO.
Reads three SAC files defined in file kcc.par (radial,transverse,vertical components) and makes the binary file kcc.data.

makegreen5

NAME: makegreen5
VERSION 5
AUTHOR: Yuehua Zeng
LANGUAGE: Fortran 77
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:

PROGRAM EXAMPLE:

You can execute the program makegreen5 at the command prompt:

makegreen5 kcc.par

The program will compute the 10 Green's functions in output file kcc.gbin. This program is run in C-shell script MAKE_SYN. The input format is text and fairly straight forward. Parameters can be edited using a text editor. The input can be in any order and fields and input are seperated by and equal sign. The input parameter list:

sacfile - name of which files *.gbin *.data *.par will be identified. I usually use the station id name.
dfile - set same as sacfile name
radfile - name of SAC file for the raidal component.
verfile - vertical component.
trafile - transverse component.
rlat - latitude of receiver (decimal degrees)
rlon - longitude of receiver
rdep - elevation of receiver
slat - latitude of source
slon - longitude of source
sdep - depth of source
dt - sampling rate (not used but is set by npts and twin)
t0 - start time of synthetics reference to origin time
npts - number of points
twin - time lenght of synthetics in seconds
fmax - maximum frequency to compute (usually set to 0.25 Hz above the band of interest)
damp - controls wrap around problem created by FFT. Not needed unless twin is very small... basically slow waves arriving later in the record are wrapped around and apparently arrive early and contaminat the body waves. (set to 1 for no damping, set > 1 for damping)
al - multiple sources (Set to 0 for now)
kmax - maximum wavenumbers. Need to increase when distances become very large (>10,000km)
smin - Noise tolerance for source? set smin to 0.0001, there is no reason to change.
eps - tolerance? set eps to 0.0001, there is no reason to change.
# - comments
nlay - number of layers in velocity model
lay01 - first layer
lay02 - second deepest layer, thickness=0 halfspace
lay{2-digit layer number}={thickness},{p-wave velocity km/s},{Qp},{S-wave velocity km/s},{Qs},{rho g/cc}

Example file kcc.par

sacfile=kcc dfile=kcc radfile=kcc.r.bp verfile=kcc.z.bp trafile=kcc.t.bp rlat=37.32 rlon=-119.3 rdep=0.0 slat=37.8463 slon=-118.215 sdep=10 dt=0.3 t0=0 npts=256 twin=76.8 fmax=0.1 damp=1.0 al=0 kmax=10000 smin=0.0001 eps=0.0001 ######### Western US ritsema and lay 1995 JGR ############### nlay=4 lay01=4.00,4.52,500.0,2.61,250.0,2.39 lay02=28.0,6.21,500.0,3.59,250.0,2.76 lay03=20.0,7.73,1000.,4.34,500.0,3.22 lay04=0.00,7.64,1000.,4.29,500.0,3.19

makegreen6

NAME: makegreen6
VERSION 6
AUTHOR: Yuehua Zeng
LANGUAGE: Fortran 77
LIBRARIES REQUIRED: none
SYNOPSIS:
DESCRIPTION:
NOTE:
BUGS
REFERENCES: Jost, M. L. and R. B. Herrmann, (1989). A student's guide to and review of moment tensors, Seismological Research Letters, Vol 60, 37-57.
SEE ALSO:

PROGRAM EXAMPLE:

makegreen6 kcc.par

This program is the same as makegreen5 but outputs the Green's functions in the fundamental faulting orientations...

kcc.g01 - RSS - radial component vertical strike-slip
kcc.g02 - RDS - radial component dipping strike-slip
kcc.g03 - RDD - radial component dipping dip-slip
kcc.g04 - REP - radial component explosion
kcc.g05 - ZSS - vertical component vertical fault strike-slip
kcc.g06 - ZDS - vertical component dipping strike-slip
kcc.g07 - ZDD - vertical component  dipping dip-slip
kcc.g08 - ZEP - vertical component explosion
kcc.g09 - TSS - transverse component vertical fault strike-slip
kcc.g10 - TDS - transverse component dipping dip-slip

This is for use with the moment tensor inversion program mtinvl2

mpf

NAME: Make Parameter File (mpf)
VERSION 0.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED: none
SYNOPSIS:
DESCRIPTION:
NOTE:
BUGS
REFERENCES:
SEE ALSO:

PROGRAM EXAMPLE: The C-shell script PROCESS_MACRO uses this program to look up values and make a parameter file.

mpf sta=kcc dt=0.3 nt=256 mod=wus fmax=0.1 vmin=2 vmax=10 t0=0 ela=34.00 elo=-118.0 edp=10 db=../../stalist

These values are returned to stdout

kcc KCC 34 -118 10 37.32 -119.3 0 256 0.3 76.5

kcc  - station name small caps
KCC  - station name capitals
34   - event lat
-118 -  event long
10    - event depth
37.32  - station lat
-119.3  - station lon
0       - station elev
256     - nt
0.3     - dt
76.5    - twin - dt (minus dt so that SAC will cut exactly 256 points in PROCESS_MACRO)

A file kcc.par will be written

sacfile=kcc 
dfile=kcc 
radfile=kcc.r.bp
verfile=kcc.z.bp
trafile=kcc.t.bp
rlat=37.32
rlon=-119.3
rdep=0.0
slat=34
slon=-118
sdep=10
dt=0.3
t0=0
npts=256
twin=76.8
fmax=0.1
damp=1.0
al=0
kmax=10000
smin=0.0001
eps=0.0001
######### Western US ritsema and lay 1995 JGR  ###############
nlay=4
lay01=4.00,4.52,500.0,2.61,250.0,2.39
lay02=28.0,6.21,500.0,3.59,250.0,2.76
lay03=20.0,7.73,1000.,4.34,500.0,3.22
lay04=0.00,7.64,1000.,4.29,500.0,3.19

The file ../../stalist has the format

kcc 37.32000 -119.30000 Kings Canyon cmb 38.03500 -120.38300 Columbia sao 36.76500 -121.44500 San Andreas Observatory

plotgrid

NAME: Plot Grid (plotgrid)
VERSION 1.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED: PgPlot
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE:
To execute plotgrid at the command prompt:
```
plotgrid if=gridsearch.out range=0.5 logplot=1 gif=0 ps=0
```
This example executes plotgrid which reads input from file gridsearch.out from output by GridSearch4. It can also plot residuals output by program sawm. range=0.5 only plots 50% of the residuals to zoom in on region of the lowest residuals. The default is range=1.0 100%. logplot=1 will plot the residual on a log axis. The default is logplot=0, no log axis. If gif=1 or ps=1 then plotgrid will make a GIF or Postscript graphics plot for printing. The default is gif=0 or ps=0 and output goes to a X-windows screen. The graphics output looks like this

myjulday

NAME: My Julian Day (myjulday)
VERSION
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE: This program can be used to convert year month day from earthquake catalogs to jullian day automatically in shell scripts.
```
myjulday 1 1 2000
```
returns the jullian day without a carrage return
```
001
```
For executing in a C-shell script:
#!/bin/csh -x set year=2000 set month=01 set day=01 set jday=`myjulday $day $month $year` echo $jday

plplotsyn

NAME: Plot Synthetics (plplotsyn)
VERSION 1.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED: PgPlot
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE:
To plot the stations in file sawm.in execute on the command promp:
```
plplotsyn par=sawm.in azm=342 dip=78 rak=-174 Mw=5.15 gif=0 ps=0
```
The file sawm.in has the same format as gridsearch.in and can be used for both programs sawm and GridSearch4. If gif=1 or ps=1 then plotgrid will make a GIF or Postscript graphics plot for printing. The default is gif=0 or ps=0 and output goes to a X-windows screen.
The output looks like:

sac2gbin

NAME: SAC to Gbin (sac2gbin)
VERSION 1.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE:
You can execute sac2gbin at the command prompt by
```
sac2gbin kcc.par
```
but this program usually executes in a C-shell scripts MAKE_SYN and FILTER_SAC.
Reads the ten SAC files of the Green's functions named kcc.g01 kcc.g02 kcc.g03 kcc.g04 kcc.g05 kcc.g06 kcc.g07 kcc.g08 kcc.g09 kcc.g10 and writes a binary file kcc.gbin.

sacsun2linux

NAME: SAC SUNOS to LINUX sacsun2linux
VERSION 1.0
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:
PROGRAM EXAMPLE: To convert by byte-swapping the binary SAC file kcc.z.dis from "Motorola Processor-SunOS, SolarisOS" Big Endian byte format to "Intel Processor-LinuxOS, MSWindowsOS, SolarisOS" Little Endian byte format type at the command prompt:
```
sacsun2linux kcc.z.dis
```
sacsun2linux will also accept wild cards. All of the header information is also transfered. Conversely:
```
saclinux2sun kcc.z.dis
```
will convert back to Big Endian byte.

sawm

NAME: Simulated Annealing Waveform Modeling (sawm)
VERSION 1.0
AUTHOR:
LANGUAGE:
LIBRARIES REQUIRED:
SYNOPSIS:
DESCRIPTION: Uses simulated annealing and downhill simplex optimization.
NOTES:
BUGS
REFERENCES:
Press, W H, S A Teukolsky, W T Vetterling and B P Flannery, Numerical Recipies in C, 2nd Ed, Cambridge Press, 1992.
Nelder, J A and R Mead, 1965, Computer Journal, vol 7, pp. 308-313.
SEE ALSO: plotgrid, plplotsyn

PROGRAM EXAMPLE:

To execute sawm at the command prompt:

sawm -iff par=sawm.in > sawm.out

The -iff switch makes sawm use a predefined set of intial fault orientations. The user can also specify an intial fault orientations in the input file. The -norm switch will normalize all the amplitudes from -1 to +1

sawm -norm par=sawm.in > sawm.out

inital magnitude and magnitude increment still needs to be defined in the input file.

Create the input file by using a text editor or by cloning and editing a previous version. The format of file sawm.in is the same as that used by plplotsyn and GridSearch4 :

{event comment}
{number of stations}
{number of iterations}  {temperature}   {tolerance}
{initial strike}        {initial dip}   {initial rake}   {initial Mag}
{increment strike}      {increment dip} {increment rake} {increment Mag}
{bin file}              {data file}     {parameter file} {lag tolerance}  {% weight} {rise time} {pulse width}
2nd station
3rd station
.
.
.
{bin file}              {data file}     {parameter file} {lag tolerance}  {% weight} {rise time} {pulse width}

event comment - optional text string.
number of stations - for now you can use up to 11 stations.
number of iterations - maximum number of iterations, somewhere between 500 100 iterations should work.
temperature - This control parameter is like temperature in simulated annealing. A larger value allows random fluctuations to be accepted more often to downhill optimzation and avoids local minima in hopes to find the global minima. A value of zero reduces the algorithm exactly to the downhill simplex method and converges to the closest local minima. A number of 0.1 for simulated annealing works well.
tolerance - In the downhill simplex method, the iterations stop when the residual=tolerance is reached. The program may also stop if the number of maximum iterations is reached or if the last residual equals the next residual. A value of 1.0E-5 works well.
initial strike - set to whatever you like, not really used
initial dip - set to whatever you like, not really used
initial rake - set to whatever you like, not really used
initial Magnitude - Use something close to the catalog magnitude.
increment strike - set to whatever you like, not really used
increment dip - set to whatever you like, not really used
increment rake - set to whatever you like, not really used
increment Magnitude - Use small value, like 0.001 when temperature=0 or 0.0 when temperature not 0
bin file - name of Green's function file e.g. kcc.gbin
data file - name of data file e.g. kcc.data
parameter file - name of parameter file e.g. kcc.par
lag tolerance - the maximum allowed number of points to shift the synthetic for maximum correlation with the synthetics. A value of zero means that the synthetics will not be shifted. A value of 10 points works well
% weight - the percent factor to use for this station in determining the average residual over all of the stations. Use a value of 1.0 for 100%.
rise time - in seconds. This is the rise time for a trapazoid source time function which is convolved with the synthetics.
pulse width - in seconds. The width of the top flat portion of the trapazoid source time function.

Example of file swam.in:

02 NOV 1997 08:51 UTC Fish Lake Valley, Nevada 4 500 0 1.0E-12 70 70 10 5.1 15 1 15 0.01 kcc.gbin kcc.data kcc.par 10 1.0 0.0 0.0 tuc.gbin tuc.data tuc.par 10 1.0 0.0 0.0 cor.gbin cor.data cor.par 10 1.0 0.0 0.0 anmo.gbin anmo.data anmo.par 10 1.0 0.0 0.0

mtinvl2

NAME: Moment Tensor Inversion L2 norm
VERSION 5
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
```
mtinvl5 par= 
```
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:

PROGRAM EXAMPLE:

The input file format is similar to sawm and GridSearch4 input format.

comment
number of statons
Mxx Myy Mzz Mxy Mxz Myz
data gbin par shift rise_time pulse_width
.
.
.
data gbin par shift rise_time pulse_width

For Example

1990 OCt 24 Lee Vining
2
0 0  0 0 0 0
gsc.gbin gsc.data gsc.par 1 1.0 0.0 0.0
pas.gbin pas.data pas.par 1 1.0 0.0 0.0

mtfwd

NAME: Moment Tensor Forwad Model
VERSION 1
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:
SYNOPSIS:
```
mtfwd par=
```
DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO:

PROGRAM EXAMPLE:

The input file format is similar to sawm and GridSearch4 input format.

comment
number of statons
Mxx Myy Mzz Mxy Mxz Myz
data gbin par shift rise_time pulse_width
.
.
.
data gbin par shift rise_time pulse_width

For Example

1990 OCt 24 Lee Vining
2
-6.93651e+07 7.0955e+07 48310.1 -2.8499e+06 1.01042e+07 6.43817e+06 
gsc.gbin gsc.data gsc.par 1 1.0 0.0 0.0
pas.gbin pas.data pas.par 1 1.0 0.0 0.0

Copy the output from the terminal obtained from running mtinvl5 and add it to the third line in the above parameter file.

makedisp

NAME: Make Displacement Seismograms
VERSION 1
AUTHOR:
LANGUAGE: C
LIBRARIES REQUIRED:

SYNOPSIS:

makedisp par= azm= dip= rak= Mw= yr= mo= da= hr= min= sec= msec= tz= tt= tr=

DESCRIPTION
NOTE:
BUGS
REFERENCES:
SEE ALSO: mpf MAKE_SYN

PROGRAM EXAMPLE:

This example will compute the greens functions, bandpass filter and compute 3 component seismograms in (R,T,Z) reference frame.


mpf sta=kcc dt=0.03 nt=512 t0=0 fmax=1.0 mod=sn ela=34.0 elo=-118.0 edp=10 db=stalist

MAKE_SYN kcc 0.01 0.5

makedisp par=kcc.par azm=360 dip=90 rak=-180 Mw=4.0 yr=2000 mo=1 da=1 hr=1 min=1 sec=1 msec=1 tz=0 tt=0 tr=0

Command line arguments:

par=kcc.par  - parameter file
azm=360 - fault strike degrees from 0 north
dip=90  - fault dip from horizontal
rak=-180  - fault rake (-180 or 180 degrees) right lateral, (0 degrees), (pure normal -90 degrees), (pure reverse 90)
Mw=4.0   - Moment Magnitude
Origin Time
yr=2000  - 4 digit year
mo=1  - 2 digit month
da=1  - 2 digit day
hr=1  - 2 digit hour
min=1 - 2 digit minutes
sec=1 - 2 digit seconds
msec=1 - 4 digit milliseconds
tz=0  - time zone in hours from UTC,GMT tz=0 
tr=0 - rise time in seconds of trapazoid source time function (default: point source tr=0, tt=0)
tt=0 - pulse width in seconds of flat top portion of trapazoid source time function