One of the most widely used analysis packages for regional and teleseismic seismic data is SAC (the Seismic Analysis Code). It was developed in the 1980s by nuclear test monitoring agencies in the United States, who freely made the source code and paper documentation available to academic users. From this initial distribution, the analysis package became popular in academic circles due to its ease of use and suitability for research data analysis in seismology and geophysics.

SAC's documentation was on ring-bound paper shipped along with the nine-track half-inch source code tapes that you received in the post. The academics in receipt of them generally made copies from the master document and distributed them to colleagues and students. They also tutored new users on the use of SAC, guiding them through their first session and then left them to the documentation. Consequently, much of the knowledge of SAC was passed tutorially from an experienced user, not unlike trade apprenticeship. Those intrepid enough to find and read the documentation usually exceeded their tutors' ability. The far fewer who delved into the source code learned of undocumented features of great utility. The usual reasons for this puzzling knowledge gap apply: in software development, documentation always lags feature development and, for SAC, new releases were sporadic and focused on new capability.

The original SAC documentation consisted of: (1) tutorial guide; (2) command table; (3) detailed command descriptions; (4) SAC file structure internals; (5) auxiliary program guide for programs to turn graphics to hard-copy form.

PHILOSOPHY AND STRUCTURE

SAC file format

A seismic data trace is a set of data points that is continuous in time but that may not have been sampled at an even rate. SAC's simple approach to seismic data deems that there is one seismic trace per file. Each file contains a header that describes the trace (also known as metadata) and a section that contains the actual data. The header occupies a fixed-length position at the beginning of each file, followed by the variable-length data section.

Header data is of mixed type: integer and real values, logical (true/false) values, categorical values (distinct properties like explosive source, nuclear source, earthquake source), or text (station code, event identification, wave arrival type). The seismic data in a trace is a sequence of real-valued numbers representing the sampled physical property.

Alphanumeric and binary forms

There are both binary and alphanumeric (character) formats for a SAC file. The binary version is a more compact format that is efficiently read and written, while the alphanumeric format is easier for user programs to read and write.

The alphanumeric data format is intended for ease of reading and writing and for transfer between different machine types. In this format, the header data is organized into sections based on the type of data, with each section subdivided into lines. While not an intuitive organization, this makes it easier to read and write the header data because all of the data on each line of the file is of the same type: integer, real, etc.

SPECTRAL ESTIMATION

Spectral estimation is the task of taking a time series and decomposing it into its component frequencies. The total length of the time series and, if it has been sampled at a fixed time interval, the inter-sample spacing control the minimum and maximum frequencies that the time series contains. Different methodologies may be used to estimate the power at each frequency at constant intervals between these bounds. This so-called power spectrum provides a way to quantitatively characterize the frequency content of the time series. The information is usually presented in the form of a graph of power versus frequency.

The power spectrum informs further processing avenues for the time series. Usually this involves discrimination of any signal in the spectrum from noise. If the signal is of high quality, the signal's power will dominate the noise. Thus the power spectrum will be peaked in a frequency band containing the signal. If the goal is to design a filtering strategy to minimize the noise, this analysis will suggest the type of filter and the corner frequency to use. Another application might be to seek tidal resonances at a coastal site based on repeated sea level measurements, or a marigram. In this case, the frequencies of the spectral peaks are the desired information, and perhaps their widths or positional uncertainties.

Spectral estimation is technically complex due to the characteristics of the signal under study.

Despite the broad range of utility that SAC provides, at times it may be necessary to use external applications to augment its capabilities. In addition, it is sometimes desirable to access the functionality of SAC without interacting manually with the program, for example to include it in a longer processing workflow. In this chapter we will describe techniques and give examples of how to achieve both these ends.

Note that details of the languages and applications for which examples are shown are beyond the scope of this book. The reader should seek more specialized books for that material.

AUTOMATING SAC EXECUTION

Running SAC from the shell

Executing SAC

While a decent amount of batch processing is possible within SAC using its built-in macro language (see Chapter 5), it is often useful (for example, to access functionality built into the operating system) to run SAC using scripting languages provided by the shell (under UNIX-like environments, examples include bash and tcsh).

One way is by using startup files (see Section 4.10). After the commands in it are executed, SAC then enters interactive mode. Because (usually) there is no need for an interactive phase when scripting, it is helpful to make the last line a QUIT command to terminate SAC and allow control to return to the shell.

PLOT ANNOTATION

Seismograms

For SAC, the basic display element is a seismogram plotted with time along the horizontal axis and the sampled quantity along the vertical axis. SAC's basic plotting commands, PLOT, PLOT1 and PLOT2, supply axis tick marks numbered to show the relevant scales and connect data points with lines. SAC draws tick marks around all four sides of the plot to facilitate quantitative use of the data in the plot.

Tick marks, axes and borders

SAC uses a good default algorithm that chooses tick mark intervals and numbering based on the range of values to be plotted. If the result is unsatisfactory, tick mark intervals may be explicitly set using the XDIV and YDIV commands for the horizontal and vertical axes, respectively. The number or spacing of the tick marks may be specified using the keywords INCREMENT or NUMBER. To reset the default settings, use the NICE option.

The TICKS command selectively turns off the tick marks on each of plot's four sides using the keywords TOP, BOTTOM, LEFT and RIGHT, respectively. The omission of tick marks results in clearer plots, but they lose their quantitative clarity. TICKS is mostly used to make publication-quality plots for inclusion in a talk or a publication.

Grid lines across the plot (like graph paper) may be useful in some circumstances but are distracting due to the resulting clutter. To draw grid lines on plots, use the GRID, XGRID and YGRID commands; they are off by default.

COMMAND STYLE

SAC commands are typed from the command line or read from a file. After each command is processed, SAC reads another command from its input source until it is told to stop or the input is exhausted.

Commands are single verbs (e.g., READ, WRITE), or a compound phrase (e.g., FILTER-DESIGN). Abbreviations exist for the longer or commonly used command names. A series of options that control the command's actions follow the command name. Command names or options may be typed in upper or lower case. However, when file names appear in commands, case does matter, and SAC preserves it.

White space separates options and the command name, and can even precede the command name. This is useful for indenting groups of commands for documentation purposes.

Multiple commands may be placed on the same line separated by the “;” (semicolon) character and will be processed left-to-right as they appear on the command line. Any command whose first character is * is a comment and is ignored. Thus the string ; * introduces a comment in the command listings that follow.

SAC supplies default command options if they are not specified. Command options, once set, stay in force for future uses of the same command. This provides a way to tailor personal command defaults. SAC can read a file of commands setting your personal defaults before it reads the input. They will be described in detail in Section 4.10.

WHAT IS SAC?

SAC is an acronym for the Seismic Analysis Code, a command line tool for basic operations on time series data, especially seismic data. SAC includes a graphical interface for viewing and picking waveforms. It defines a standard type of file for storing and retrieving time series and also reads files written in other data formats (SEG-Y, MSEED, GCF) used during field collection of seismic data. SAC is self-contained and does not rely on network access for any of its capabilities, including documentation, which makes it useful for field data quality control.

SAC reads data formats (CSS and GSE formats) used by nuclear test monitoring agencies. It also contains programming language constructs that provide basic methods for developing elaborate, multi-step analysis methodologies. Collectively, these features make SAC a useful interactive platform upon which customized analytical methods may be built and prototypical procedures may be developed.

SAC is widely known. The IRIS Data Management Center (DMC), one of the largest whole-Earth seismological data repositories in existence, allows data to be requested in SAC form. The instrument response information provided by the DMC's SEED reading program, rseed, is usable by SAC in pole-zero or evalresp form. Owing to SAC's longevity, a rather large and well debugged software tool ecosystem has evolved around its file format. One such tool, jweed, searches for and retrieves data held by the DMC. SAC data and SAC-compatible instrument response information are among its output options.

THE CONCEPT OF 3D DATA

SAC has facilities for handling a basic type of three-dimensional (3D) data, which are function values evaluated on a regular grid of (X,Y) positions. This type of data can be viewed as evenly sampled in space and fits naturally into SAC file types. The file consists of the string of samples at each grid point, and file header information that flags it as 3D gives the X and Y dimensions of the grid. SAC calls data of this type XYZ.

The simplest visualization of 3D data is map elevation data. Imagine a section of land whose elevation is specified on a grid every 100 m eastwards and 100 m northwards. Over a plot of land 1 km by 1 km, the elevation could be described by a string of 121 values, the first 11 values along the westernmost traverse and the final 11 along the easternmost traverse.

Another type of data that can be represented in 3D is a spectrum as a function of time. If equal lengths of time are cut from a seismogram, the number of points in their FFTs will be the same. If the interval between spectral samples is constant, these values will also form agrid.

A final type of data inherently 3D is a misfit function evaluated over a two-dimensional (2D) grid of points. The misfit minimum within the sampled parameter space yields the optimum combination of parameters that best fit some observations, and the shape of the minimum provides a measure of the joint uncertainty of the parameter estimates.

Of the many data types that are used in programming – numbers, strings, symbols, lists – numbers are perhaps the most familiar. You can work with all kinds of numbers in Mathematica, but what distinguishes it from traditional programming languages and other computational systems is that with it you can operate on numbers of any size and to any degree of precision. In this chapter we will explore some of the issues related to working with numerical quantities and show how you can incorporate these ideas into programs that involve numerical computations to gain greater control over the precision and accuracy of your results as well as to improve the efficiency of your numerical computations and programs.

Numbers in Mathematica

One of the first things you will notice as you start using Mathematica is the manner in which it treats numbers compared with other systems such as calculators, traditional programming languages, and other technical computing software. In most traditional programming languages, you must declare the type of number your functions can take as an argument. Although Mathematica automatically handles such details for you an understanding of the different number types and how they invoke different algorithms is helpful for taking full advantage of Mathematica's numerical capabilities and for writing efficient programs.

Mathematica operates differently depending upon the type of input you give it. For example, the following two inputs each compute sin(π / 3) but something quite different results.