.. _sac2:

*********************************
Seismic Analysis Code 2
*********************************

This lab focuses on additional material relating to graphics in SAC, and how to write SAC macros (which are essentially scripts to carry out a series of SAC operations, allowing you to conduct more complex computations using SAC).

===========================
Details about SAC Plots
===========================

---------------------------
Quick and Dirty Plots
---------------------------

Run SAC on your computer, and load one of the seismograms that you downloaded for the previous class, and plot the function. Your plot should look something like this:

.. image:: sac2_quickdirty.png

Notice the box in the lower right hand corner-- for my plot, it says "26" but depending on the length of the signal that you plotted you will see a different number. This number signifies that SAC is only plotting every 26th point in this graph.

Why does SAC do this? On the original Tektronix system, the data transfer rate is painfully slow, to the point that displaying a plot of seismic data at a typical sampling frequency would take several hours. Thus, by default SAC desamples long time signals before plotting them in order to reduce the time needed to transfer the data to the display.

Unfortunately, this is not implemented very well. When SAC performs this operation, it simply desamples the data, taking very few points. This is unlikely to give accurate information about the amplitude of the signal (because the probability that the maximum value of a particular wiggle is displayed is only 1/26), and means that the plots are not particularly useful for analysis (particularly for very long signals with many data points, which would require even more aggressive desampling). A better way to do this would be to take the maximum value within every desampling window and plot that value, making the signal amplitudes useful (referred to as decimation). While it requires an additional processing step, the time required to do this extra processing step is minimal compared to the data transfer time on the original Tektronix display.

Fortunately for us, since we have modern computers with modern displays, we do not have to worry about long times to display signals. To prevent this behavior in SAC, we need to turn off the "quick and dirty plot" option (which is set to "on" by default. If you type ``qdp off`` into the command line and then plot your data again, you will see a plot with no number in the corner:

.. image:: sac2_fullplot.png

This plot has information that accurately reflects the signal amplitude (and notice that the display time was not noticeably different than before), though for us the amplitudes aren't markedly different. If you are going to use SAC plots for any real scientific purpose, turn off the quick and dirty plots!

--------------------------
Labeling Axes and Titles
--------------------------

To label axes in SAC, use the ``XLABEL`` or ``YLABEL`` command. The syntax is ``xlabel on '<label>'`` for the :math:`{x}`-axis (similar for the :math:`{y}`-axis. For example, if I wanted my :math:`{x}`-axis to be labeled as 'Time (s)' I would type ``xlabel on 'Time (s)'`` into the terminal. As you may remember, this will have no effect on the current plot until you issue the ``plot`` command again, so you need to issue these commands before plotting for them to take effect.

To add a title, use the ``TITLE`` command. It is issued in the same way as the ``xlabel`` command, so to title a plot, type ``title on 'Seismogram'`` to add the title 'Seismogram' to your plot. Again, you will not see any effect from the ``title`` command until you issue the plot command.

---------------------------------
Changing the Plot Appearance
---------------------------------

We have already discussed how to change the limits on a plot and how to make a color plot, but there are additional customizations that are possible. To change the spatial size of a plot, use the ``XVPORT`` and ``YVPORT`` commands. The syntax is ``xvport <min> <max>``, where ``<min>`` and ``<max>`` are numbers between zero and one that specify the minimum and maximum values of the axes as a fraction of the entire plot window. The default values of ``xvport`` are 0.1 and 0.9, and the default values of ``yvport`` are 0.15 and 0.9. To make a smaller plot, I would issue the commands ``xvport 0.1 0.5`` and ``yvport 0.15 0.55`` and then the plot command.

To put multiple plots in the same figure, we use the ``xvport`` and ``yvport`` commands in conjunction with the ``BEGINFRAME`` command (shortened to ``begfr``). When you enter the ``beginframe`` command, SAC will add a new axes to the existing figure when subsequent plot commands are issued, rather than making a new figure. Thus, once we enter the ``beginframe`` command, we can add plots in specific locations by setting ``xvport`` and ``yvport`` and then giving the plot commands. When the ``ENDFRAME`` command is issued, the default behavior resumes. An example to make a plot with two axes is: ::

    beginframe
    xvport 0.1 0.9
    yvport 0.15 0.475
    funcgen impulse
    title on 'Impulse'
    xlabel on 'Time (s)'
    plot
    xvport 0.1 0.9
    yvport 0.575 0.9
    funcgen boxcar
    title on 'Boxcar'
    xlabel on 'Time (s)'
    plot
    endframe

.. image:: sac2_subplot.png

----------------
Saving Plots
----------------

So far, we have been making plots and displaying them on the screen. What if we want to save a plot for use in another document?

The simplest way to save figures from SAC is to use ``SAVEIMG`` to directly output a figure window to a standard graphics format. This command is flexible and will let you save the currently displayed figure window to PDF, PS, or PNG formats without doing any additional work. Once you have a figure up and displayed, simply enter ``saveimg <filename>.<format>``, where ``<filename>`` is your desired name for the file and ``<format>`` denotes the extension of the graphics format that you wish to use. Thus for example, to plot an impulse function and save as PDF, you would do the following: ::

    funcgen impulse
    title on 'Impulse'
    xlabel on 'Time (s)'
    plot
    saveimg impulse.pdf
    
This has the advantage over the more complicated procedure described below in that the labels are rendered as text rather than as lines. Because of this, I would generally advocate using ``saveimg`` to produce plots.

However, if you are interested in more details about producing graphics in SAC, there is another procedure that gives you a bit more control. SAC has its own graphics format, known as ".sgf" (SAC graphics format). This format is analogous to the MATLAB ".fig" files as they only readable by SAC, and is not useful if we want to include the figure in other documents like presentations, posters, or papers. To save a plot in the ".sgf" format, we need to follow a few steps.

SAC has two different output devices. A device is just a fancy name for the destination where the graphics will be displayed. So far we have used the screen as our output device, but we have the ability to tell SAC to send our plots to a different device. To use the screen as an output, we use the ``XWINDOWS`` device (can be abbreviated ``x``). You may have noticed so far that when we issue a plot command, the plot window shows up in an application called "X11" which is the xwindows application for a Mac.

The other available output device is ``SGF``, or **SAC Graphics Format**. If we choose this output device, plots will not be displayed on the screen, but instead will be written to file. To change the output device, we use the ``BEGINDEVICES`` command (``bd`` in shorthand). By default, the device is set to ``xwindows``, but by typing ``begindevices sgf``, all future plots will be saved to file as an ".sgf" file instead of being displayed on the screen. To switch the device back to the screen, type ``begindevices xwindows``.

When you use the sgf device, each plot command that you execute will save to a new sgf file. The first plot after you start SAC will be saved to "f001.sgf" in your current working directory. Subsequent commands will save as a file with the name "f002.sgf," "f003.sgf," etc. If any file already exists (say from a previous SAC session), then it will be overwritten without a warning, so if you want to save your sgf files, you will need to rename them after each SAC session.

To convert an sgf file into a graphics format that is useful outside of SAC, we use the ``sgftops`` utility (note that I called ``sgftops`` a utility, as it is not an official command in the same sense that ``plot`` and ``read`` are commands). This converts an sgf file to PostScript, which is a standard graphics format that can be opened with Applications such as Preview on the Mac (when Preview opens a PostScript file, it automatically converts it to PDF), Adobe Illustrator, or on the command line using ghostscript. To convert a file, use the syntax ``sgftops f001.sgf sacplot.ps`` which converts the sgf file "f001.sgf" to the file "sacplot.ps" which can then be viewed in another application.

There are two additional options when using ``sgftops`` that you can include in the command. The first specifies the line thickness in the resulting PostScript file (the default value is 1), and the second allows you to specify the page layout or add an id that specifies some details about the plot. To add an id, make the second option ``i``, to change the page layout specify ``s``, or ``si`` for both. The ``s`` option is useful if you want your plot in portrait rather than the default landscape mode -- this option will prompt you to enter specific values for the :math:`{x}` and :math:`{y}` translations (relative to the lower left corner of the page, assumed to be US Letter size in portrait orientation), rotation (measured counterclockwise from horizontal), and scale. An example, to obtain a plot with thicker lines type ``sgftops f001.sgf sacplot.ps 1.5`` into the command line. To rotate and scale the plot into portrait mode, use ``sgftops f001.sgf sacplot.ps 1 s`` and then at the prompt, type the following four numbers, each followed by a carriage return: 0.5 0.5 0 0.75. This will result in a plot in portrait mode.

As you can see, this requires considerably more effort than using ``saveimg`` without a huge difference in the final outcome. Because of this, I do not use the SGF procedure much any more, and simply use ``saveimg`` to save figures in SAC.

**Exercises:** Try out some of the plot commands described above. First, plot a seismogram with the quick and dirty plot option off. Then, plot a seimsogram with labeled axes and a title, save it as a PDF file. Compare with the results using an SGF file converted to PostScript. Finally, use ``beginframe`` along with ``xvport`` and ``yvport`` to illustrate four of the different functions created by ``funcgen`` in a 2x2 grid of plots.

==================
SAC Macros
==================

Like with MATLAB .m files, we can write scripts to carry out a series of commands, called "macros" in SAC. SAC macros are text files containing a series of commands that can be called either from within SAC or directly from the command Terminal.

For example, here is a series of commands. Enter them into a text file (using an editor of your choice), and then launch SAC. You can name the file whatever you like, but I prefer to include ``.macro`` at the end so that I know that the file is a SAC macro. From inside SAC, type ``macro test.macro`` (using the name that you chose for your file) to execute the macro. For example, ::

    # SAC macro to generate three functions, write to disk, and plot
    funcgen impulse
    write impulse.sac
    funcgen boxcar
    write boxcar.sac
    funcgen triangle
    write triangle.sac
    read impulse.sac boxcar.sac triangle.sac
    plot1

``#`` and ``*`` serve as comment characters in SAC, so that line will not be executed. However, you cannot put inline comments into SAC, so comments must be on their own line. Alternatively, you can run the macro from the Terminal by typing ``sac test.macro`` directly from the command line.

--------------------------
Input Arguments
--------------------------

While this type of macro is easier than typing a series of commands into the command line (where you are likely to make a mistake if you are not careful), macros are more useful when they can take input arguments like a function in MATLAB or Python. SAC macros do not have an explicit way to write functions, but you can pass an arbitrary number of arguments to any SAC macro, and the parameters will be available using the variables ``$1``, ``$2``, etc. As a simple example, if we wanted to write a macro to read three SAC files and plot them, change your macro to the following: ::

    # SAC macro to read three files from disk and plot
    read $1 $2 $3
    plot1

If you type ``macro test.macro impulse.sac boxcar.sac triangle.sac`` into the command line. You should get the same plot as before. If you call the macro without specifying any arguments, you will get a prompt asking you to input the arguments.

What if you want to read an arbitrary number of SAC files? You can make the macro more flexible by specifying keywords, which can take an arbitrary number of arguments. Change your macro to the following: ::

    # SAC macro to read an aribtary number of files from disk and plot
    # files to read are set using the "files" keyword
    $keys files
    read $files
    plot1

By typing ``macro test.macro files impulse.sac boxcar.sac`` into the command line, you can run this macro, and you can change the number of input files to whatever you like without changing the macro. If you do not give any values for ``files``, SAC will prompt you to enter the values (Note: I have found that this does not always work on SAC installed in the Mac Lab, so you may end up getting an error if you do not provide the input arguments).

You can also specify default parameters. If you want your macro to plot impulse.sac, boxcar.sac, and triangle.sac by default, then change the macro to: ::

    # SAC macro to read files from disk and plot
    # files to be read are set using the "files" keyword
    # if no files provided, default files are impulse.sac boxcar.sac triangle.sac
    $keys files
    $default files impulse.sac boxcar.sac triangle.sac
    read $files
    plot1

Now you can run your macro with any number of arguments, or with no arguments to plot those three functions by default. (Note that I have also found that setting default parameter values does not always work correctly in the Mac Lab for some reason, so you may sometimes get an error when doing this in the Mac Lab.)

--------------------
Variables
--------------------

One of the valuable aspects of any programming language is the ability to assign values to variables. This can be done in SAC using what is known as a "blackboard" variable. You can set a blackboard variable using the ``SETBB`` command, followed by the variable name and the value. For example, ``setbb xmin 40`` will define a variable ``xmin`` with the value 40. To access the variable, use ``%xmin``.

You can also calculate other parameters from variables, which you can do using the ``EVALUATE`` command. Let's say that we want a macro that will plot an arbitrary number of functions from an arbitrary starting time for 100 seconds. We will have keywords for the entries for xmin and files, or ::

    # SAC macro to read files from disk and plot with a time range of 100 s
    # files to be read are set using the "files" keyword
    # if no files provided, default files are impulse.sac boxcar.sac triangle.sac
    # also must provide argument xmin to set start time on plot
    $keys xmin files
    $default xmin 0
    $default files impulse.sac boxcar.sac triangle.sac
    setbb range 100
    evaluate to xmax $xmin+%range
    read $files
    xlim $xmin %xmax
    plot1

Run this by typing ``macro testmacro.txt xmin 25 files impulse.sac boxcar.sac triangle.sac`` to plot the three functions from 25 to 125. The first line defines that input arguments for ``xmin`` and ``files`` should be provided, then it determines ``xmax`` by adding 100 to ``xmin``. Then, the macro reads data from the files, sets the :math:`{x}` limits, and plots the signals.

Sometimes, you might want to use an input argument or a blackboard variable to represent a part of a file name. A common example might be if you are reading data froma a seismic network, and the file names are identical with the exception of one piece representing the station name or the component in question. To use a variable within another string, you need to append the variable name with the same character that is used to access its value. So for example, let's say you have two stations, S1 and S2, with components N, E, and Z, and the filename convention from the network operator is stationS1N.sac, stationS1E.sac, stationS2Z.sac, etc. The following macro would correctly read the file specified by the ``stat`` and ``comp`` inputs: ::

    # SAC macro to read a station component from disk and plot
    # must input the station name and component when calling
    $keys stat comp
    read station$stat$$comp$.sac
    plot
    
This can be called using ``macro test.macro stat S1 comp N``. Similarly, if ``comp`` was stored as a blackboard variable instead of a keyword, the equivalent read command would be ``read station$stat$%comp%.sac``. This is very handy when complicated, but predictable, filenames are used that only differ in a couple of arguments.

------------------------
If statements
------------------------

Macros can use if statements, with a syntax very similar to that of MATLAB. Let's say that in the above macro, you only want to change the :math:`{x}` limits if ``xmin`` is greater than zero. In that case, we could use an if statement: ::

    # SAC macro to read files from disk and plot with a time range of 100 s
    # files to be read are set using the "files" keyword
    # if no files provided, default files are impulse.sac boxcar.sac triangle.sac
    # can also provide argument xmin to set start time on plot (default 0), only works if positive
    $keys xmin files
    $default xmin 0
    $default files impulse.sac boxcar.sac triangle.sac
    setbb range 100
    if $xmin gt 0
        evaluate to xmax $xmin+%range
        xlim $xmin %xmax
    else
        xlim 0 %range
    endif
    read $files
    plot1

This will set the limits differently if xmin is negative (``gt`` stands for greater than). Other comparisons available inlcude ``lt`` (less than), ``le`` (less than or equal to), ``ge`` (greater than or equal to), ``ne`` (not equal), and ``eq`` (equal). The ``else`` (or ``elseif``) is optional.

------------------------
Loops
------------------------

There are several ways to execute a loop in SAC. The first option is to loop over numeric values (like in MATLAB or using ``range`` in Python). Let's say you have a series of 10 .sac files, with names "signalx.sac", where x is a number from 0-9. To read all 10 files in using a loop, we can do so as follows: ::

    # SAC Macro to read in 10 files and plot
    do i = 0, 9
        read more signal$i$.sac
    enddo
    plot1

This will read in all 10 files and plot them all on the same axes. Optionally, you can add an increment using ``do i = 2, 8, 2`` to increment by 2 each time. Another way to write this is to use the syntax ``do i from 2 to 8 by 2``.

If you want to loop over a specific set of items, the syntax is a bit different. If we wanted to read all three components of a particular station and plot them, we might modify the above example to be the following: ::

    # SAC macro to read 3 components of a station from disk and plot
    # must input the station name when calling
    $keys stat
    do comp list N E Z
        read more station$stat$$comp$.sac
    enddo
    plot1
    
This version explicitly loops over all three components specified by whatever follows list.

If you want to get some files from the shell using a particular pattern, you can loop over files matching a certain pattern, which is very useful for reading in files. In particular, if you wanted to load the vertical component of any station in the current directory, you can do ::

    # SAC macro to plot vertical components in the current directory
    do file wild station*Z.sac
        read more $file
    enddo
    plot1

You can also specify a directory other than the current one if you use ``do <variable> wild dir <directory> expr``.

Finally, you can do a standard ``while`` loop in SAC. The syntax is ``while <condition> ... enddo``, and ``break`` will work as we have seen previously in Python and MATLAB to exit loops.

=====================
Summary
=====================

While programming in SAC is a bit less refined than programming in Python or MATLAB, it still lets you get the job done fairly well if your work requires processing .sac data files.

**Exercise:** Write a SAC macro that generates three sine functions with the same frequency and a fixed phase difference between them. Plot all three functions on the same axes, using a different color for each function. The frequency and the phase offset should be specified through command line variables, as should the number of points and the time spacing used in generating the signal. Default values should be 100 time points, time spacing of 0.05 s, and a frequency of 1 Hz (there should not be a default for the phase offset).