SPHYPIT90 MANUAL

Table of Contents

• AUTHOR
• INTRODUCTION AND OVERVIEW
• THEORETICAL BACKGROUND
• I/O
• INPUT FILES
• ARRIVAL TIME FILES
• THE STATION LIST FILE
• THE 1D VELOCITY MODEL FILE
• THE 3D VELOCITY MODEL FILE
• STATIC INPUT FILES
• THE PARAMETER FILE
• DAMPER AND INVERSE COVARIANCE MATRIX FILE
• OUTPUT FILES
• SPHREL3D NOTES
• DIMENSIONS
• SUBROUTINES
• KNOWN AND SUSPECTED BUGS
• ACKNOWLEDGEMENTS

AUTHOR

INTRODUCTION AND OVERVIEW

The tomographic method used is based on that of Aki et al. (1976) but with many extentions and modifications. These modifications have been documented in various publications, mostly in JGR (I'll try to get a bibliography together sometime) - most are reviewed in this document. The inversion procedure is based on that presented by Tarantola and Valette (1982).

The purpose of this document is to provide the potential user of SPHYPIT90 some useful information about the program. It is by no means complete, nor is it likely to win any awards for clarity or conciseness. It is, however, better than nothing, which I'm afraid is the only other option. Anyone who uses this program is encouraged to not treat it as a black box, but to look at and understand the code before committing a lot of time and energy (and funds) towards running it.

SPHYPIT90 works in conjunction with other programs, most importantly SPHREL3D90 which relocates earthquakes in the 3D models produced by SPHYPIT90. Many of the conventions and subroutines are the same for both these programs, so there will be some overlap in their description. In particular, both programs use the same I/O routines even though a number of parameters and files are exclusive to one or the other routine. There is at present no separate manual for SPHREL3D90 but most of what you need to know is discussed here.

SPHYPIT90 also requires a couple of auxilliary routines in order to perform the inversion of complicated covariance matrices. I'll write that up later.

Note that 2001.0404 is a beta version of this routine, meaning that it has undergone only rudimentary testing. The procedures used in this program have a long heritage, however, so the bugs are likely to be simple ones (and easily fixed). Let me know if you run into trouble.

Some History

This program grew out of a travel time inversion project started in 1977, and, as I am a real softy at giving out software, there are many different versions of it floating around. The two most important varieties are the flat earth program (called HYPIT or HYPIT89) and the extension of this program to a spheroidal earth (spherical HYPIT or SPHYPIT). For all practical purposes the spheroidal program supersedes the flat earth program, so all future development will be directed towards it.

Much of the original development was done on a SUN 3/160, although all programing has been done in simple, non-machine specific FORTRAN that allows it to be easily portable to different machines.

THEORETICAL BACKGROUND

The earth is parameterized by a mesh of constant velocity blocks, each with its own P and S velocity. All boundaries between blocks are horizontal or vertical in spherical coordinates, and an earth-flattening transformation is applied to reduce the problem to a flat earth system.

The algorithm follows the formalism of Tarantola and Valette (1982), and solves for the k+1^th iterate of the problem

where d is the arrival time data, p is the model consisting of the hypocenter coordinates and the velocity model (and occasionally station corrections) and g is the functional relation between the two (in this case, the times provided by the ray tracer). We determine p by applying the following formula:

where

po	The initial a priori estimate of the model (p).
pk	The estimate of p at the kth iteration
Gk	The m x n matrix of partial derivatives (dg/dp)
Cdd	The a priori data covariance matrix
Cpp	The a priori model covariance matrix
do	This is essentially the same as d (the data).

At every iteration SPHYPIT90 solves for the perturbation p_k+1 - p_k. In early versions of SPHYPIT90, p_k was fractional slowness but now is slowness (1/velocity).

Calculation of G and g(p)

Local Arrivals

The general equation for describing the arrival time for a ray from earthquake i to station j (T_ij) in the linear hypocenter-velocity problem is

or

where h refers to the hypocenter coordinates (origin time, latitude, longitude, and depth) and s can be either the slowness or the fractional slowness change (this version uses slowness exclusively).

Ray tracing is always done in a one-dimensional structure. If a three dimensional structure is used as a starting model, a best-estimate one-dimensional structure is calculated in the subroutine RSLOWSPH, and this structure is used to calculate a ray-path. Travel time is accumulated along this ray path in the programs TRACKSPH and MOVEITSPH. Initial ray paths are calculated in the subroutine SRTIMESPH, which calls the program DIRECT to calculate the direct ray path, and uses the output from TTERM or TTERMS for the refracted ray times. See the procedure description of these routines for details on how all this is done. (NB: this documentation is a work in progress so only a few may be included with this release).

The hypocenter terms in the first summation are eliminated using parameter separation. If A is an m x 4 matrix containing the hypocenter derivatives, then it can be shown that

Because the last m-4 rows of [[Lambda]] are zero, only the first 4 columns of U contribute to forming A. Because the columns of U are orthonormal, the last m-4 columns thus form an annihilator for the A matrix. That is, because

Then

and the joint problem reduces to

What this means is that SPHYPIT90 does not calculate new hypocenter locations for you. That is done in a separate step, by any relocation program such as SPHREL3D90.

Teleseismic Arrivals

The arrival time for a ray from a teleseismic earthquake i to station j (T_ij) can be expressed as:

where

OTi	is the origin time of the earthquake
Bij	is the travel time from the source to the base of the model
tijl	is the time through block l
gijl	is either the time or the raylength through block l
sl	is either the fractional slowness or slowness perturbation in block l

The travel time residual is then

The average residual for all rays from earthquake i is thus

And the deviation from the mean is

This procedure eliminates the origin time from consideration. If the initial model is one-dimensional and the rays can be assumed to be parallel, then the terms involving the time through the medium vanish as well. Early versions of this program used parallel rays through the structure based on one calculated to the center of coordinates. The time to the base of the model (B_ij) was also calculated relative to the center of coordinates. Because these approximations can lead to substantial errors when regional size arrays are used, all rays in SPHYPIT90 are now traced individually using the station coordinates and the ray parameter as boundary conditions (note that there is still some redundant code left over in the routine PARSEDATA from the old days). The travel time to the base is calculated directly from the travel time table.

Note that if the variable s_l is the fractional slowness, then g_ijl is the travel time through the block, but if s_l is the slowness itself, then g_ijl is the length of the ray. In the current version of the program, slowness is used as the variable.

Weights

There are many ways to weight observations to reflect their relative uncertainty. The first is to weight by phase, using the relationship

where [[sigma]]_ij is the a priori estimate of the uncertainty in the arrival time of the phase. If phase weights are used, they can be input in either an integer or real format. If integer, some formula must be used to translate these values into real standard deviations. There are a number of ways of doing this by using the variable wttype. In particular, setting wttype to 3 or 4 allows the user to specify his/her own weight conversion function. Integer weights may be further modified by multiplicative constants that are specific to phase type (using the w(2) array) or to station (read in with the station list). They are also modified by truncation values used to eliminate outliers. The two truncation options are residual size (set by the variable trnc) and station-epicenter distance (set by the variable delmax).

If the user uses real (as opposed to integer) weights, then no additional modification is performed on the weights. The real weights that are read in are 1/[[sigma]]_ij and are squared in the program. Part of reason for this is that some datasets may be noisy and have a number of observations that are on the borderline of being labeled outliers, and can thus change their status from iteration to iteration. Such behaviour can lead to instability, and so it can be best to try to identify (and remove) outliers at the outset.

Application of weights to locally recorded data is straightforward (simple multiplication - additional discussion appears below), but the application to teleseismic data needs a bit of explanation. An appropriate way to weight teleseismic arrivals is to employ the formula:

in other words, we weight the difference between the residual and the weighted mean of the residual. It is easy to show that the origin time still drops out, as does any term that is common to all observations. In long form, we have:

Again, the t_ijl term will drop out if it is common to all rays.

Spherical Earth Calculations

All calculations in SPHYPIT90 are performed in spherical, rather than cartesian, coordinates. Spherical earth calculations can be divided into 3 parts: ellipticity corrections, earth-flattening transforms for radial variations, and spherical geometry calculations for computing horizontal distances. Because the actual ray tracing is done in 1D structures (in DIRECT), and lateral variations are calculated by following paths horizontally within each layer (in MOVEITSPH and TRACKSPH), the horizontal and vertical aspects of the ray tracing and inversion can be treated separately.

Ellipticity corrections are calculated when arrival times are read in by the subroutine PARSEDATA. Subroutine BJDAZ2 is called to determine true distance and azimuth from event to station, and ELPCR determines the appropriate correction (stored in the array elpcor(i)),from the tables of Dziewonski and Gilbert (GJRAS, 44, 7-17, 1976). The tables are read in the main routines, by a call to REDTAB. The corrections are applied in SRTIMESPH, by adding elpc to the time correction tcor.

Special note when using the relocation program: Because the ellipticity corrections are calculated only once, when the data is read, if the event moves a great deal on later iterations the corrections will be slightly wrong. For distances of a couple thousand km or less, corrections are generally a couple tenths of seconds and are relatively insensitive to reasonable hypocenter adjustments.

All depths and velocities are converted with earth flattening transforms (e.g., Chapman, 1973) upon input, to the velocity structure in INBLKMOD and to event depths in PARSEDATA. These transforms have the form

where z and v are the true depth and velocity in a spherical earth, z_f, v_f are the equivalent depth and velocity in a flat earth, and R is the earth radius (6371 km). For radially varying structure, these transforms have the property that travel times and effective ray paths between any two points are preserved; i.e. as long as the ray path is determined in a 1D structure, the same travel times are calculated if one uses spherical geometry and (z,v) variables or one uses cartesian geometry and (z_f, v_f ) variables. Thus, one can use subroutine DIRECT without modification to determine raypaths with (z_f, v_f ) variables and obtain the correct locations along the raypath for layer intersections. The transforms are calculated by functions FLATZ and FLATVEL, and the inverse transforms are calculated for locations in ROTBACK (together with horizontal coordinate transforms) and directly in SPHYPIT90 for velocities. For velocities within finite-thickness layers (all except for the basal half-space) velocities are transformed so that constant v_flyr used for the layer is the integrated average of the transformed velocity v_f within the layer:

where v_flyr, v_lyr are the transformed and true average layer velocities, and z₁ and z₂ are depths to the top and bottom of the layer. Except for within the input routines and output after the inversions, all depths and velocities are in the equivalent flat-earth coordinates, and care must be taken if any intermediate output is introduced.

All horizontal distances and intersections of raypaths with block boundaries are calculated using horizontal spherical coordinates. All raypaths are assumed to be great circles between event and station. The coordinate system and all length variables are defined differently than in the flat-earth versions: lines parallel to the Y axis are great circles, and lines parallel to the X axis are treated as small circles. The line parallel to the X axis passing through the center of coordinates is a great circle (the equator in internal coordinates). Thus, the internal set of coordinates is effectively a set of latitudes (Y) and longitudes (X) with the 0 longitude point on the equator set to the center of coordinates specified in the parameter file. If the coordinate system is rotated, the azimuth is defined as the angle between the central meridian and the true longitude line at the center of coordinates. This coordinate transform is just a series of 3 rotations stored as the double-precision coordinate transform matrix TRCO(3,3) , which is established in SETORGSPH. The transform is applied in DISTSPH, and its inverse,which is just the transform of the orthonormal rotation matrix TRCO, is applied in ROTBACK to hypocenters. Internally, all horizontal coordinates are stored as psuedo-km, equal to the local latitude/longitude coordinate in radians times the earth radius (6371). Neglecting ellipticity, these coordinates should be equal to true km at the center of coordinates and the X-coordinate will progressively deviate from true km as abs(Y) gets large, and longitude lines converge. These coordinates are scaled this way so that variables such as location errors, accumulated change in hypocenter, etc. will be approximately correct without complicated additional transforms. Horizontal hypocenter derivatives are corrected for the difference between true and psuedo-km and should give correct adjustments to hypocenters.

NOTE: for various arcane reasons, the Y variable is defined positive south (X is positive east), even though azimuth variables and variables that are sin(y) with names like sy1, sy0 are defined as if Y were positive north. Extreme caution is advised if you feel compelled to alter these internal variables.

Starting with version 2001.0404, the mode of input has changed dramatically from previous versions. This upgrade (so to speak) came as a result of a combination of having to reformat/reprogram every time a new data set came along and having to relearn how the program worked each time I applied it. The current version allows for the addition of voluminous comments so the that information-retention challenged (such as myself) can remember what is going on.

Files that contain parameter settings, file names, and format descriptors adhere to the following rules:

1. Lines beginning with a "#" are ignored.
2. Blank lines are ignored
3. Parameters/filenames/etc. are specified by a line containing the reference name (usually the same name as is used in the code) followed by a delimiter and a value.
4. Delimiters between quanities on the same line are either spaces or tabs. Multiple tabs and spaces are ignored. Strings with spaces can be treated as a single quantity by enclosing them in double quotes (").

SPHYPIT90 now sets "reasonable" defaults for variables, so it often is not necessary to specify everything. However, it is always a good idea to appreciate the meaning of all variables since your idea of reasonable and mine might not be the same.

SPHYPIT90 starts by calling a routine called SETDEFS that assigns defaults for most of the key variables and files. (details discussed below). It then asks you to specify a file that contains names of all the relevant input files, which it reads using a call to PARSEFILE. Depending on what you want to do and on the options specified, there are several possible file attachments as shown in the table below:

INPUT FILES
Reference Name	Reference File#	Actual File#	Default Filename	Description
descfile	lundes	3	default.desc	Description of the data format
parmfile	lunprm	4	default.parms	General parameters for controlling program operation
stanfile	lunsta	8	default.sta	Station locations
vel1file	lunv1d	9	default.vel1d	One-dimensional Velocity model
blckfile	lunblk	10	default.blck	Block description
slowfile	lunslo	12	default.slow	Slowness model
haltfile	lunalt	14	default.halt	Alternate Headers
datafile	lundat	20	default.data	Input data (local events)
telefile	luntel	21	default.tele	Input data (teleseismic events)
shotfile	lunsht	22	default.shot	Input data (shots)
repmfile	lunrep	30	default.reparm	Reparameterization file
icovfile	luncov	31	default.invcov	Inverse covariance matrix
dampfile	lundam	32	default.damper	Dampers for individual (reparameterized) blocks (P and S)
aprifile	lunapr	33	default.apri	A priori model (for sphypit = po)
prayfile	lunpry	40	/ernie/roecker/sphypit99/data/P.rayparms	IASP P ray parameters
srayfile	lunsry	41	/ernie/roecker/sphypit99/data/S.rayparms	IASP S ray parameters
ptimes	lunptm	42	/ernie/roecker/sphypit99/data/P.times	IASP P travel times
stimes	lunstm	43	/ernie/roecker/sphypit99/data/S.time	IASP S travel times
herrd	lunhrd	44	/ernie/roecker/sphypit99/data/herrin.dist	Herrin travel times
herrpv	lunhrv	45	/ernie/roecker/sphypit99/data/herrin.pvel	Herrin wavespeed model
velmod2	lunvl2	46	/ernie/roecker/sphypit99/data/velocity.mod2	IASP wavespeed model
elpcorf	lunelp	47	/ernie/roecker/sphypit99/data/elpcor	Elipticity corrections
OUTPUT FILES
Reference Name	Reference File#	Actual File#	Default Filename	Description
sumofile	lunsum	15	default.sumout	Summary output (SPHREL3D90 Only)
genofile	lungen	16	default.genout	Generic output
dumpfile	lundum	17	default.dump	Debugging dump
tresfile	luntrs	18	default.telres	Teleseismic residual output
syndfile	lunsyn	19	default.synout	Synthetic data output
nhedfile	lunnhd	23	default.nhead	Revised headers from sphrel3d90
pickfile	lunpck	24	default.picks	Xphase pick file
nslofile	lunnsl	34	default.nslow	Revised Slowness model
hitmfile	lunhit	35	default.hitmap	Hitmap
fracfile	lunfrc	36	default.frac	Fractional slowness output
reslfile	lunres	37	default.resol	Resolution output
totlfile	luntot	38	default.total	Total Error output
raysfile	lunray	50	default.rays	Output ray paths
proffile	lunpro	51	default.prof	Hypocenter output in Xprofile format

The number of files listed above may look daunting, but is actually quite simple. For input you need to have some arrival time data (variously in descfile, haltfile, datafile, telefile, shotfile, nheadfile, pickfile), station locations (stanfile), and a model (vel1file, blckfile, slowfile). There are a number of files that contain information that is more or less independant of the particular dataset, mostly having to do with the background Earth model (prayfile, srayfile, ptimes, stimes, herrd, herrpv, velmod2, elpcorf). A few files are needed to control how the program does the inversion (parmfile, repmfile, icovfile, dampfile, aprifile). The remaining files are all output files containing different types of information that help you to evaluate the results. Let's discuss each of these files in turn, starting with the files used to specify arrival time data.

The "filename file" is the only input file Sphypit90 requests at the standard input. All other attachments are made by Sphypit90 by associated the reference file names in the table above with the associated file. Each of these associations appears in a line containing a Reference Name and the associated file, separated by a delimiter (either a tab or a space). The following is an example of a "filename file":

####################################################################
# This is an example of a file containing the names
# of files read in by sphypit or sphrel3d90 for input 
# and output.
#
# Each line contains two fields:  the parameter name expected
# by the program and the name of the file the program should open.
# Anything after that is considered a comment.  Fields are delimited
# by spaces or tabs.
#
# Blank lines and lines beginning with "#" are ignored.  Any text
# after the second entry is ignored and can be used for comments.
#
# Output files:

sumofile	test.sumout		! Summary Output file
genofile	test.genout		! Generic Output
nhedfile	test.nhead		! New locations headers
proffile	test.prof		! Xprofile Output file
pickfile	test.pick		! Xphase Output file
# dumpfile	test.dump		! Data dump file
# tresfile	test.telres		! Teleseismic Residual file
# syndfile	test.syndata		! Synthetic data file

# Input files:

parmfile	merge.parms	! File containing the program control variables
descfile	merge.desc	! File containing the description of the data format
stanfile	park.sta	! File containing the stations list
datafile	merge.data	! File containing the data
# slowfile	test.slow	! File containing the slowness values
# haltfile	grade.heads	! Alternate location headers
# shotfile	test.shot	! Shot data
# telefile	test.tele	! Teleseismic data
# blckfile	test.block	! 3D block parameterization
vel1file	iasp.vel1d	! 1D wavespeed model
# repmfile	test.reparm	! Reparameterization file
# icovfile	test.invcov	! Inverse Covariance matrix
# dampfile	test.damper	! Dampers for meta-velocities


# Files used to calculate travel times

prayfile	/ernie/roecker/sphypit99/data/P.rayparms	! P ray parameters
srayfile	/ernie/roecker/sphypit99/data/S.rayparms	! S ray parameters
ptimes		/ernie/roecker/sphypit99/data/P.times		! P travel time tables
stimes		/ernie/roecker/sphypit99/data/S.times		! S travel time tables
herrd		/ernie/roecker/sphypit99/data/herrin.dist	! Herrin Times (No longer used)
herrpv		/ernie/roecker/sphypit99/data/herrin.pvel	! Herrin Ray Parameters (No longer used)
velmod2		/ernie/roecker/sphypit99/data/velocity.mod2	! IASP velocity model
elpcorf		/ernie/roecker/sphypit99/data/elpcor		! Ellipticity corrections
###################### END OF FILE #################################

NOTE:It's a good idea to have the first filename in this file be your preferred preferred summary file, so that the filename attachments are written to it.

There are three types of arrival time files:

1. P and S Arrivals from local earthquakes
2. P and S Arrivals from shots
3. P arrivals from teleseisms

Arrival time files consist of at least two parts: a format descriptor file (descfile) and a file containing the actual observations (datafile for local earthquakes, shotfile for shots, and telefile for teleseisms. Data for each event begins with a header line containing the location and origin time of the event, followed by a series of lines containing the observations (generally one observation per line). The last observation for a given event is followed by a blank line.

The format descriptor file (descfile) tells the program how to interpret the information in the arrival time files. The first part of the descriptor file constructs the format of the header, the second part constructs the general arrival time format. These are the header variables:

Variable Name	What it Means	Possible Values	Default
ottype	Specifies how time is represented in the header	= 0 year, month, day, hour, min, sec	1
		= 1 year, julian day, hour, min, sec
		= 2 epoch (as used in Antelope).
loctype	How Locations are specifed in the header	= 0 decimal degrees	0
		= 1 degrees and minutes
lonsign	Sign on Longitude	= 1 keep sign on longitude as is	1
		= -1 reverse sign (used for Western Hemisphere, where longitudes may be positive)
latsign	Sign on Latitude	= 1 keep sign on longitude as is	1
		= -1 reverse sign (used for Southen Hemisphere, where latitudes may be positive)
evid	Event identifier	An optional CHARACTER*6 variable used to identify the event. I use this to associate the arrivals with the headers to be sure there is no mistake about which observations go with which event. It also can be used to check if the alternate header file is in proper sequence.	EVID
vers	Location Version Number	An optional integer variable used to keep track of the hypocenter version. SPHREL3D90 increments this number by "1" when it writes a new header file.	1

These are the arrival time variables:

Variable Name	What it Means	Possible Values	Default
ttype	Specifies how an arrival time is represented	= 0 absolute year, month, day, hour, min, sec	1
		= 1 absolute year, jday, hour, min, sec
		= 2 absolute epoch (as used in Antelope).
		= 3 seconds relative to origin (i.e. a travel time; real time is found by adding to origin time)
		= 4 seconds relative to sec = 0 in header (i.e., year, jday, hour, min are the same, but secs are different)
ptype	Specifies how an Phases appear in the file	= 1 one phase per line	1
		= 2 two phases per line (P and S for a single station)
		> 2 compressed format (more than two phases per line; ptype is the number of phases per line)
fldlen	In the case that ptype > 2, this is the length of the field for a single phase.	= n length of the field	1
psign	Specifies how different phases are identified	= 0 no indicator; assume all P waves	2
		= 1 P and S on same line; P first, S second (used only if ptype=2)
		= 2 alphanumerical (e.g. 0 = P, 1 = S or just P, S)
psymb	In the case of psign = 2, specifies the symbol used to identify a P phase.	any character in a1 format	"0"
ssymb	In the case of psign = 2, specifies the symbol used to identify an S phase.	any character in a1 format	"1"
wttype	Specifies how to interpret integer or real weights	= 0 no phase weights set [phase weight = 1.0]	2
		= 1 interpret as real valued time uncertainty [phase weight = 1.0/(weight*weight)]
		= 2 interpret as real valued true weight (1/variance) [phase weight = weight]
		= 3 interpret as integer argument to the weghting function "wtfunc"[phase weight = 1./(wtfunc(iweight)*wtfunc(iweight))]
		= 4 interpret as integer argument to wtfunc, but add "1" to the integer first (to allow for a value of "0"). As above, with iweight => iweight+1.
evid	Event identifier	An optional CHARACTER*6 variable used to associate the arrivals with the proper event. Note that this evid is read in but not currently used in the software - it's a bookeeping device for the user.	(NONE)

As with the location header, first specify each of these main format variables, followed by the positions and formats of the parameters themselves (like sta, time, weight, etc. as required)

After these definitions, we tell the program where to read in the variables it needs. This is done by specifying the variable name, the column where the number begins, and the format to use. For example, suppose that the latitude of the event begins in column 20 and should be read in as f8.4, then the line

xlat 20 f8.4

Variable Name	What it Means	Default
lat	Integer Degrees of Latitude (if loctype = 1)	(None)
lon	Integer Degrees of Longitude (if loctype = 1)	(None)
xlat	Decimal degrees of Latitude (if loctype = 0) or minutes of latitude (if loctype = 1)	34.0
xlon	Decimal degrees of Longitude (if loctype = 0) or minutes of longitude (if loctype = 1)	-121.0
depth	Depth of Event (usually in kilometers)	0.0
year	Year of Origin Time (if otype = 0 or 1)	2001
mon	Month of Origin Time (if otype = 0)	1
jday	Day of Origin Time (if otype = 0) or Julian Day of Origin Time (if otype = 1)	1
ihr	Hour of Origin Time (if otype = 0 or 1)	1
min	Minute of Origin Time (if otype = 0 or 1)	1
sec	Second of Origin Time (if otype = 0 or 1)	0.0

The list of header variables is terminated with the single line beginning with the word "END", e.g.,:

END

Variable Name	What it Means
year	Year of Arrival Time (if otype = 0 or 1)
mon	Month of Arrival Time (if otype = 0)
jday	Day of Arrival Time (if otype = 0) or Julian Day of Arrival Time (if otype = 1)
ihr	Hour of Arrival Time (if otype = 0 or 1)
min	Minute of Arrival Time (if otype = 0 or 1)
sec	Second of Arrival Time (if otype = 0 or 1)
sta	Station code (character)
phase	Phase type
weight	Weight

#
#  Data format descriptor
#
#	Location headers must contain, at minimum, the
#	origin time and 3d location of the hypocenter.  
#
#	A line of data must contain, at minimum, a station
#	code and a time.  It can contain additionally a weight
#	and a phase identifer.  If the identifier is missing, it
#	is presumed to be a P wave.  If the weight is missing, it
#	is given the default value of "1".  Lines of data may also
#	contain left and right uncertainties, but this has not
#	been incorporated into the code as yet.
#
#
# Section for describing the format of the location header
#
#	ottype	0	year, month, day, hour, min, sec
#	ottype	1	year, jday, hour, min, sec
#	ottype	2	epoch

ottype	1

year	1	i4
jday	6	i3
ihr	10	i2
min	12	i3
sec	16	d8.4

# Location part of header
xlat	25	f9.5
xlon	34	f11.5
depth	46	f6.2

# Event id - must be a character
evid	53	a6
vers	59	i4
END
#
# Example of a header described by the above (preceded by character counter line).
# 2000 324  1 28  42.7807  36.00000 -120.50000   1.00    802   1 
#
# Example of a phase line described by the data format specs below:
# MUST  2000 324  1 28  52.781 0  -0.200   0.200    802 
#
# Section for describing the format of the data line
#
#	ttype	0	absolute year, month, day, hour, min, sec
#	ttype	1	absolute year, jday, hour, min, sec
#	ttype	2	absolute epoch
#	ttype	3	relative year, month, day, hour, min, sec
#	ttype	4	relative year, jday, hour, min, sec
#	ttype	5	relative epoch
#
#	ptype	1	one phase per line
#	ptype	2	two phases per line (P and S for a single station)
#	pytpe	>2	compressed format (more than two phases per line)
#
#	psign	0	no indicator; assume all P waves
#		1	P and S on same line; P first, S second (only if ptype=2)
#		2	alphanumerical (e.g. 0 = P, 1 = S or just P, S)
#			In this case follow the number 2 with the two symbols used
#
#	psymb	a1	The character used to identify P waves
#	ssymb	a1	The character used to identify S waves
#
#	wttype	0	no weights set; use a uniform weight
#		1	interpret as time uncertainty
#		2	interpret as weight (1/variance)
#		3	interpret as argument to wtfunc
#		4	interpret as argument (+1) to wtfunc
#
ptype	1
psign	2	
ttype	1
wttype	1
psymb   0
ssymb	1

sta	1	a4
year	7	i4
jday	12	i3
ihr	16	i2
min	19	i2
sec	23	d6.3
# Note phase should be character format
phase	30	a1
weight	40	f7.3

evid	48	a6


END OF DESCRIPTOR

NOTES:

1. If you use the epoch option, you might want to specify the format for "sec" as double precision ("d" format). Otherwise you could suffer some roundoff errors.

Seismograph station codes and locations appear in a separate file. Unlike the arrival time files, however, the format descriptor is prepended to the station list. Conventions for reading in the format descriptor are the same as those for the arrival times. The following table summarizes the pertinant variables:

Variable Name	What it Means	Possible Values	Default
loctype	How Locations are specified	= 0 decimal degrees	0
		= 1 degrees and minutes
lonsign	Sign on Longitude	= 1 keep sign on longitude as is	1
		= -1 reverse sign (used for Western Hemisphere, where longitudes may be positive)
latsign	Sign on Latitude	= 1 keep sign on longitude as is	1
		= -1 reverse sign (used for Southen Hemisphere, where latitudes may be positive)
ieltype	Elevation convention	= 0 elevation is in (integer) meters above sea level	0
		= 1 elevation is in (real) kilometers meters above sea level

These are the individual station variables:

Variable Name	What it Means	Default
sta	Station code (4 characters max)	UNKW
iw	Integer station weight. The interpretation of this integer weight depends on the parameter iform (see also discussion in Parameter File below) as follows: iform = 0; weight = 1.0/(2.0**float(iw)) iform = 2; weight = 1.0 iform = anything else; weight = 1.0/real(iw*iw)	1
lat	degrees of latitude (loctype = 1)	(None)
lon	degrees of longitude (loctype = 1)	(None)
xlat	minutes of latitude (loctype = 1) or decimal degrees of latitude (loctype = 0)	34.0
xlon	minutes of longitude (loctype = 1) or decimal degrees of longitude (loctype = 0)	-121.0
ielev	Elevation in meters above sea level	0.0
ptc	P time correction (seconds)	0.0
stc	S time correction (seconds)	0.0
fl	The thickness of a station specific layer (in km) beneath the station. Note that this layer should not be so thick that it goes into the second layer. Not used if set equal to 0.	0.0
ps	P velocity in sediment layer. If this is zero, it is set equal to the velocity in the first layer.	vp(1)
ss	S velocity in sediment layer. If this is zero, it is set to ps times the Vs/Vp ratio in the first layer.	ps*vs(1)/vp(1)

Notes:

1. Locations less than one degree west of the prime meridian should be entered as 0 degrees and negative minutes. For example, 0'45"W should appear as 0-45.00.

2. Locations less than one degree south of the equator should be entered as 0 degrees and negative minutes. For example, 0'45"S should appear as 0-45.00.

3. It is possible to use more than one format in a single station list file by separating the different parts of the file with a "!" in the first column. Schematic example:

   [Description of 1st format]
   [Stations listed in 1st format]
   !
   [Description of 2nd format]
   [Stations listed in 2nd format]
   

   ... and so on

############################################################
# This is an example of a station list.  The format
# in this case is based on the UCSD stapar file

# Station Code
sta	1	a4
# Station Weight
iw	6	i1
# Latitude (North Positive)
xlat	13	f9.5
# Longitude (East Positive)
xlon	24	f10.5
# Elevation (Meters above sea level)
ielev	41	i4
# P correction
# ptc
# S correction
# stc
# Substation Layer thickness
# fl
# Substation Layer P wave speed
# ps
# Substation Layer S wave speed
# ss
END

#23456789012345678901234567890123456789012345678901
BUMP 1    1  35.98058  -120.55225        751
CANS 1    2  35.90038  -120.56625        669
################### END OF FILE ###########################

The 1D Velocity Model file serves two purposes:

1. It can be used to specify 1D model as a starting point for an iterative 3D inversion.

2. It is used to specify the depths of horizontal boundaries in a 3D model. Note that if we input a 3D model, the velocities specified in the 1D file are not used for anything (but the depths of the interfaces are). The use of 1D velocities is decided by the value of the variable iasn (see discussion below).

The 1D Model file has the same makeup as the station list file: start with a format descriptor followed by the specification of variables. End the descriptor with a line containining the work "END". There are three variables that need to be read in:

Variable Name	What it Means
vp	P wavespeed (usually in km/s)
vs	S wavespeed (usually in km/s)
ht	Depth to the top of the layer

Important Note on "ht": Because the station elevations are usually given with respect to sea level, the origin "zero" is usually at sea level as well. Thus, negative ht puts interfaces above sea level (postive depth is down). In order for the rays to reach any given station, the station needs to be imbedded in the model. Thus, it is important to be sure that the initial value of ht puts it above the elevation of the highest station. For example, if the highest station is at 4.3 km elevation, then ht must be less than -4.3 (absolute value > 4.3). Moreover, ht should appear in the input file in lines from least (shallowest) to greatest (deepest) values. At present, the program does not check for out-of-sequence depths, so be sure to check this before starting.

The three variables appear in a one-line-per-layer form. The actual number of values read in depends on the parameter nl (number of layers) specified in the parameter file (see below). If the value of nl exceeds the number of lines in this file, nl is redefined to be this number of lines.

######### This is an example of a 1D wavespeed file.  ##########

# P wavespeed
vp	16	f5.3
# S wavespeed
vs	26	f5.3
# Depth to top of layer
ht	 1	f8.2
END
# The following is from the IASP91 file
#23456789012345678901234567890123456789012345678901
   -5.00       5.800     3.360     2.600    1752.1     778.7
   20.00       6.500     3.750     3.626    1412.2     627.6
   35.00       8.040     4.470     3.388    1449.1     601.4
   71.00       8.044     4.483     4.149    1464.1     607.1
  120.00       8.050     4.500     4.186     195.0      80.0
  210.00       8.300     4.518     4.254     195.8      80.3
  271.00       8.523     4.628     4.299     365.2     143.0
  410.00       9.030     4.870     4.404     365.8     143.0
###################### END OF FILE #############################

The 3D model is described by a set of constant velocity volumes that fill the model space. In their simplest form these volumes are rectangular prisms, but through a procedure called "reparameterization" we can join several smaller volumes into larger ones, thus producing regions with a complicated shape. We therefore need to speak of two types of volumes: the elementary volumes (the mesh of small rectangular prisms) and the reparameterized volumes (made of groups of elementary volumes).

Specifying the Elementary Volumes

The elementary volumes are specified by three sets of perpendicular interfaces: two vertical and one horizontal. The position of the horizontal interfaces is specified in the 1D velocity model file (described above) and values are usually in km w.r.t. sea level. The positions of vertical interfaces are specified w.r.t. a center of coordinates (variables olat and olon set in the parameter file; see discussion below) and can be different in each layer. As discussed above the default positive directions are south and east, but this can be modified by specifying an azimuth of rotation in the parameter file. Not that the rotation is defined in a counterclockwise sense, so that 20^o rotation means that an interface originally pointing to the north is now pointing to N20^o W.

The vertical interface positions appear in the Block Description File (blckfile). This file contains three entries for each layer. The first entry has four numbers:

nxs	The number of N-S oriented interfaces
nys	The number of E-W oriented interfaces
pnit	The layer scaling factor for P data
snit	The layer scaling factor for S data

The next two entries contain the actual positions of the interfaces:

bx(i,j)	The displacement of the N-S interface i in layer j (East positive)
by(i,j)	The displacement of the E-W interface i in layer j (South positive)

The following is an example of a block description file for a three layer model:

  31  19    1.0      1.0         layer 1
 -459.23 -396.61 -313.11 -292.24 -271.36 -250.49 -229.61 -208.74 -187.87 -166.99
 -146.12 -125.24 -104.37  -83.50  -62.62  -41.75  -20.87    0.00   20.87   41.75
   62.62   83.50  125.24  166.99  208.74  271.36  292.24  313.11  354.86  396.61
  438.35 
 -277.65 -222.12 -194.36 -166.59 -138.83 -111.06  -83.30  -55.53  -27.77    0.00
   27.77   55.53   83.30  111.06  138.83  166.59  194.36  222.12  277.65
  31  19    1.0      1.0         layer 2
 -459.23 -396.61 -313.11 -292.24 -271.36 -250.49 -229.61 -208.74 -187.87 -166.99
 -146.12 -125.24 -104.37  -83.50  -62.62  -41.75  -20.87    0.00   20.87   41.75
   62.62   83.50  125.24  166.99  208.74  271.36  292.24  313.11  354.86  396.61
  438.35 
 -277.65 -222.12 -194.36 -166.59 -138.83 -111.06  -83.30  -55.53  -27.77    0.00
   27.77   55.53   83.30  111.06  138.83  166.59  194.36  222.12  277.65
   19  17    1.0      1.0         layer 3
 -459.23 -396.61 -354.86 -313.11 -271.36 -208.74 -166.99 -146.12 -104.37  -62.62
  -41.75  -20.87   20.87  150.00  271.36  313.11  354.86  396.61  438.35 
 -277.65 -166.59 -138.83 -111.06  -83.30  -55.53  -27.77    0.00   27.77   55.53
   83.30  111.06  138.83  166.59  194.36  222.12  277.65

Specifying the Reparameterized Volumes

Once we have specified the elementary volumes, we can join subsets of them together by reparameterizing. We do so by creating a table that appears in the repmfile file. The idea is this: assign a distinct number to each elementary volume that you wish to join together, in the order that the elementary volumes appear (starting in the NW corner to the NE corner, with the SE corner being the last entry). An example that uses the above interfaces appears below. Note that in each case there is a "dummy" line for inserting the layer number (formated as i5 - this is to help you keep track of where you are). As above, the formats here are hard wired and the reparmeterization numbers are read in as (20i5).

    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    2    2    2    2    2
    2    2    2    2    2    2    2    2    2    2
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    4    4    4    4    4
    4    4    4    5    5    5    5    5    5    5
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    3    3    6    6    6
    4    4    4    5    5    5    5    5    5    5
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    3    3    3    3    6
    6    6    6    7    7    7    8    8    8    8
    9    9   10   10   10   10   10   10   10   10   10   11   11   11   11   12   12   12   12   13
   13   13   14   14   15   15   16   16   16   17
    9    9    9    9   10   10   10   10   10   10   10   11   11   11   11   18   18   18   18   19
   19   19   20   20   21   21   22   22   17   17
    9    9    9    9    9    9    9   10   10   10   10   24   24   24   24   25   25   25   25   26
   26   26   23   23   27   27   22   17   17   17
   28   29   29   29   29   29   29   30   30   10   10   24   24   24   24   31   31   31   31   32
   32   32   33   33   33   17   17   17   17   17
   28   29   29   29   29   29   29   30   30   30   30   30   35   35   35   35   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   29   29   29   29   37   37   38   38   30   30   30   60   60   35   35   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   39   40   40   37   37   38   38   38   60   60   60   60   60   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   40   40   40   40   37   37   37   38   38   38   60   60   60   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   40   40   40   41   41   41   41   42   42   42   63   63   64   64   64   45   45
   45   45   45   45   45   45   45   45   45   45
   28   48   48   41   41   41   41   42   42   42   42   63   63   63   63   64   64   64   45   45
   45   45   45   45   45   45   45   45   45   45
   47   48   48   49   49   49   49   42   43   43   43   43   64   64   64   64   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   51   51   52   52   52   52   43   43   44   44   44   44   50   50   50   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   51   51   52   52   52   52   44   44   44   50   50   50   50   50   50   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   53   53   53   53   53   53   54   54   54   54   54   54   54   54   54   54   54   54   54
   45   45   45   45   45   45   45   45   45   45
    2
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    2    2    2    2    2
    2    2    2    2    2    2    2    2    2    2
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    4    4    4    4    4
    4    4    4    5    5    5    5    5    5    5
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    3    3    6    6    6
    4    4    4    5    5    5    5    5    5    5
    1    1    1    1    1    1    1    1    1    1    1    3    3    3    3    3    3    3    3    6
    6    6    6    7    7    7    8    8    8    8
    9    9   10   10   10   10   10   10   10   10   10   11   11   11   11   12   12   12   12   13
   13   13   14   14   15   15   16   16   16   17
    9    9    9    9   10   10   10   10   10   10   10   11   11   11   11   18   18   18   18   19
   19   19   20   20   21   21   22   22   17   17
    9    9    9    9    9    9    9   10   10   10   10   24   24   24   24   25   25   25   25   26
   26   26   23   23   27   27   22   17   17   17
   28   29   29   29   29   29   29   30   30   10   10   24   24   24   24   31   31   31   31   32
   32   32   33   33   33   17   17   17   17   17
   28   29   29   29   29   29   29   30   30   30   30   30   35   35   35   35   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   29   29   29   29   37   37   38   38   30   30   30   60   60   35   35   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   39   40   40   37   37   38   38   38   60   60   60   60   60   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   40   40   40   40   37   37   37   38   38   38   60   60   60   35   35   35   35
   35   36   36   36   36   17   17   17   17   17
   28   39   39   40   40   40   41   41   41   41   42   42   42   63   63   64   64   64   45   45
   45   45   45   45   45   45   45   45   45   45
   28   48   48   41   41   41   41   42   42   42   42   63   63   63   63   64   64   64   45   45
   45   45   45   45   45   45   45   45   45   45
   47   48   48   49   49   49   49   42   43   43   43   43   64   64   64   64   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   51   51   52   52   52   52   43   43   44   44   44   44   50   50   50   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   51   51   52   52   52   52   44   44   44   50   50   50   50   50   50   50   50   50   50
   45   45   45   45   45   45   45   45   45   45
   47   53   53   53   53   53   53   54   54   54   54   54   54   54   54   54   54   54   54   54
   45   45   45   45   45   45   45   45   45   45
    3
    1    1    1    1    1    1    1    1    1    1    2    2    2    2    2    2    2    2
    3    4    5    5    5    5    5    5    5    5    6    6    6    6    6    6    6    6
    3    4    4    5    5    5    5    5    5    5    6    6    6    6    6    6    6    7
    3    4    4    4    5    5    5    5    5    5    6    6    6    6    6    6    7    7
    3    4    4    4    4    8    8    8    8    8    8    8   18   18   18    7    7    7
    3    4    4    4    4    4    8    8    8    8    8    8   18   18    7    7    7    7
    3    9    9    9    9   29   29   29    8    8    8    8   18   18    7    7    7    7
    3    9    9    9    9   29   29   29   29   10   10   10   10   10    7    7    7    7
   11   12   12   12   12   29   29   29   29   29   10   10   10   10    7    7    7    7
   11   12   12   12   12   12   12   29   29   29   10   10   10   10    7    7    7    7
   11   12   12   12   12   12   12   23   23   23   23   14   14   14   14   14   14   14
   11   12   12   12   12   23   23   23   23   23   23   14   14   14   14   14   14   14
   15   13   13   13   13   23   23   23   23   23   16   16   14   14   14   14   14   14
   15   13   13   13   13   23   23   16   16   16   16   16   14   14   14   14   14   14
   15   13   13   13   13   16   16   16   16   16   16   16   14   14   14   14   14   14
   15   16   16   16   16   16   16   16   16   16   16   16   14   14   14   14   14   14

Notes:

1. The numbers used to label common reparameterized volumes are arbitrary, but can depend on whether reparameterization is done in a layer-specific or volume-specific manner. If we choose layer-specific (only blocks within a given layer are grouped together) then the numbers used in reparameterization need be unique only within each layer. If we choose volume-specific (blocks can be grouped across layers) then the numbers are unique for the entire model and are not reset in each layer. The choice of layer- vs. volume-specific reparameterization is made with the parameter ireparm (see discussion below). Note also that in the case that iasn = 1, each meta-block in a given layer is assigned the same velocity.

2. Even though the outer interfaces are specified, the program extends them out to infinity so that there is no "edge" to the model (only a bottom and top). The reason for including them is for auxilliary plotting routines that use these same files. Note that the program will not allow a ray to enter from the side of the model, which generally occurs when the ray is from a not-too-distant teleseism and bottoms out at a depth shallower than the bottom of the model.

3. 1D inversion can be accomplished by setting the number of lateral interfaces to 2. The actual distances used don't matter as in this case each block extends to infinity in all horizontal directions (and hence is 1D).

Specifying the Velocities

There are two files where models can be read in: one is the slowness model for the current iteration (p_k) and the other is the a priori slowness file (p₀). Both of these are in the same format and specify a slowness for each global (reparameterized) element. Currently, P and S waves must have the same structure, so an equal number of slownesses must be read in for each (P is first, S is second). Note that if you use "reparameterization" you should specify only as many velocities as you have "meta" variables.

The values of the velocities (actually reciprocal velocities or slownesses) in the current 3D model appear as a series of numbers in the Slowness Model File (slowfile). The format is (8f10.7). Here is an example of a few lines taken from one of these files:

  .1863099  .2051207  .2296407  .1992282  .1968406  .1891249  .2135800  .2111019
  .1992003  .2034132  .1952003  .1957237  .2334384  .2216818  .1946827  .2342241
  .2054453  .1878531  .1966237  .2262344  .2082787  .1921463  .1960050  .1894153

It is important to remember that these numbers correspond to the reparameterized model, and not the elementary model, and the order is in the order of reparameterization. Generally speaking this is the order in which new reparameterized volumes are encountered by the program. All the P slownesses are specified first, followed by the S slownesses.

SHYPIT90 reads in a number of tables; most of these are for teleseismic inversion which require travel times and ray parameters from the source to the base of the model. Historically there are eight of these files:

prayfile	IASP P ray parameters
srayfile	IASP S ray parameters
ptimes	IASP P travel times
stimes	IASP S travel times
herrd	Herrin travel times
herrpv	Herrin wavespeed model
velmod2	IASP wavespeed model
elpcorf	Elipticity corrections

Each of these tables is included as part of this distribution. You only need to put them somewhere on your system and then specify the path in the "filenames" file (see below for a discussion of this file). The Herrin tables are no longer used in this version but are included anyway.

Note that because these files are static they may be hardwired in SETDEFS to some very machine-specific location. You may want to edit these locations in SETDEFS prior to compilation.

The Parameter File (parmfile) is a set of switches and initial values that control the operation of the program. The values of parameters are specified one per line, typically as the parameter name and value separated by a delimiter (space or tab). The actual order in which the parameters are specified is inconsequential, and most have default values. In the current version there are 58 of these that the user can specify. For clarity I will group these together by function.

Parameters Related to Coordinate Reference Frame
Variable Name	What it Means	Default Value
olat	Decimal Latitude of Center of Coordinates	35.00
olon	Decimal Longitude of Center of Coordinates	-121.00
oz	Depth of Center of Coordinates	0.0
az	Rotation of Coordinate system, counterclockwise w.r.t. north	0.0

Notes:

1. It is generally a good idea to pick the center of coordinates to be somewhere near the areal center of the volume you are trying to image.

2. Generally you should choose oz to be 0. (= sea level)

3. az is the counterclockwise rotation of the earth model with respect to north. Thus a 20^o rotation means that an interface originally pointing to the north is now pointing to N20^o W.

Parameters Related to Weighting
Variable Name	What it Means	Default Value
w1	Relative weight assigned to all P arrivals	1.0
w2	Relative weight assigned to all S arrivals	1.0
wfac	Relative weight assigned to all arrivals	1.0
jwt	If this is 1, all teleseismic weights are reset to 1.0. Otherwise, phase specific weights can be assigned.	0
wt1, wt2,..., wt8	These are the 8 values of wtfunc, the eight element array that tells the program how to interpret integer weights as standard errors. For example, the sequence (0.1, 0.2, 0.3, 0.4,...) tells the program to interpret a zero phase weight as a reading with an uncertainty of +/-0.1s, a phase weight of 1 as +/-0.2s and so on. This function is used if wttype = 3 or 4.	wt1 = .01 wt2 = .02 wt3 = .04 wt4 = .06 wt5 = .08 wt6 = .10 wt7 = .15 wt8 = .20
delmax	This is the truncation value for long distances. Data associated with rays that travel more than this distance are weighted down severely. Note that if delmax = 0 then it is not used at all as a weighting factor.	50.0
trnc	This is the truncation value for big residuals. Data with residuals above this amount are weighted down severely	10.0
trnc2	residual truncation value used in Sphrel3d90	1.0
iform	Controls how station weights are used, and causes phase weights to be interpreted literally (i.e., not dynamically weighted) if set to 2: = 0 integer station weights (iw) interpreted as weight = 1.0/(2.0**float(iw)) = 2 real phase weights are used, i.e., the weights are taken literally and no "dynamic weighting" is applied. In this case the station weight is set to 1.0. = anything else; station weight = 1.0/real(iw*iw)	2
iwgt	This is a weighting option that allows for weights to be applied on a layer by layer basis. If iwgt = 0 no weights are applied by layer, if iwgt = 1, then weights are read in from the parameter file and applied, and if iwgt is anything else, weights are calculated according to the thickness of the layer.	0

Notes:

1. All of the weights applied are simply multiplied together unless real phase weights are used, in which case none of these extra weights are applied. Application of weights is "dynamic" meaning for example that when locating an event the weight of a phase can change as its residual or distance from source to receiver changes.

2. There are two residual weights: trnc and trnc2, which weight down residuals that exceed these values. These numbers are interpreted differently by SPHYPIT90 and SPHREL3D90. In SPHYPIT90, trnc is applied to local events and trnc2 is applied to teleseisms. In SPHREL3D90, trnc is applied for the first several iterations, and this value is gradually replaced by trnc2 as the number of iterations approaches the maximum. The philosophy here was to allow a lenient truncation value when the hypocenter was poorly known, and then gradually require the residuals to be smaller as iterations progressed. This part of the code was written when CPU time was expensive. Presently I generally just set trnc and trnc2 to be the same value and gradually decrease both of them on subsequent runs of SPHREL3D90. THE USER SHOULD BE MINDFUL THAT SPHYPIT90 and SPHREL3D90 TREAT trnc2 DIFFERENTLY WHEN DOING A COMBINED LOCAL/TELSEISMIC INVERSION

3. The actual weighting applied by trnc, trnc2, and delmax is an exponential decay near the value of the parameter (as opposed to a sharp cutoff). The weight function wf is 1.02/(1. + .02*exp(rx)) where rx = res*res/trc. res is the residual (or distance) and trc is t*t/(alog(50.*(exp(1.0) - 1.0)), where t is either trnc, trnc2, or delmax.

Parameters Related to the Data Specification
Variable Name	What it Means	Default Value
neqs	The number of local earthquakes in the data file. If this number is greater than or equal to the the actual number in the data file, the program will use all of the local data	1
nshot	Number of shots to use in the inversion	0
ntel	The number of teleseismic earthquakes in the data file. If this number is greater than or equal to the the actual number in the data file, the program will use all of the teleseismic data	0
ihead	= 1 input starting headers from another file	0
mnd	If this parameter is 0, then the program uses the origin time in the data file as the real origin time, and travel time is calculated relative to this time. With any other value the origin time is presumed to be 0 (i.e., the data are in travel times)	0

Parameters Related to the Inversion
Variable Name	What it Means	Default Value
idmptp	This parameter controls the type of damping applied, as follows: = 0 Use the same damper for all parameters (= vthet) = 1 Use the same damper for all variables in a layer. In this case, layer specific standard deviations of slowness ([[sigma]]) are read in after the station list. These numbers are converted to dampers by determining 1/[[sigma]]2 = 2 Use a different damper for each variable. In this case, parameter specific standard deviations of slowness are read in from a separate file (dampfile) and are converted to dampers in the same way as described above.	0
vthet	This is a global slowness damper, and is applied only if the parameter idmptp is equal to 0. Note that vthet is added directly (it is not redefined) and thus should be the inverse of the model covariance.	1.0
iasn	If this equals 1, the one-dimensional structure in the parameter file is used as the starting model, if it equals anything else, a three-dimensional slowness model is read in from slowfile.	0
kfirst	This allows layers near the surface to be excluded from the inversion. For example, if kfirst = 2, then the first two layers are not included.	0
klast	This allows the layers near the bottom to be excluded from the inversion. For example, if klast = 5, then layers 5 and greater are not included. If klast is 0 or set less than or equal than kfirst it is disabled by redefining it as one more than the number of layers.	0
nmin	the minimum number of hits required for a block to be included in the inversion	1
irflag	This parameter allows a variety of inverse covariance matrices to be used. If irflag = 0, then no inverse covariance matrix is read in. If irflag = 1, then layer-specific (no coupling between layers) covariance matrices are read in. If irflag = 2, then the "full" covariance matrix is read in. These matrices can be generated by the utility programs rmm.f and rmmfull.f	0
nkr	This parameter allows the user to define which blocks to perturb, and which to keep constant. If nkr > 0, a list of nkr selected blocks is read in. If nkr = 0, then all eligible blocks are solved for.	0
mstor	If nkr > 0, then this array stores the numbers of the blocks that will be included in the inversion. The numbering scheme is the same as the global numbering scheme. mstor is the only multi variable input parameter. To input the values to the mstor array, simply list them after the variable name. Here is an example: nkr 5 mstor 100 101 203 300 301 There should be nkr numbers appearing after mstor.	(None)
iusepo	If = 1, then the program will read in the starting model and apply the Po constraint.	0

Notes:

1. Some of these parameters were introduced when CPU was slow and memory was small, but nevertheless can serve some useful functions. nmin was used to exclude blocks that would not have been well resolved anyway, but now I generally just set it to 1 (include all blocks) and ignore the poorly resolved blocks after the fact. kfirst, klast, and nkr were also introduced to reduce the size of the inverse problem, but now can be used advantageously to test the robustness of the model. For example, one might want to test the influence of crustal structure on a teleseismic inversion by fixing the velocities in the crust.

2. The value of the dampers, such as vthet have a theoretical meaning based on the ratio of the data uncertainty to the model uncertainty. For example, if the data is poor and the model is well known, the damper will have a large value (which tends to keep the corrections small). Practically speaking, dampers will have little effect if they are not of comparable magnitude to the size of the diagonals of the normal equations. The user is advised to examine the sizes of these values when choosing an appropriate damper. It is also a good idea to use large dampers at the first iteration or so to quell the influence of nonlinearities.

Parameters Related to the Model specification
Variable Name	What it Means	Default Value
nl	Number of layers (horizontal interfaces). If this number is greater than the number of layers in the model file, it is redefined accordingly.	4
ireparm	This allows small blocks to be grouped together in "meta-blocks", and thus allow for complicated geometries without introducing a lot of parameters. = 0 Don't reparameterize (each block is separate) = 1 Reparameterize by layer = 2 Global reparameterization	0

Parameters Controlling the Output
Variable Name	What it Means	Default Value
idatdmp	= 1 dump a variety of data specific information to a file = 0 don't dump	0
iray	= 1 write out raypaths to file raysfile.	0
iprof	= 1 output the headers in a file that can be interpreted by the program Xprofile.	1
ipick	= 1 output the data file in a file that can be interpreted by the program Xphase.	0
ioutflg	= 0 The default SPHYPIT90 output (simple) = 1 Include resolution and covariance of adjoining blocks, and percent of power in the diagonal = 2 Find block number and magnitude of largest off-diagonal term for resolution and covariance. = 3 All of the above = Above +4 Also output teleseismic residuals = Above +8 Also output teleseismic residuals and weights actually used	0
maksyn	If this parameter is 1, then the program will generate an output file containing the calculated travel times for all of the local and teleseismice data that was used in the analysis. If set to 0, no such file is created. Note that for cases when the user wants only the synthetic data file to be created (as opposed to doing a full inversion) it would be most efficient tu set ires = -1 in conjunction with maksyn = 1.	0
ires	This controls the output of the program, as follows: = 0 Don't calculate resolution and error = 1 Print resolution and error summary maps (diagonals only) = -1 Stop after printing a hit map = -2 Stop after printing diagonals of normal equation matrix	0

Miscellaneous Parameters
Variable Name	What it Means	Default Value
ichkhd	If = 1, Check the alternate header file against the one in the data file for consistency.	0
comment	Output a comment to the summary output file.	0
iformt	This parameter appears in the code but is now obsolete	(None)
ihypflg	This parameter appears in the code but is now obsolete.	0

Parameters Used by Sphrel3d90 (but not Sphypit90)
Variable Name	What it Means	Default Value
isw1	= 2 for skipping correction output = 3 collect global statistics on events	2
isw2	= 1 for adding a damper	1
isw3	= 1 to skip output of residual map	0
isw4	= 0 to free the depth = 1 to fix the depth by resetting the correction to zero = 2 to fix the depth by removing the depth eqn. > 2 fix the depth as above for n-2 iterations	0
isw5	= 1 to calculate station stats = 2 for full phase (data) file output to header file = 3 make a phase file with absolute weights	0
isw6	(Not used in this version)	1.0
ittmax	Maximum number of iterations.	10
depset	Minimum depth allowed. If the depth is corrected to some value less than depset, it is set equal to depset.	-2.0

Notes:

1. ISW normal settings : 121000 for normal relocation and 220000 for a residual map

2. SPHREL3D90 also checks to make sure that a corrected depth does not go above the highest altitude station, so can supercede depset if it is set to too high an elevation.

There are various options for adding on an inverse covariance matrix (or damper) to the weighted normal equations. The simplest one is to add a constant (vthet) to all diagonals. We can make the damper a vector instead of a scalar to allow for different uncertainties in variables. There are two other options: the first is to assign dampers on a layer by layer basis, the second is to allow a different damper for each variable. To exercise the first option a line of dampers is read in (the first nl dampers for P velocities, the second nl for S velocities) in the parameter file. For the second option a separate file is opened up (file 34) and the dampers read in.

More complicated damping is done by inverting the covariance matrix with off-diagonal elements. This can also be done on a layer-by-layer basis, or a full, general sense. Because of the time involved, the matrix inversion is now done in a separate step (in programs like RMM and RMMFULL) and the inverse matrix is read in from a separate file (icovfile).

How the Data Weights are Used

Theoretically, all data weights appear as part of the data covariance matrix (C_dd ). Current versions of SPHYPIT90 do the matrix multiplications

and

by assuming that the data errors are uncorrelated (C_dd is diagonal, and the elements are equal to the variance of the observation in question), and multiplying the matrices G, G^T, and the vector d_o - g(p_k) by one over the square root of the diagonal. In order to do this correctly then, the weights in the wt array should be equal to one over the presumed variance of the observation (the square root is taken in the subroutine BUILDG2 in SPHYPIT90 or MATRIX in SPHREL3D90).

Again, the makeup of C_dd depends on the value of iform and wttype. If iform = 2, then only the real weights in the data file are applied. If wttype = 2 and iform != 2, then the data is subject to station weights, global and phase-type weights, and cut-off weights for distance and residual. If iform = 1, then individual phase weights (converted from integer) are applied as well. All of these weights are scalars that multiply the terms in C_dd^-1.

There are several output files. Most of the output is self-explanatory. It is either in the form of regurgitated variables or in maps of various things like velocities, resolution, error, and so on.

Reference Name	What it is	Description of Contents
sumofile	Summary output	This is where most of the regurgitated input information winds up.
genofile	Generic output	where most of the newly created output goes. For SPHREL3D90 the iterative solutions appear here.
dumpfile	Debugging dump	Basically this is a lot of output that could be helpful in figureing out what is going wrong. Too much to describe here.
tresfile	Teleseismic residual output	Just like it sounds.
syndfile	Synthetic data output	The Synthetic Data file is generated when the variable "maksyn" is set to 1. This file contains calculated travel times for all of the local and teleseismic data that was used in the analysis. This file can be very useful in testing the program or for hypothetical tests of the final model.
nhedfile	Revised headers from SPHREL3D90 are writen to this file.
pickfile	Xphase pick file	Used by Xphase to plot arrivals on waveforms.
repmfile	Reparameterization file	Contains the reparameterization table
nslofile	Revised Slowness model	This is the new slowness model that can be used in SPHREL3D90 for relocation.
hitmfile	Hitmap	automatically created by SPHYPIT90 and contains integer number of hits/metablock in the same ordering scheme as the slowness file.
fracfile	Fractional slowness output	Just like it sounds.
reslfile	Resolution output	The resolution matrix is calculated in the usual way, and to my mind the best interpretation of its elements is that suggested by Jackson and Matsu'ura (1985). Generally, resolution is a measure of how much information you are getting from your observations relative to how much you are getting from your a priori constraints. Diagonals greater than 0.5 mean that more information is coming from the data, less than 0.5 means that more is coming from the a priori constraints. Poor resolution basically means you aren't learning anything new from your observations.
totlfile	Total Error output	The a posteriori model covariance matrix is calculated from equation (44) of Tarantola and Valette (1982): In older versions of SPHYPIT90 two error matrices were calculated. The first one was called "standard error" and was derived from the formulation of Aki and Richards. It is not consistent with the above equation and has therefore been dropped in more recent versions. The second matrix is called the "total error" and corresponds to the a posteriori model covariance shown above.
rayfile	Ray Path output	Ray segments that can be plotted to show ray coverage. Output if irayp=1.
proffile	Hypocenter output in Xprofile format.	Used by sphrel3d90 only.

Finally, several scalars are output into the Generic Output file that describe the data fit:

residual norm	Sum of the squared residuals for the entire data set
standard variance	The variance of the entire data set
original variance	Like standard, only not demeaned
residual variance	Estimate of the variance left over
mod residual norm	Sum of the squared residuals for the reduced data set
mod standard variance	The variance of the reduced data set
mod original variance	Like mod standard, only not demeaned
mod residual variance	Estimate of the reduced variance left over
percentage improvement	Estimated variance improvement for the entire data set
mod percentage improvement	Estimated variance improvement for the reduced data set

By "reduced data set" I mean the data left over after the hypocenter decoupling has been done.

As mentioned above, corrections to hypocenters are not calculated explicitly in SPHYPIT90 but are relocated in a separate step using the program SPHREL3D90. There are a number of reasons for doing things this way, most of them based on the notion that hypocenter parameters are the least well known and therefore should be allowed to absorb as much residual as possible. They also should be reviewed at every iteration and outlier events (those that have excessively high residuals or sensitivity to noise or model) should be removed. Bad hypocenters provide the classic example of the GIGO principle (garbage in, garbage out).

I/O to SPHREL3D90 is mercifully nearly identical to that of SPHYPIT90 and the conventions and procedures for coordinate system, ray tracing, weighting, and even variable names are pretty much the same for both routines. The substantial difference is that SPHREL3D90 forms the hypocenter matrix and residual vector and solves iteratively for adjustments to hypocenter. The final locations are written to a header file that can then be read in by SPHYPIT90.

SPHREL3D90 nominally runs for a user-specified number of iterations (ittmax) but will terminate before this limit is reached for the following reasons:

1. The corrections to the hypocenter are very small.

2. The hypocenter matrix is too close to singularity (small eigenvalues).

3. There are too few observations with non-zero weight.

4. Small changes in hypocenter increase the misfit to the observations. This is judged by stepping back from a correction that increased the misfit. If the number of steps back exceeds a limit (currently 5), then the iterations terminate.

5. A ray tracing error was encountered.

SPHREL3D90 allows the user to put restrictions on the locations to help stabilize them in the early going (when the information about location is very poor). The more useful of these are:

1. Fixing the depth. Epicenters are more stable than hypocenters, so generally it is a good idea to fix the depth for a few iterations at the outset. Use the parameter isw4 to do this.

2. Damping the inversion. Large residuals and small eigenvalues are a bad combination, so it can be a good idea in the early going to damp the hypocenter matrix. Use the parameter isw2 to do this.

3. Raise the residual truncation value. Large residuals in the early stages could be caused more by a bad location than by bad data (outliers) so we want to be more lenient with large residuals at the outset. Set the parameters trnc and trnc2 (see discussion of these guys above) to large values, and then decrease them with subsequent iterations.

Several kinds of information in the SPHREL3D90 output files are very useful in deciding which hypocenters to keep for inversion (or, in general, to believe). Most of this is summarized in the output header, the columns of which are:

NOTES:

1. The condition number is the ratio of the largest to smallest singular (eigen) values and is a meaure of a solution's sensitivity to noise. Generally I toss events where the condition number is > 100, but this is somewhat subjective.

2. The actual location quality is a subjective measure but these bits of output provide some guidance. I use a routine called GRADLOC to apply some grading criteria that applies grades 1-8. For tomography I keep only the grade 1 events.

It is quite possible that the program's dimensions will have to be changed from time to time to allow it to run bigger problems or to run on a smaller machine. The following discussion is to summarize what is important to know about dimensions in SPHYPIT90.There are several types of dimensioned variables in the program, some of which are defined by the user, and others that are intrinsic to the program. I have tried to streamline the dimensioning by putting all the relevant parameters in the include file parameter.h and adding some warning lines in the code, but I have not been systematic at checking every part of the code where a dimension could be exceeded. For this reason, I STRONGLY ADVISE YOU TO ALWAYS USE A SUBSCRIPT RANGE CHECKER WHEN RUNNING THIS SOFTWARE at least at the outset. For many fortran compilers this means compiling with the "-C" option.

Here are the contents of parameter.h:

c---MAXSTA is the maximum number of stations allowed
	parameter (MAXSTA=300, MAXSTA2=2*MAXSTA)
c---MAXOBS is the maximum number of observations per event
	parameter (MAXOBS=300)


c---MAXNXS is the maximum number of X interfaces
	parameter (MAXNXS=80)
c---MAXNXS is the maximum number of Y interfaces
	parameter (MAXNYS=80)
c---MAXLYR is the maximum number of Layers
	parameter (MAXLYR=20, MAXLYR1=MAXLYR+1, MAXLYR2=2*MAXLYR)
c---MAXNBL is the maximum number of elementary blocks
	parameter (MAXNBL=5000, MAXNBL2=2*MAXNBL)
c---MAXPRM is the maximum number of meta blocks
	parameter (MAXPRM=5000, MAXPRM2=2*MAXPRM)

c---MAXVAR is the maximum number of Variables
	parameter (MAXVAR=7600, MAXKVR=MAXVAR*(MAXVAR+1)/2)
c---MAXNKR is the maximum number of Specified Variables to Solve for
	parameter (MAXNKR=7600)


c---MAXKBL is the maximum number of blocks intercepted by all of the rays from one local event
	parameter (MAXKBL=5000)
c---MAXKBT is the maximum number of blocks intercepted by all of the rays from one tele event
	parameter (MAXKBT=5000)
c---MAXNBK is the maximum number of blocks intercepted by a single ray
	parameter (MAXNBK=100)
c---MAXWRK is the maximum number of elements for AUX array
	parameter (MAXWRK=MAXOBS)

Even with the above dimensions, the limits of the following are obscure and difficult to nail down in advance (this is why you should use the subscript range checker).

Variable	Meaning
kbl	The total number of blocks intercepted by all of the rays from one event
kblt	The total number of blocks intercepted by all of the rays from one teleseismic event
nbk	The total number of blocks intercepted by a single ray
nwrk	A dimension for aux, a "working array" in MFD, with a vague meaning. It should be similar to MAXOBS.

Most arrays appear in common blocks, which now are usually in include files (*.h).

NOTES:

1. The rmm matrix is currently set up for the problem of long-wavelength smoothing on a layer by layer basis. Thus, if the maximum number of variables in a given layer in nlvar, then the first index should be dimensioned as klvar = nlvar*(nlvar+1)/2. The second index is the number of different layer smoothing functions that are used.

2. The arguments to MFD must be changed if the number of observations is changed (or, more specifically, if the dimensions to a() and p() are changed) because a and p are treated as one-dimensional variable in the routine.

The following subroutines/functions are called by SPHYPIT90:

addcov

Adds the inverse covariance matrix to the normal equations.

baz2az

From two points on a sphere, the distance between them, and the azimuth from the first to the second, calculates the azimuth along the same ray at the second point (pi - back-azimuth). DOUBLE PRECISION; Uses internal psuedo-km coordinates with Y negative; OVERWRITES OLD AZIMUTH. Called by MOVEITSPH.

bjdaz2

Calculates distance, azimuth, and back-azimuth between two points on a sphere.

buildg2

Builds the matrix of normal equations

calctraysph

Used by TELTIMSPH to calculate the raypath of the teleseismic ray.

cppkpo

Applies the a priori model constraint.

delcal2

Calculate distance from base of model to surface (for teleseisms).

deploc

Locates an event in the appropriate layer.

direct

Calculates the direct ray path

diff

Calculates the distance deficit for direct

distsph

Determines that distance of any point from the center of coordinates

distsph2

Determines that distance of any point from the center of coordinates (spherical version).

errcorsph

Applies spherical corrections to error estimates.

elpcr

Looks up an ellipticity correction from a table.

etoh

Converts epochal time to "human" time. This and other routines in the file "timesubs.f" are fortran adaptations of the DSAP C code.

find

Called by ELPCR

flatz

Converts from radial to earth-flattened depth

flatvel

Converts from radial to earth-flattened velocity

getline

Reads a line of text from a file, skipping over blanks and comments (lines beginning with #).

getfield

Recovers variables from a line of text.

glat

Converts geographic to geocentric latitude

glatinv

Inverse of GLAT; converts geocentric to geographic latitude

inblkmod

Reads in the block model, reparameterization table, dampers, slowness model, in inverse covariance matrix (as required).

htoe

Converts "human" time to epochal time. This and other routines in the file "timesubs.f" are fortran adaptations of the DSAP C code.

interp

General linear interpolating routine.

julcon

Converts calendar month and day to Julian day.

lat2del

Calculates the intersection of a ray leaving a point at a given azimuth with a specified latitude line; i.e., calculates distance along a ray to its intersection with a E-W block boundary. DOUBLE PRECISION; Uses internal psuedo-km coordinates with Y negative. Called by TRACKSPH and RSLOWSPH. NOTE: this operation is inherently multi-valued, as a great circle will in general cross a given latitude line twice or not at all. The output value is a particular root of a quadratic equation that seems to be correctly chosen, but there have been problems with early versions of this subroutine.

lng2del

Calculates the intersection of a ray leaving a point at a given azimuth with a specified longitude line; i.e., calculates distance along a ray to its intersection with a N-S block boundary. DOUBLE PRECISION; Uses internal psuedo-km coordinates with Y negative. Called by TRACKSPH and RSLOWSPH.

mfd

Decompose the hypocenter matrix using singular value decomposition.

mfss3

Decomposes the matrix of normal equation

mlsh

A shortened form of mlss3.

mlss3

Completes matrix inversion by multiplying G^-1d

mlsh A short form of mlss3 designed to return only diagonal elements

moveitsph

Called by srtime to move ray from one block to another

neighb

Calculate resolution and covariance for neighboring blocks

newpos

Determines coordinates of a point on a sphere from those of another point, and the azimuth and distance between them. DOUBLE PRECISION; Uses internal psuedo-km coordinates with Y negative. Called by MOVEITSPH.

offbig

Determine which off-diagonal blocks have the most resolution power.

openold

Utility routine for opening existing files.

outputsph

General output routine - used for creating output maps of various arrays.

parsedata

General arrival time input routine.

parsedesc

Parses the data description file.

parsefile

Parses the file containing all the file names for attachments.

parseparms

Parses the file containing all the job control parameters

parseparms

Parses the file containing the station locations.

parvel1d

Parses the file containing the 1D model.

readrayp

Reads in the IASP ray parameter file.

readtbl2

Reads in the IASP travel time files.

readtbl2

Reads in the IASP 1D velocity model.

redtab

Reads in the Dziewonski. and Gilbert ellipticity table from file elpcor

reorderp

Subroutine called by inblkmod to sort out the reparameterization.

revdaz

Finds lat and lon given distance and azimuth from a reference point. Does the inverse of bjdaz2.

rotback

Converts back from internal coordinates to global coordinates by inverse of TRCO matrix. Also, undoes earth-flattening for depth.

rslowsph

Determines an approximate 1D structure from the 3D structure

rtable2

Don't remember at the moment. Will look up later.

setdefs

Assign default values to major parameters and files.

setorgsph

Set up short distance conversion factors

slocal2

Look up slowness in model and interpolate as necessary.

srtimesph

Traces a ray through an arbitrary 1D structure

srtimesph

Traces a ray through an arbitrary 1D structure

tdirectsph

Vesion of direct used for tracking teleseismic ray paths. Serves same purpose as direct for long distances.

tracksph

Calculates time and partial derivatives in the 3D structure

tterm

Calculates the time terms for the 1D model

tterms

Calculates the time terms for the approximate 1D model

ttable2

Reads in IASP travel time tables.

turncheck

Checks for funny rays that cross E-W block bounds that have Y coordinates that are not between source and receiver, and finds the min/max Y value 'turning points' of these rays. Called by MOVEITSPH and RSLOWSPH.

track2

For rays that were flagged by TURNCHECK, this routine breaks the ray into two segments at the 'turning point' and calls TRACKSPH separately for each segment.

velder

Calculates 1D velocity derivatives as required.

u2

Determines the hypocenter matrix annihilator.

1. The current version of the software is a beta version, meaning that not all of the features have been fully tested since the last "upgrade". As most of the changes have been in I/O and are not procedural, I suspect that there might be some simple bugs, but only time (and trial) will tell.

2. The routines that take care of the inverse covariance matrix business have not been tested in quite a while and are almost certainly not compatible with reparameterization. Because reparameterization enforces a desired scale length on the model it is not clear that you would want to do this anyway. There is no conceptual barrier to making the two compatible, but it certainly has not been tested.

ACKNOWLEDGEMENTS

Many persons contributed to this program over the process of its (ongoing) development. Principal among these were Geoff Abers, Bill Ellsworth, Craig Jones, Paul Friberg, and Kaye Shedlock. If this program works for you, let me know. If it doesn't, talk to the other guys.

Go to Table of Contents