Introduction

The distributed codes have two purpose. First to perform the tomography and second to access and display the tomography results. The compiled programs are then invoked from within tailored shell scripts. Since the organization of the distribution is discussed in another link, the discussion here is on the compile codes.

src/

The contents of the src directory are

alltogc.c    compTimesSpk.c  General	geoggeoccnvc.c	invsubs.h  Makefile	 nmulti_lsqr.c

compTimes.c  doquad.f	     geoconv.f	geoggeoccnvc.h	libget.h   Makefile.gcc  ctou96.f

and the contents of src/General are

cfft.c		 GetMisfit.c	      in_rect.c       Makefile	    nreadobs.c	 wtminc_lsqr.c
delta_az.c	 getpar.c	      invsubs.h       Makefile.cc   ran2.c	 
gcarc.c		 get_perturbation.c   latlon_to_xy.c  Makefile.gcc  readobs.c
get2Dsmoother.c  get_random_number.c  libget.h	      minc_lsqr.c   row_io.c
get_fatray.c	 get_rough.c	      libtomo.a       mnrutil.c     test_lsqr.f
get_lengths.c	 image_smooth.c       LocalSearch.c   normalize.c   write_model.c

These codes were originally written by C. J. Ammon. Slight modifications to the original code were made for this distribution.

Compile

The Makefile is set up to use the gcc/gfortran compilers. To compile the codes, perform the following steps:

cd src
cd General
rm -f libtomo.a       [This ensures that the proper library is used. An old one may be from another compiler/system]
make clean
make all
cd ..                 [move up to the parent src directory]
make all

As a result, the following executables are placed in the bin directory, one level above:

alltogc
doquad
compTimesSpk
compTimes
nmulti_lsqr

The use of these codes are as follow:

alltogc

This program reads the file of observations and converts the coordinates from a geodetic (geographic) system to a geocentric system. This is because the inversion code assumes that the Earth is a sphere. This is a simple program that reads a line, changes the coordinates, and then writes a line in the same format.

The input line consists of 10 columns:

Column	Value
1	Identifying string
2	Geodetic latitude of starting point
3	Geodetic longitude of starting point
4	Geodetic latitude of end point
5	Geodetic longitude of end point
6	Dispersion velocity for the path
7	Period
8	Error in velocity
9	L or R to represent Love or Rayleigh data
10	Distance of path in km   - this is used together with the velocity to get the travel time for the path

The output is the same format, except that the latitudes geocentric. The following displays a sample set of input and output for this program.

Input:
MALL.DSP   49.4504  -82.5079   47.9154  -89.1514  4.0983 27.000000 0.001000 L 517.910000
MALL.DSP   49.4504  -82.5079   47.9154  -89.1514  4.0712 26.000000 0.001000 L 517.910000
MALL.DSP   49.4504  -82.5079   47.9154  -89.1514  4.0470 25.000000 0.001000 L 517.910000


Output:
MALL.DSP   49.2602  -82.5079   47.7239  -89.1514  4.0871 27.000000 0.001000 L 516.492737
MALL.DSP   49.2602  -82.5079   47.7239  -89.1514  4.0601 26.000000 0.001000 L 516.492737
MALL.DSP   49.2602  -82.5079   47.7239  -89.1514  4.0359 25.000000 0.001000 L 516.492737

doquad

The output of the tomography code is a large file that consists of lines such as

45.922 -94.005 3.589 633.94
45.922 -93.682 3.645 1057.00
45.922 -93.359 3.677 510.06

which are the geocentric latitude and longitude of each cell, the dispersion velocity and the sum of all ray paths through that cell in km. The program doquad is controlled from the command line arguments:

 Usage: doquad -LAT lat -LON lon  -F filename -h [-3]
  Interp[olate grid results to give dispersion for LAT,LON
  -LAT lat  (no default) geographic latitude
  -LON lon  (no default) geographic longitude
  -3        (no default) Normally read lat,lon,v,distance
                         from tables. With -3 read lat,lon,v
  -C        (no default) Debug: output results for empty cells

In normal operation, converts the geodetic latitude/longitude on the command line to geocentric coordinates, and then it searches through the output of the tomography code to define the 4 cells surrounding the desired latitude and longitude. It then interpolates to output the estimated velocity at the desired coordinate. Internally in the source code a variable dl=50 which means that the interpolation should only use those cells for which the sum of all ray paths through the neighboring cells is greater than 50km. This value is acceptable since the tomography grids used are 25 km x 25 km. The -3 flag does not use that raypath sum.

compTimes

compTimesSpk

These programs use the observed data set to make a synthetic data set for a visual resolution analysis. The compTimes creates a checkerboard data set, which the compTimesSpk alternates +- spikes on an 8x8 grid. Because of the data sets and the smoothing, the output of the tomography will spread the spike over neighboring cells, and that spread is a rough estimate of the spatial resolution.

nmulti_lsqr

This is the inversion program. It performed the inversion by an iterative process because of the size of the unknowns. thus a formal error analysis is not possible. The use of the code is described somewhat in the processing scripts. This code has the ability to invert many periods at the same time. However in this distribution, the program is invoked separately for each period in the data set. The results are similar, but the computation time is much faster.

ctou96

If tomographic inversions are performs for group velocity and phase velocity, the dispersion at a given coordinate must be related because of the definition of the group velocity in terms of the phase velocity. This program attempts to test this statement given the results of the independent inversions. Since we have found that the group velocity measurements are not as robust as those of the phase velocity, this code provides a way to weight the group velocity values against the phase velocity for structure inversion.

The relation between group velocity, U, and phase velocity, c, as a function of frequency f is

Now assume that we have modeled the 1/c with a set of cubic Bézier curves as a function of ln f as

from which one obtains

and hence

These are in the form to use a least squares technique to estimate the P_i coefficients.

the only choice to be made is the number of breakpoints to use. The minimum is 4 for a cubic spline and the maximum permitted by the code is 20.

In operation select either the Love or Rayleigh wave observations from the observed dispersion, determines the unique periods of the data set, and only use the group velocities within the same period range as the phase velocities. The data are then spline fit and a residual is computed. The final output will be for the number of break points giving the minimal residual.

Caveat: As with any curve fitting, a higher order will give a better fit, but there may be many excursions in the predicted dispersion because a higher order polynomial is being fit.

The program ctou96 is controlled from the command line arguments:

 Usage: ctou -DI dispin -DI dispout -NBMAX nbreak  -L -R
  -DI dispin  (required) Input dispersion file in surf96 format
  -DO dispout (required) Output dispersion file in surf96 format
                         containing least squares smooth phase 
                         and group velocity curves
                         These are at selected periods between the extremes in the data set and not at all observed points
  -NBMAX  nbreak (default=4) Maximum number of breakpoints for spline
  -L          (no default) Process Love-wave dispersion
  -R          (no default) Process Rayleigh-wave dispersion
    NOTE IT IS ASSUMED THAT ONLY THE FUNDAMENTAL MODE IS PRESENT
    NOTE The group slowness is included in the inversion, but
    NOTE is weighted 0.01 of phase slowness since the question
    NOTE is whether the group velocity is compatible with the phase velocity
 -?                   Display this usage message
 -h                   Display this usage message

An example of its use is given in CUSExamples/index.html with the script FDODISPspline.