## Enhancements Bug and Fixes 330

### Introduction

The purpose of this page is to notify users of the Computer Program in Seismology (CPS) package of repairs or extensions made to the code. Whenever this page is updated, the down loadable version of CPS is updated at the download FTP site.,

### Enhancements to Version 3.30 codes

2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014 2015 , 2016, 2017, 2018 For each year, the changes are (listed in in order from newest to oldest)

#### December 31, 2017

• VOLV/src/shwmod96.f -  removed extraneous NLAY and NBDY
• VOLV/src/refmod96.f -  correctly implemented the -HS source_depth for reflection multiples.
• VOLVI/src/hrftn96.f - cleaned usage text for -h command
• VOLVI/src/hudson96.f -  added -V flag for verbose output
• VOLVIII/src/sacpol.f - related the arrow length to the length of the axes
• VOLVIII/gsac.src/gsac_ch.c - Previously if the epicenter and station latitudes are defined, then the distance, azimuth and back azimuth are calculated. However when working with exploration data, oen may want to enter these coordinates as kilometers instead of as degrees and to manually specify the distance and not to compute these values. The behavior is now controlled by the LCALDA header variable. If LCALDA == TRUE, then it is assumed that the coordinates are in degrees and the distance and azimuths are computed for the Earth. Now if LCALDA == FALSE, then the EVLA, EVLO and STLA, STLO fields are not used to comptue the distance and azimuth. Instead the user must compute these separately and enter them with the ChangeHeader command.
• VOLVIII/gsac.src/gsac_corr.c - added the LCALDA == TRUE to the computation of distances and azimuths for the following reason. One use of cross-correlation is to obtain the empirical Green's functions from ambient noise. In this case it is desired to get the distance between the two stations. So the EVAL EVLO fields of the original waveforms are ignored, the coordinates of the master station  are used for the  EVLA and EVLO for all traces, and the distances are computed for the Earth. This is the case for LCALDA == TRUE.    Another use is to  focus on the great circle distance from the earthquake to each station. In this case we are interested in the difference in great circle distances between the two stations. This is accomplished by setting LCALDA = FALSE.

#### October 16, 2017

• Because of a request to estimate phase velocities using two two stations on the same great circle path, changes were made to do_mft, sacmft96, and gsac:
• gsac - the correlation command now correctly outputs the difference in epicentral distance between the master trace and the other traces
• sacmft96 -  supports the command line flag -P0 ( p zero ) for the computation of phase.  Given the CPS definition of the Fourier transform, the phase velocity is defined by c = omega * dist / (tphase + N 2 pi ).
• For interstation Greens functions from ambient noise, tphase is -phase + pi/4 + omega*dist/U, where phase is the instantaneous phase at the maximum of the envelope and U is the corresponding group velocity. For the derivation of this refer to http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/EMPIRICAL_GREEN/MFT.pdf
• For the cross-correlation of two stations on the same great cirle path, tphase is -phase + omega*dist/U
• do_mft - supports the command -P0 ( p zero ) to activate the interactive phase velocity module that calls the modified sacmft96

#### September 12, 2017

• Corrected logic error in the file PROGRAMS.330/VOLII/src/do_mft4.c which caused the phase match filter NOT to execute. The old code was
#ifdef DEBUGfprintf(stderr,"ModeSelect: %d pathname: %s\n",ModeSelect, pathname);sprintf(ostr,"%s -F %s -D disp.d -AUTO",pathname, fname);#endif
The corrected code  removes the ifdef/endif which caused this important section to be bypassed.
fprintf(stderr,"ModeSelect: %d pathname: %s\n",ModeSelect, pathname);sprintf(ostr,"%s -F %s -D disp.d -AUTO",pathname, fname);
Thanks to Leticia Duca, Universidad Nacional de La Plata, AModifiedrgentina

#### July 7, 2017

• Corrected spulse96 for use with slat2d96.  The description in /eqc/eqc_cps/TUTORIAL/SLAT2D/index.html     is correct. The implementation in spulse96 was not.  The problem was that there was a geometrical spreading term sqrt(wvno r)   where the wvno must be the wavenumber at the receiver, and not the average wavenumber over the path.  The documentation is correct. the example has been rerun.
• On 29 JUL 2017 corrected sregn96 slegn96 tregn96 tlegn96 to handle the case of a source in the halfspace. The code splits a layer to place the source at a layer boundary. However this was not done properly in the halfspace when the thickness was 0 or less that required to place the source depth there. The correction involved making sure the the halfspace thickness (which is of course infinite) can define a layer to the split.

#### May 27, 2017

• Modified the routines ramat.f,  jramat.f, rftndr96.f and  modls.f and the program rftndr96.f to add the command line option -SAC or -sac to be able to run the program from the command line to output both the observed and predicted receiver functions as sac files. The sac files will be of the format CCM_____2017122212232.5, e.g., SSSSSSSS_YYYYMMDDHHMMalp  where alp is the Gaussian filter parameter. The purpose is to permit the user another way to compare the observed and predicted receiver functions instead of using the graphics in the joint96 or rftn96 programs. Basically normal operation of the inversion is performed, the inversion is paused, and then one executes rftndr96 -sac from the command line.

#### May 1, 2017

• The current release is NP330.Apr-12-2017.tgz
• CRITICAL: The subroutine varsv used by  sdisp96, scomb96, sregn96 did not handle correctly the case of vertical wavenumber = 0 in the propagator matrices. The problem arose with the evaluation of sin $\mathrm{\left(\nu z \right)/}$ z.
• VOLIII:
• Subroutine wrdisp.f now works.
• sdisp96, scomb96, sregn96:
• tdisp96, tregn96, tcomb96:  The problem here was that whereas the vertical wavenumber ν   for isotropic media is either real or imaginary, it can also can be imaginary.   for some VTI media. Credit is due to Prof. S. N. Bhattacharya of New Delhi for pointing this out. This realization required a major code rewrite to handle the complex vertical wavenumbers. Prof. Bhattacharya showed that the propagator matrix and its compound version are always real, however the boundary conditions for a halfspace may introduce a complex condition. For surface waves, the period equation can be complex, but it seems is if the real part has a zero, the complex number is also a zero for a dispersion value. The surface wave eigenfunctions and dispersion values are real though.;

To start the search for dispersion curves in isotropic media, one can search between  0.9 βmin ≦ c ≦ βmax. For a model for which velocity increases with depth, the first value would be slightly less than the classical Rayleigh wave phase velocity in a halfspace, while the second limit would permit higher modes. For a general TI media, the following is done:

• For the minimum value, the topmost solid layer is used to find the zero of the Rayleigh-wave period equation  at a frequency of 1 Hz. Since this is treated as a halfspace, there is no dispersion.
• For the maximum value, we start at the highest S velocity, and then search for lower phase velocities until both the quasi-P and quasi-SV vertical wavenumbers have a positive real part. There is no requirement that the imaginary part be zero.   This must be tested more.

• sdpdsp96, sdpegn96, tdpegn96 tdpsrf96, sdpsrf96: Changed label from sec to s.
• VOLII:
• do_mft: To speed processing of a lot of waveforms, the AutoFeed feature has been extended so that much less user interaction is required - previous processing choices are remembered. The original sequence, which is still the processing default, is to pick the Sac file, then set the units, then set the processing parameters, and  interactively pick the group velocities or phase velocities, and then return to the first page. When AutoFeed is selected, then it is assumed that the units and processing parameters are to be the same, the next file is then automatically processed. This speeds the processing by not having to display three menu pages with the required interaction for meny selection.
• VOLV:
• mtinfo: in subroutine mteig, the array real*8 z(3,3) was not defined. This cause a memory problem
• VOLVI:
• tspec96: Accounted for  the problem in the ratio   sin (nu z ) / nu as nu → 0 .
• VOLVIII/src
• sacpol:  Introduced -ARROW flag to indicate direction of increasing time
• VOLVIII/gsac.src:
• Changed the following source code for gsac.
• gsac_prs.c: If the reduction velocity is used, e.g., ' p pval' then the PRS???.CTL will covert that to the proper VRED for refmod96
• gsac_trans.c:  Changed NCHAR from 200 to 300 to permit long names for files
• gsac_map.c: Change the gmtset line to read
gmtset BASEMAP_TYPE FANCY _DEGREE_FORMAT -D MEASURE_UNIT cm
to be compatible with GMT4.3.1 and to set the degree format to vary from -180 to + 180 for longitude.
• CALPLOT/Utility:
• genplt: Make the length of the line segments slightly longer for the NO and DA legends.

#### May 8, 2016

• The current release is NP330.May-08-2016.tgz
• New program: mtinfo   This program is an implementation of Jost, M. L., and R. B. Herrmann (1989). A student's guide to and review of moment tensors, Seism. Res. Letters 60, 37-57. SRL_60_2_37-57.pdf. This program decomposes a moment tensor into different representations, e.g., isotropic, major and minor double couple, etc.   The program is invoked as
mtinfo -h
Usage: mtinfo
-XX Mxx -YY Myy -ZZ Mzz -XY -Mxy -XZ Mxz -YZ Myz
[-dc|-mc|-cc|-vd|-clvd|-a|-crack]

-XX Mxx -YY Myy -ZZ Mzz  Moment tensor elements in units of
-XY Mxy -XZ Mxz -YZ Myz    dyne-cm
Note -xy Mxy works as well as -yx Mxy
-dc               3 double couples
-mc               major couple
-cc               double couble and  CLVD
-vd               3 vector dipoles
-clvd             3 CLVDs
-crack            opening crack
-a                all of the above
Note: if none are set, -a is the default
-?                   This online help
-h                   This online help


As an example,
mtinfo -XX 0.264E+22 -YY -0.219E+22 -XY -0.195E+21 -XZ 0.188E+21 -YZ 0.100E+22 -ZZ -0.454E+21 > mtinfo.txt

produces the file mtinfo.txt.

#### April 28, 2016

• Corrected sac2eloc which got information from the header by using the call to brsac in sacsubf.f to a call to brsach. The problem was that brsac required a data array dimension that was not used since no data were required. This caused the program sac2eloc to fail for a file with many data points. Thanks to George França from Brazil.
• Fixed error in sacsubf.f n which the character header was not read using a call to brsach.

#### March 17, 2016

• gsac writesp now supports the options am, ph, rl, im to output the amplitude and pahse spectra and the real and imaginary parts of the complex spectrum. Previously only the amplitude spectrum was written. The result is placed in the same directory as the original file, but with .am, ph, .rl or .im appended to the name.
• Corrected the routine errelp in elocate. The code now follows  Flinn,  E. A. (1965). Confidence regions and error determination for seismic event location, Rev. Geophysics 3, 157-185. Note that the code does not multiply the error estimate by the F-statistic for the confidence region or by the t-statistic for the error in latitude, longitude, depth and origin time.

#### February 17, 2016

• Corrected an error in the FORTRAN routine etoh in SUBS/sacsubf.f and VOLVIII/src/elocate.f. This routine converts epoch time (time measured with respect to 1970/01/01 00:00:00.00 ) to a human form of Year, Month, Day, Hour, Minute, Second, Day Of Year. The routine did not give the proper conversion for times prior to the reference epoch.  This was discovered by writing MATLAB/Octave routines. The code snippets for testing the conversion are as follow:
• FORTRAN - using sacsubf.f for the time routines
        real*8 epoch        character str*40        integer year,month,day,hour,minute,date,doy        real second        year = 1969        month = 01        day = 01        hour = 18        minute = 27        second = 23.408        WRITE(6,*)year,month,day,hour,minute,second        call  htoe(year,month,day,hour,minute,second,epoch)        write(6,'(f20.3)')epoch        call etoh(epoch,date,str,doy,     1      year,month,day,hour,minute,second)        WRITE(6,*)year,month,day,hour,minute,second,doy,str        end        1969           1           1          18          27   23.4080009           -31469556.592        1969           1           1          18          27   23.4080009               1 1969/001 1969/01/01 18:27:23.408
• C - using csstim.c in PROGRAMS.330/SUBS
#include "csstim.h"main(){        double epoch ;        char str[40] ;        int year,month,day,hour,minute,date,doy ;        int sec, msec;        float second ;        year = 1969 ;        month = 01 ;        day = 01 ;        hour = 18 ;        minute = 27 ;        sec = 23;        msec = 408;        second = 23.408 ;        printf("%4.4d %2.2d %2.2d %2.2d %2.2d %2d %3d\n",year,month,day,hour,minute,sec,msec);        htoe2(year,month,day,hour,minute,sec,msec,&epoch);        printf("%20.3f\n",epoch);        etoh(epoch,&year,&doy,&month,&day,&hour,&minute,&sec,&msec);        printf("%4.4d %2.2d %2.2d (%3.3d) %2.2d %2.2d %2.2d %3.3d  \n",year,month,day,doy,hour,minute,sec,msec);}1969 01 01 18 27 23 408       -31469556.5921969 01 01 (001) 18 27 23 408
• C - using csstime.c in PROGRAMS.330/SUBS
#include "csstime.h"struct date_time T;main(){        double epoch ;        char str[40] ;        int year,month,day,hour,minute,date,doy ;        float second ;        T.year = 1969 ;        T.month = 01 ;        T.day = 01 ;        T.hour = 18 ;        T.minute = 27 ;        T.second = 23.408 ;        printf("%4.4d %2.2d %2.2d %2.2d %2.2d %6.3f\n",T.year,T.month,T.day,T.hour,T.minute,T.second);        mdtodate(&T);        htoe(&T);        printf("%20.3f\n",T.epoch);        etoh(&T);        timeprintstr(&T,str);        printf("%4.4d %2.2d %2.2d %2.2d %2.2d %6.3f %s \n",T.year,T.month,T.day,T.hour,T.minute,T.second,str);}1969 01 01 18 27 23.408       -31469556.5921969 01 01 18 27 23.408  1969001 Jan  1,1969 18:27:23.408
• Corrected error in computing dT/dh for direct arrivals in elocate.f.  This dT/dh is the change in travel time for a change in source depth. For the direct arrival for a general layered structure, there is no simple T(X) relationship. Instead we must compute an X(p) and T(p) where p is the ray parameter.  Formally

T(p) =

$X\left(p\right) = \sum _i \frac\left\{H_i p V_i\right\}\left\{ \left( 1 - p^2 V_i^2\right) ^ \frac\left\{1\right\}\left\{2\right\} \right\}$
When we try to compute a dT/dh for the direct arrival to a fixed distance, we  must note that the ray parameter changes for the ray to the station.  Thus we cannot just take a partial with respect to the source depth. To resolve this consider the general case when both the source depth, h, and the ray parameter, p, both change. Thus to first order

$T\left(p+\Delta p,h+\Delta h\right)=T\left(p,h\right)+\left(\partial T/\partial p\right)\Delta p+\left(\partial T/\partial h\right)\Delta h$
$X\left(p+\Delta p, h+\Delta h\right) = X\left(p,\Delta \right) + \left(\partial X/\partial p\right) \Delta p + \left(\partial X/\partial h\right) \Delta$

Since we do not want the X to change, the change in p is related to the change in h, and thus
$\left(\frac\left\{\partial T\right\}\left\{\partial h\right\}\right)_\left\{fixed X\right\} = \left(\frac\left\{\partial T\right\}\left\{\partial h\right\}\right)_\left\{fixed_p\right\} - \left(\frac\left\{\partial T\right\}\left\{\partial p\right\}\right)_\left\{fixed_h \right\} \left( \frac\left\{\partial X\right\}\left\{\partial p\right\}\right)_\left\{fixed_h\right\}^\left\{-1\right\} \left(\frac\left\{\partial X\right\}\left\{\partial h\right\}\right)_\left\{fixed_p\right\}$
• Added new feature to gsac cut command.  If the origin time and the distance are defined in the sac header, then the window can be defined using the syntax
VEL vref offset , where absolute time of the cut is given by Origin time + Distance/vref + offset.  One can also use a group slowness with P pref  offset
I have been doing group velocity cuts for moment tensor inversion by using awk to compute start and end times by the commands
	CUTL=echo $DIST$VEL $TB | awk '{print$1/$2 -$3}' 	CUTH=echo $DIST$VEL $TE | awk '{print$1/$2 +$3}' and then	CUT o ${CUTL} o${CUTH}
within gsac to implement something like cut o DIST/3.3 -30 o DIST/3.3 +70.   The other reason for introducing this option is to be able to look at segments of empirical Green's functions without having to also include the noise.  An example of the use of this command, some Green's functions were read and displayed.

 Use of cut to select a 30 second window about a group velocity of 3.5 km/s

• Added new feature to gsac xlim command.  If the origin time and the distance are defined in the sac header, then the window can be defined using the syntax VEL vref offset , where absolute time of the cut is given by Origin time + Distance/vref + offset.  One can also use a group slowness with P pref  offset. The next figure gives and example of the use of this command with the same traces displayed above. Note the markt on command was used to add group velocity markers to each plot.

NOTE: This only works with plot relative and will not work if the computed beginning of the trace occurs before the fist sample point.

 Use of the xlim command to display a 30 second window about a group velocity of 3.5 km/s

#### September 11, 2015

• The current release is NP330.Sep-11-2015.tgz.
• CALPLOT/Utility/genplt.f - Now have a -3 option where the data file has 3 columns, x, y, pen_color
• VOLIII/src/sdpegn96.f - Now has -E flag to plot Rayleigh wave ellipticity. Must be called as -R -E
• VOLIII/src/sdpspc96.f   - cleaned up the code.  If the -O observed -STA sta -COMP cmp  are given on the command line, the observed and theoretical spectral amplitudes are compared. The theoretical is attenuated to the distance of the observation, and then both are correct for geometrical spreading to the -DIST dist reference

If -O observed is not declared, then the theoretical spectra is attenuated using the  Q model (gamma) and corrected for geometrical spreading to the distance in -DIST dist.

In all case, the standard output shows all command lines and gives a tabulation of period spectral amplitude

If one does not want the Q effect, just use  -DIST 1.0 to get the spectra at 1 km for earthquake studies. The spectral at another distance is found by multiplying by sqrt ( 1 / desired_distance)
• VOLVI/src/hudson96.f - fixed error in subroutine modmerge in which the number of layers for the source and receiver models were not specified correctly.  Also reorked the logic on modrange for k < 0 return
• VOLVIII/src/elocate.f - added -C flag to output covariance matrix. Also documented and corrected the determination of the rotation angle for the confidence ellipse. Note the confidence ellipse is 1 sigma and the result must be multiplied by sqrt (2 var F(2,2,N-2))
• VOLII/src/do_mft4.c VOLII/src/sacamft96.f   - fixed XOR problem when displaying phase velocities. Also added a yellow background to the phase velocity pick
• VOLVIII/gsac.src - added new features
• Implemented the command map5. This has the same syntax as map. The two programs create the shell scripts map5.sh and map.sh, respectively. The former uses GMT5 (href="http://gmt.soest.hawaii.edu/"> http://gmt.soest.hawaii.edu/) and the latter used GMT3/4. There is a difference in the syntax of the codes. Both shell scripts are heavily commented to permit the user to modify the presentations. The following shows the output of scripts:
 GSAC› map kstnm on GSAC> sh map.sh GSAC› map5 kstnm on GSAC> sh map5.sh

The comments in map5.sh indicate how to change the precision of the latitude/longitude markers. The size of the PNG file is different since GMT3/4 attempts to create a BoundingBox based on the plot, while GMT5 does not, and instead relies on the media size.
• Added Smooth [ON|OFF] option to psp or plotsp
• Implemented the command psppk or plotsppk to permit interactive picking on the spectra. This was required to plot the amplitude spectrum of instrument responses and then to select the gain at a given frequency. This will be used to check the instrument responses of networks. The syntax is a subset of that used for the command ppk oe plotpk. One difference with respect to plotpk is that pick information is written onto the terminal window for other use.
The help psppk give the following:
SUMMARY:     Interactive spectra pickPlotSPPK     The options are   AMplitude          : Plot amplitude spectrum (default)   PHase              : Plot phase spectrum   PErplot [n|OFF]    : Plot n spectra per frame (default off)   Overlay [ON|OFF]   : Overlay all spectra (default off)   SMooth [ON|OFF]    : Apply 5 point smoothing to spectra                        (default OFF)   XLIn               : X-axis is linear   XLOg               : X-axis is logarithmic (default)   YLIn               : Y-axis is linear   YLOg               : Y-axis is logarithmic (default)   FMIn               :  Minimum  frequency for plot                         (default: DF for XLOG 0 for XLIN)   FMAx               : Maximum frequency for plot (default: Nyquist)   AMIn               : Minimum spectral amplitude to plot (default:                             0 for YLIN  and 0.0001 Amax for YLOG)   AMAx               : Maximum spectral amplitude to  plot  (default:        maximum)   Default            : Do not put up a phase menu, turn off markingDESCRIPTION:     The cursor responds to the following commands:     -         : compress frequency scale by factor of 2, recenter trace     _         : compress frequency scale by factor of 2, recenter trace     +         : expand frequency scale by factor of 2, recenter trace     =         : expand frequency scale by factor of 2, recenter trace     (space) : recenter trace     B         : move to the previous page of traces     L         : give frequency and amplitude of point beneath cursor     N         : move to next set of traces     O         : return to original trace scaling     Q         : end interactive trace picking     X         : Define trace window by entering X two times     Z         : define frequency range for linear regression to get t*SEE ALSO     PLOTSP
You will notice that the commands are a subset of the ppk or plotpk commands. To use this command, read the traces, fft, and then psppk. The output of the 'L' command in interactive cursor mode, is the following:
	0.030174  2.8778e+01 W60ABHRS55ABHR.WSTK
where the first column is the frequency in Hz, the second is the spectral amplitude, and the last is the name of the file. If the 'overlay' mode is sued, then the spectral amplitude for all spectra displayed is given, e.g.,
	0.030174  2.8778e+01 W60ABHRS55ABHR.WSTK	0.030174  9.3662e+01 W60ABHTS55ABHT.WSTK	0.030174  2.8613e+02 W60ABHZS55ABHZ.WSTKAn example of the use of the Z-Z option of pkppk is in the following image

This image was made with the shell script commands:
#!/bin/sh#####
#    Alaska event 2015/09/27 Mw=6.3 H=110, recorded#####gsac << EOF######    cut about pP#####cut a 20 a 40r *Zfileid namebackground on color 1color rainbowfileid namesort up dist#p relativertrtaper w 0.25######    differentiate ground velocity to get#    acceleration. For large earthquake#    acceleration spectrum is flat, any trend#    is due to tstardiffftpsppk xlin fmin 0 fmax 4 overlay on  smooth onqEOF


The terminal output is

lnA = a + bf : a=-10.794666 b=-1.466224 t*=0.466714  TPHLB02HHZ SMOOTH YES
lnA = a + bf : a=-11.126728 b=-0.823330 t*=0.262074  TPNVUS00BHZ SMOOTH YES
lnA = a + bf : a=-11.227097 b=-1.524203 t*=0.485169  SHPNN__HHZ SMOOTH YES
lnA = a + bf : a=-11.111622 b=-1.161991 t*=0.369873  GSCCI__BHZ SMOOTH YES
lnA = a + bf : a=-10.466206 b=-1.192656 t*=0.379634  LDFCI__BHZ SMOOTH YES
lnA = a + bf : a=-9.526356 b=-1.321453 t*=0.420632  OK029GS00HHZ SMOOTH YES
lnA = a + bf : a=-9.342464 b=-1.841290 t*=0.586101  WMOKUS00BHZ SMOOTH YES
lnA = a + bf : a=-9.339035 b=-1.845053 t*=0.587299  WMOKUS10BHZ SMOOTH YES
lnA = a + bf : a=-9.279791 b=-1.345738 t*=0.428362  OK025GS00HHZ SMOOTH YES

which gives the regression line, the t* and infomation about the trace and whether it was smoothed.

TO DO: Clean up the output of map/map5 so that event and stations are not repeated. This will make the scripts and resulting PostScript smaller in size

#### November 11, 2014

• VOLIII/src/sdpspc96.f  - made sure it worked. Also added -AZ az option so that a plot can be made without a data set

• VOLIII/src/sdprad96.f - At request of Rachel Noriega implemented -V   to output theoretical amp(azimith) radiation pattern

• VOLII/src/do_pom.c -  checked with gcc -c -Wall

• VOLIV/src/trftn96.f  - checked with gfortran -c -Wall

• VOLII/src/do_mft.c - Fixed so that if distance is negative, e.g., not defined if -12345, the program will not crash, but will rather flag the entry with a RED instead of BLUE color in the routine Pick_file()_

• VOLIV/src/jntpre96.f - Put in a safety check if the number of dispersion data exceed the defined limit, currently parameter (NM=12000)

• VOLIV/src/srfpre96.f - Put in a safety check if the number of dispersion data exceed the defined limit, currently parameter (NM=12000)

• VOLII/src/do_mft1.c - check with gcc -c -Wall

• VOLII/src/fmenunod.c - checked with gcc -c -Wall

• VOLII/src/do_pom1.c - checked  with gcc -c -Wall

• VOLII/src/do_pom4.c - checked  with gcc -c -Wall

• VOLIV/src/rftndr96.f - changed the screen output to provide sufficient space for the fit. Previously if the fit was negative, we might have something like ...T-99.5 . Now there is ...T -99.5. I wanted to be able to parse the 8'th or 11 columns of output to automate an external script to eliminate bad RFTN's from the inversion

• VOLIV/src/ramat.f - checked with gfortran -c -Wall

• VOLVI/src/jamat.f - hecked with gfortran -c Wall

• VOLII/src/do_mft4.c - checked  with gcc -c -Wall

• VOLVIII/gsac.src - Implemented QDP capability to plots  created using the refr or prs commands at the request of  Marcus Vinicius Lima, Federal University of Pampa, Brazil. This speeds up the display for very long time series. On many displays the number of pixels is about 1000 vertically. Thus trying to plot 30000 points, for example, requires a lot of computation but does not lead to anything useful on the display. This option can be invoked using the QDP command or by clicking on the QDP button under the refr display.

• Modified the gsac_prs.c and gsac_refr.c codes to work correctly when the distances are negative.  We can examine exploration data by setting lcalda false and then entering negative distances for the spread. When the ray parameter is now given for a reduced travel time plot, the plots look better

The images below show the effect of the change and also the "QDP" button.

 Normal display showing split profile Split profile after using a reduction velocity of 6 km/s
• The current release is NP330.Nov-11-2014.tgz.

#### June 10,2014

• Corrected the VOLII/src/nfmenu.h and VOLII/src/fmenunod.c  to handle the initialization of strings correctly. Under MacOS-X, the gcc 4.9 compiler led to SIGABRT after the phase match filter was used. This was traced to improper memory allocation for strings in the linked list. I read Kernightan and Ritchies, The C Programming Language 2nd Ediation, 1988, pages 139-143 and use strdup().  I had the pointers messed up.

• The current release is NP330.Jun-11-2014.tgz.

#### May 24, 2014

• Modified sacmft96 work around problems with improper station and component specifications in Sac files.  sacmft96 is called by do_mft for interactive analysis of group velocities and spectral amplitudes. The output file of sacmft96, mft96.disp, is used by do_mft for picking. Since a fscanf is used to read mft96.disp, the column entries must be exact and not corrupted.

The Sac KSTNM and KCMPNM header fields are 8 characters long. Usually these are of the form SLM and BHZ, respectively. Recently there were two examples of Sac files that had more complicated entries in these header fields, which cause the KSTNM to be effectively two columns, and this broke the scanf in do_mft. To see what happened, the command od -c was used to look at the binary contents of the Sac file:
od -c  COR_RNAK_RNSG.ZZ_p
         0000640  307 317 377 377 001  \0  \0  \0  \0  \0  \0  \0 001  \0  \0  \0
0000660  001  \0  \0  \0  \0  \0  \0  \0   R   N   S   G  \0   s
0000700    R   N   A   K  \0   5          \0   -   1   2   3   4   5
0000720    -   1   2   3   4   5           -   1   2   3   4   5        
Note that KSTNM is R N S G \0 s    and KEVNM is R N A K \0 5
Note the Null Characters in KEVNMC
To make everything work, sacmft96 now reformats the KSTNM and KCMPNM to replace non-printing character by an underscore, _, and to replace spaces between printing characters by an underscore.  Now mft96_disp contains RNSG_s for KSTNM. The previous line of mft96.disp looked like
MFT96 R U -1     0.10000    0.43319    0.00967    1.94096     36.9 0.7675E-05   39.531769 -119.850060   39.545757 -119.836517 0 10  0.1074       50.00     COMMENT: RNSG s  -12345   1970   1  0  0
and the new, corrected line, now is
MFT96 R U -1     0.10000    0.43319    0.00967    1.94096     36.9 0.7675E-05   39.531769 -119.850060   39.545757 -119.836517 0 10  0.1074       50.00     COMMENT: RNSG_s  -12345   1970   1  0  0
The change affects the mft96.disp and phv96.disp files created by sacmft96 as well as the selected dispersion contained in the XXX.dsp and YYY.phv files created by do_mft. The number of columns is now correct.  Thanks to J. Bonner and R. Noriega for noticing a problem and sending me the Sac files.
• Fixed problems in do_mft routines due to bad programming. The problem was a crash on MacOS-X on return from match filtering

• The current release is NP330.May-24-2014.tgz.

#### May 6, 2014

• Larry Baker, USGS, traced down the problem in compiling rdseed using gcc4.8 under MacOS-X  that we noticed August 22, 2013. This was in PROGRAMS.330/IRIS/rdseedv5.3.slu/Utilities/strlst.c.    rdseed now creates Sac files  using gcc4.8 and 4.9 compilers on on the current MacOS-X and Ubuntu 14.04.

• VOLVIII/gsac.src/csstim.h   - fixed the funding on the milliseconds.

• VOLVIII/src/pltsac.f - fixed error originally documented 09 AUG 2013

• VOLVIII/src/sacpol.f - added character sta2*8, comp2*8, cdate2*12. This is a new program to plot polarizations. A description will be posted.

• Many of the C programs required int main() to return the status.

• Programs changed are given in 20140506/bugfix.txt

• The current release is NP330.May-06-2014.tgz

#### April 28, 2014

• Larry Baker, USGS, found two problems in the low level graphics:

• PROGRAMS.330/CALPLOT/src/cmd/plotxvig.c line 2167 } else if(iy0 = iy1){  should be  } else if(iy0 == iy1){

• PROGRAMS.330/CALPLOT/src/clib/grd.c lines 496 and 516: At line 496 the old code is commented out and replaced:

•         /*
for( ; *s = (char )getc(fin); s++)
*/
for( ;; s++){
*s = (char )getc(fin);
if(*s == '\n')
break;
}

at line516

•                  /*  for(nread=0 ; s[nread] = (char )getc(fin);  nread++, nread < n){
*/
for(nread=0 ; nread < n ;  nread++){
s[nread] = (char )getc(fin);
if(s[nread] == '\n'){
break;
}
}
s[nread] = '\0';

• Corrections were also made PROGRAMS.330/CALPLOT/src/cmd/plotdriv.c and reframe.c

#### March 28, 2014

• Fixed bug in griding for plotsp. The lower part of the x-axis often had colored dashes instead of the solid black line.  Fixed by moving the gbox(...) after the axes were plotted.

#### March 13, 2014

• Added safety checks in the routine sacdatchk().  I created a simple sac file using wsac1() in sacsubc.c. The file had "-12345" for the station name, component name, etc. This broke the call to sprintf();

#### February 7, 2014

• Fixed minor error in fmplot.  When plotting first motions, fmplot lists the observations and compares the theoretical P-wave amplitude for that ray to the input first motion. If there is a sign difference, the phrase "Inconsistent" is printed.  This did not work correctly for  upper hemisphere projections.  The problem only affected the listing and did not affect the plots.

#### February 6, 2014

• Makefiles for VOLV/src for the program ttime96 use the isotropic model specification igetmod.f instead of tgetmod.f  This meant that only isotropic models could be used. By changing the entry in the Makefiles the code will actually computer travel times for TI and Isotropic models.

#### November 8, 2013

• redodate in VOLI/src changed. The offset for time computation was changed from  integer seconds to floating point seconds for more precision

• PROGRAMS/VOLIV/src:  Changed rftndr96.f output format to be able to display more precision in sample interval. surf96.f, rfth96.f and joint96.f have had the density relation modified to be able to go to about 900 km in depth. srfpre96.f changed to accomodate the relation between velocity and density to greater depths.

• PROGRAMS/VOLV/src/fmlpr.c:   Updated with latest reference to the NDK format.

• PROGRAMS/VOLIV/src/hspec96.f:   Corrected code at line 452 to with with Solaris f77 compiler

• PROGRAMS/VOLVIII/gsac.src/gsac_markt.c:  The markt  option now shows a 0.3 km/s group speed for looking at air waves

• PROGRAMS/VOLX/src/shallow96.f and shalpr96.f:  The density relation was modified to be consistent with surf96/joint96/rftn96

• do_mft and sacmft96 were significantly revised for the purpose of being able to get phase velocities from the interstation Green's functions obtained from noise cross-correlaiton. See the tutorial on this subject. There is no difference in the program output for group velocity analysis.

• The current release is NP330.Nov-16-2013.tgz

#### August 22, 2013

• Fixed a bug in the Makefile.OSX40 OSX40-32 in PROGRAMS.330/IRIS/rdseedv5.3.slu.   I had installed the latest gcc 4.8 compilers into /usr/local/bin on a new MacBook Air. rdseed seemed to work but did not create SAC files.  In the Makefiles I changed "CC = gcc " to "CC = /usr/bin gcc"

• I changed the routines gtrofa() and gtaofr() in the surface wave, receiver function inversion codes to add some higher velocities, so that the codes would work for the Earth to about 900 km

• The current release is NP330.Aug-22-2013.tgz

#### July 26, 2013

• Modified PROGRAMS.330/VOLVIII/gsac.src/gsac_plotsub.c to test if twin == 0.0 to make a cleaner plot. A one point time series is used by the new programs hstat96 hanal96 to provide the Greens functions for step source.

• Finally released the programs to create the Green's functions for static step change in the moment tensor for point sources. These programs are used in the sequence:  hprep96, hstat96/hanal96, fmech96/f96tosac in the same manner for creating synthetics, e.g., hprep96, hspec96, hpulse96. The new programs create an output file with the name file96 which can be converted to Sac files using the command cat file96 | f96tosac -G

• hsanal96 gives the static deformation (cm) for a moment 1.0E+20 dyne-cm step function.  This a uses an analytical expression for halfspace and a wholespace.

• hstat96 solves the problem of a source and receiver any where in a stack of elastic layers by numerical integration of an expression derived using propagator matrices. The program hsanal96 was written to provide an analytic solution for testing this wavenumber integration code.

• f96list reads a file96 waveform file and outputs the value of the first sample. This was useful for the static deformation programs. There is no on-line help. This program is run as

	f96list < file96
Units: CM/SEC
Created by None
Comment: ak135Q
R = 964.00000 AZ = 0.0000000
BAZ = 180.00000
ZDD 1.17266125E-12
........................
RHF -9.01416372E-11
THF 9.20795995E-10


#### July 22, 2013

• Change dimension in two parameter statements of spulse96 from NLAY=100 to NLAY 200 to be compatible with other surface-wave programs

#### May 22, 2013

• Corrected an error in do_pom (VOLII/src/do_pom.c) in procedure Inquire_file which to replace the line for(i=0;i< aargc; i++){ by for(i=0;i<=aargc; i++){.  This mistake caused the file list to be truncated by not showing the last file on the command line.  This error was detected when I attempted to us do_pom with only two traces.

• Because I was frustrated using the gsac * and / commands of the ppk (plotpk) command, I now permit the mouse scroll wheel to change the amplitudes.   Basically, because of the detault trace scaling, it was difficult to see the true first arrival and to pick the first motions without rescaling the amplitudes.  For regional earthquakes, I band pass the waveform between 1 and 5 Hz using the commands hp c 1 n 2 and lp c 5 n 2.  When I have the arrival suitably placed in the window, using the x x sequence, I then move the wheel forward to increase the amplitude and backward to decrease the amplitude.  The source code changes were in the following routines:  CALPLOT/src/plotxvig.c and VOLVIII/gsac.src/gsac_plotpk.c, changing one line in each file - a trivial change for major new capability.

• Modified pltsac to support an overlay mode ( -O  flag ) so overlay predicted onto observed waveforms. Current SLU moment tensor solutions use the overlay now to show the fits.

• f96tosac has a new command line flag for the naming of the output file ( -E ) which gives the Green function a name in the pattern DDdddHhhh.grn  so that the name can characterize the distance and depth to the nearest meter. This is useful to exploration data sets.  With this format the file name can represent distances between 0.000 and 99.999 km, and source depth between 0.0 and 9.999 km.

• Luis Rivera noted an error in the computations of first arrival times in hspec96. This occurred in the iterative computation for the direct (non-refracted) arrival, which had the entry

if(iter.lt.20)then
pnew = 0.999d+00 * pupper
else
Also the command line argument to subroutine frstar was changed by adding the logical parameter dogeom, which indicates whether to compute the geometrical spreading.  This is only needed for teleseisms in the time96 and hudson96 codes. The modification was to change 0.999d+00 to 0.99d+00. Changes were made to VOLIII/src/spulse96.f, VOLVI/src/hspec96.f, VOLVI/src/hudson96.f, VOLV/src/time96.f, VOLV/src/refmod96.f and VOLV/src/timmod96.f.
• Modified the output for the -z and  -r options of hrftn96. The seismograms are now for a unit area P- or SV-pulse incident in the halfspace. This change does not affect the receiver function, but permits a comparison of surface amplitude to incident amplitude. The comparison is actually done in the frequency domain, where a factor of two free surface effect is seen for vertically incident rays at zero frequency.

• In plotnps, the initialization 1 setlinecap 1 setlinejoin newpath is changed to 1 setlinecap 2 setlinejoin newpath. When plotting a almost horizontal line, the image seems to consist of filled overlaps half circles which give a jagged edge. This  was observed when the command plotnps -W10 was given to set the linewidth. It was bad for -W5 but not for -W50.  The 1 setlinejoin indicates that a round join be used when joining lines.  The 0 setlinejoin and 2 selinejoin indicate miter and bevel joins, respectively.  Hopefully the bevel join leads to a pleasing image.

• The current release is NP330.Sep-20-2012.tgz.

#### December 15, 2012

• Corrected compile warnings in VOLIII/src/sdisp96.f, VOLIII/src/tdisp96.f, VOLIX/src/cprep96.f, VOLIX/src/cpulse96.f and VOLVIII/src/sacfilt.f

#### December 10, 2012

• Modified surf96, rftn96 and joint96. Justin Ball of the University of Colorado was using the 'negative' layer thickness flag but this choice was not maintained for the temporary output files.  The 'negative' layer thickness is a way to set the reference level and to control the assignment of depths. It introduced for the atmosphere - solid earth problem. Rather than define source and receiver positions with respect ot the top of the model, one can define the reference level as the atmosphere - solid interface, and, as a consequence, positive values of source or receiver depths indicate positions in the solid while negative ones indicate elevations above the surfece.

While fixing the code, I use the opportunity to modify the dispersion display program srfphv96 so that plots of Love and Rayleigh dispersion use the same scales:

• Modified genplt to permit placing a legend within the plot. This is implemented using two new command line flags: -L lcmdfil -LPOS position. The command line syntax now has

-C cmdfil
-A acmdfil
-L lcmdfil

where the following are examples of the different files:

 cmdfil: 'HAM.RDD.val' 1 0.01 'CI' 'CUS.RDD.val' 1000 0.01 'CI' 'R41.RDD.val' 1025 0.01 'CI' 'S42.RDD.val' 1050 0.01 'CI' 'S43.RDD.val' 1075 0.01 'CI' acmdfil: 'HAM.RDD.val' 1 0.01 'CI' 0.04 'CUS.RDD.val' 1000 0.01 'CI' 0.04 'R41.RDD.val' 1025 0.01 'CI' 0.04 'S42.RDD.val' 1050 0.01 'CI' 0.04 'S43.RDD.val' 1075 0.01 'CI' 0.04 cmdfil: 'HAM.RDD.val' 1 0.01 'CI' 0.04 'HAM' 'CUS.RDD.val' 1000 0.01 'CI' 0.04 'CUS' 'R41.RDD.val' 1025 0.01 'CI' 0.04 'R41' 'S42.RDD.val' 1050 0.01 'CI' 0.04 'S42' 'S43.RDD.val' 1075 0.01 'CI' 0.04 'S43'

The postion parameter is a string of the form 'TL', 'ML', 'BL', 'TR', 'MR' or 'BR' for top-left, middle-left, bottom-left, top-right, middle-right and bottom-right, respectively. The legends plotted to the right are indented so that the entire string fits within the image frame.

Note that one can use all three commands. To make the plot shown, the following command line was used:

genplt -XMIN 100 -XMAX 1000 -YMAX 0.0001 -YMIN 0.000001 -XLOG -YLOG -L cmdfil -TX 'distance (km)' -TY 'Filtered Velocity'  -LPOS 'BL'

• Modified do_mft to be more efficient and added additional command line flags.

The new set of command line flags are as as follow:


-G           The default name of the  final dispersion file is StaCmp.dsp.
With this flag, the output dispersion file is the original filename with
.dsp appended, e.g., NMSLMBHZ01.dsp
-T           Show the 'Tomo' button
-DMIN dmin   only use traces with distances ≥ dmin km
-DMAX dmax   only use traces with distances ≤ dmax km
-11MIN 11min only use traces with IHDR11 ≥ 11min
-11MAX llmax only use traces with IHDR11 ≤ 11max

When processing there are several steps that are performed on different pages: select the trace on page 1, set the units on page 2, set the multiple filter analysis and display parameters on page 3 and then select the dispersion on page 4. Then repeat the entire process for the next waveform. Normally the parameters on pages 2 and 3 are not changed. The new Autofeed button permits one to skip pages 2 and 3. This small change saves a lot of time and mouse movements when working with large data sets.

The second change was specifically written for use with inter-station Green's functions computed by cross-correlating ground noise. Since gsac is used, the stack function saves the number of traces in the ihdr11 variable. The command line of do_mft can be used to use only traces with a range of ihdr11 values. The command line can also select distance ranges now, a feature of use for earthquake data. Finally a new -T flag introduces a Tomo button.

The Tomo button executes the script MFTDOOVERLAY that is in the users PATH. It is currently used to determine the tomographically predicted dispersion between the station and epicenter for earthquakes or two stations for cross-correlation, and then to overlay this onto the MFT image. The two figures show the buttons and the predicted dispersion.

The purpose of making this overlay is to assist in determining the proper period range for dispersion. For earthquake data, the spectral shape is often used as a guide, but the comparison may be useful. Another obtain is to overlay model predicted dispersion curves created from the -S flag of sdpegn96. The script MFTDOOVERLAY can be programmed to do a lot of things. For example a combination of gsacGMT can be used to create the file map.eps which can be viewed using gv or other PostScript view:

Tomography results and display scrips will be provided shortly.

• The current release is NP330.Dec-10-2012.tgz.

#### October 22, 2012

• Corrected an error in the gsac interpolation routine in  gsac_in.c. The problem was a spurious value for the last point. Shaoqian Hu fixed the error. The following image compares the original trace in red to the interpolated trace in blue. In order to make the presentation, the observed trace was read using a cut, and the cut version saved. The cut trace was then interpolated and saved. Finally both traces were red using the command cuterr fillz and cut with wider limits to focus on the behavior at the end. Thus is why both traces are zero at the end. The problem has long been noticed in waveform plots from the moment tensor inversion - however since only one point was involved, the moment tensor inversion results do not change.

 Original code Corrected code

#### September 20, 2012

• Upgraded the distributed IRIS program rdseed from  5.2 to 5.3.  In order to support the CPS architectures, e.g., CYGWIN (Windows), MacOS-X, LINUX and Solaris, changes were made to these codes. The changes are given in the links in this section.Only minor changes are made. The comments in the Sac pole-zero files are a superset of rdseedv5.3 and rdseedv5.0.slu.

• Main/output_resp.c -  modify to permit extra output for use in seismic network QC

• Main/output_sac.c - modified comments in Sac pole-zero file

• Main/summary.c - changed the call to fabsl to fabs since there was no fabsl in CYGWIN and since the use of the code only required double and not long double.

• Parsers/parse_50.c - eliminate non-printing characters in the site name

• Utilities/log_errors.c - add #include <string.h>

• Detailed changes

• Modified the gsac ppk to recognize the use of the 'm' key as the same as the '*' key to make the trace amplitudes larger.

• The current release is NP330.Sep-20-2012.tgz.

#### September 9, 2012

• The include file calplot.h now has prototypes for the routines void gcont(float x, float y);
void gmove(float x, float y);
(Thanks to Shaoqian Hu)

• Fixed the bsort and dbsort subroutines in sprep96.f, sdisp96.f, tprep96.f, tdisp96.f, scomb96.f and tcomb96.f. The problem was that the bsort routine sorts and array, and then keeps only the unique values. Thus the length of the array returned can be smaller than the initial array.  The last line of the subroutinhe was written as  n=n-m instead of nn=n-m. (Thanks to ruoshan at USTC)

• The current release is NP330.Sep-09-2012.tgz.

#### August 6, 2012

• Fixed an error in the Setup script so that the rdseed now compiles under Cygwin. The current release is NP330.Aug-06-2012.tgz

#### July 26, 2012

• hpulse96 now uses list direction IO instead of formatted IO for user defined source time function

• The program srfker96 is now part of the package. This program reads the temporary files from surf96 to list all partial derivatives as a function of depth.  A tutorial on the use of this program is available at the tutorial web page http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/SWKERNELS/ .

• After creating the above tutorial, I realized that I cannot compute the partials of the anelastic attenuation coefficient γ   with respect to changes in P- and S-wave velocity or layer thickness. These partials created in subroutine PROGRAMS.330/VOLIV/src/jsamat.f   now set dgdv = 0 for each layer.

• The furrent distribution is NP330.Jul-26-2012.tgz

#### April 9, 2012

• Upgraded the distributed IRIS program rdseed from  5.0 to 5.2.  In order to support the CPS architectures, e.g., CYGWIN (Windows), MacOS-X, LINUX and Solaris, changes were made to these codes. The changes are given in the links in this section.Only minor changes are made. The comments in the Sac pole-zero files are a superset of rdseedv5.2 and rdseedv5.0.slu.

• Main/output_resp.c -  modify to permit extra output for use in seismic network QC

• Main/output_sac.c - modified comments in Sac pole-zero file

• Main/summary.c - changed the call to fabsl to fabs since there was no fabsl in CYGWIN and since the use of the code only required double and not long double.

• Parsers/parse_50.c - eliminate non-printing characters in the site name

• Utilities/log_errors.c - add #include <string.h>

• Detailed changes

#### March 26, 2012

• Changes to gsac:
• #### outcsvcommand  adds a column header which now indicates Freq(Hz) for spectral files in addition to the previous Time(sec) for time series.

• prs (plotrecordsection) permits annotation using the title command. However this is not perfect since the location of a title at the top of the page conflicts with the annotate option of prs and the location on the left conflicts with the y-axis label

• prs command now behaves better with the use of the prs title title_text command to supersede the normal axis label. The annotation is required since the Sac header only permits a few of the variables to have associated text strings, and if the USER1 happens to be set for a user defined Angle, one can arrange the record section accordingly with the command prs user1 title "Angle (deg)"
• tregn96, tdisp96, sregn96, scomb96, tcomb96: cleaned up compiler warnings

• hspec96: Fixed a problem with flagging a receiver position is a fluid to compute the pressure fields in the fluid for the different sources. This only appeared with an atmosphere model was used which required refdep to be other than 0.0

• tspec96: Fixed a problem with flagging a receiver position is a fluid to compute the pressure fields in the fluid for the different sources. This only appeared with an atmosphere model was used which required refdep to be other than 0.0. Hoever I have to check the travel times for a fluid.

#### February 8, 2012

• Modified genplt to have additional plotting controls on the size for the symbol shape. This means that the other program tgenplt is deprecated. The syntax of the current program and an example of its use is given in the line GENPLT/nindex.html. genplt will run as before. The new capabilities are invoked with the -A acmdfil command-line option.

• The current distribution of Computer Programs in Seismology is now NP330.Feb-09-2012.tgz

• #### January 14, 2012

• Increased the size of the outstr string in do_mft to prevent overflow with long file names

#### November 23, 2011

• Fixed an error in the synchronize command of gsac. The problem arose when reading just the headers and then synchronizing.  The code has been fixed. The current version of gsac is now V1.1.44 21 NOV 2011.

• The command line parsing of the program fmplot was rewritten,

• The current distribution of Computer Programs in Seismology is now NP330.Dec-01-2011.tgz

#### November 1, 2011

• In response to user requests the following changes were made to gsac. The current version of gsac is now V1.1.43 31 OCT 2011.

• At the request of H. Benz, gsac now sets a header value corresponding to the time of the maximum (depmax) and minimum (depmin) amplitudes.  The previously unused float/real header values in the 70 position real header, FHDR65 and FHDR66, are now called TIMMAX, TIMMIN, respectively. This request was made to use the times of the maximum and minimum points of a trace in external programs.

• These other codes were changed to reflect these new parameters:  scmxmn routine in SUBC/sacsubc.c, SUBS/sacsubc.h, SUBS/sacsubf.f . The programs that use these general Sac I/O routines were re-compiled.

• A new command, ricker, was added to gsac.  This permits one to convolve a time series with a Ricker wavelet of a given frequency. This was implemented because of difficulties with the spulse96 -f wavelet_file   and hpulse96 -f wavelet_file implementation.

• The current use of Sac header variables by Computer Programs  in Seismology is as follows:

• USER0 - Gaussian filter parameter, alpha, used in saciterd, saciterdc, rftn96 and joint96

• USER1 - minimum period of filter. This is introduced so that users do not falsely accept dispersion values outside the bandpass of a filter

• USER2 - maximum period of filter.

• USER4 - ray parameter for joint96 and rftn96.  Set from command line by saciterd and saciterdc.

• USER5 - percent fit value. Appears in pltsac, wvfgrd96, wvfmt96, wvfmtd96, hrftn96, trftn96, joint96 and rftn96.

• USER9 - time shift for waveform fit. Used by pltsac and set by wvfgrd96, wvfmt96 and wvfmtd96.

• TIMMAX - offset with respect to reference time indicating time of trace maximum

• TIMMIN - offset with respect to reference time indicating time of trace minimum

• IHDR20 - set by gsac as part of the ppk qc command to indicate whether a trace passed QC.  I obtain this value from a sac file in a shell script with a call to saclhdr, which is then use in an awk script.

• The current version of Computer Programs in Seismology is NP330.Nov-01-2011.tgz.

#### October 20, 2011

• C. Buffoni noted a problem with the presentation shaded area seismic trace plots in PostScript using the command plotnps. It took a while to track this down, but the problem was repaired. In addition the shaded area plot correctly for the command plotxvig -R, although this would rarely be used.  In addition plotgif -R works correctly.

 OLD (problem) NEW (Repaired)

• Modified the CALPLOT functions pend() and gend() to reset the pen width to the default 0, and the pen to 1 (black). This assists merging files.

• NOTE: The S-wave receiver function, e.g., Z/R, computed by hrftn96 is not stable.  This may be due to the problem of deconvolving with a wavelet that is not minimum phase. This is a major issue that is deferred to a future update.

• The current version of Computer Programs in Seismology is NP330.Oct-27-2011.tgz.

#### October 13, 2011

• The mt (moment tensor) command of gsac now supports the specification of the moment tensor source on the command line, e.g.,

• mt to ZRT MXX 0.0000000 MXY -1.00000002E+20 MXZ 0.0000000 MYY 0.0000000 MYZ 0.0000000 MZZ 0.0000000 AZ 25. FILE /backup/rbh/PROGRAMS.310t/GREEN/CUS.REG/0100/030000100

• The current version of gsac is V1.1.42 13 OCT 2011

#### September 23, 2011

• Major change in gsac:

• refr command - filters are 2-pole 1 pass and not 2-pole 2 pass - thus they are causal!!

• plot command - for many traces displayed sequentially using the perplot command the y/n/b query is now y/n/b/q so that one can quit

• mt command - fixed code so that the t0 and t1 arrival times of SV and SH are correctly read and then saved in the headers. This arose during conversion of horizontals to N E

• prs command - do not annotate trace if KSTNM is -12345 and add LAST kolor FIRST color to permit different colors. I wanted to compared one observed trace to many predicted traces using different colors for observed and predicted.

• csstime.c - MAJOR change. I used an incorrect rounding when going from epoch time to integers NZYEAR NZJDAY NZHOUR NZMIN NZSEC and NZMSEC. This inadvertently messed up the seconds by one unit.

• Minor changes to sacmat96, sepegn96.

• Changes to the VOLIV surface-wave inversion programs to permit computation of du/da and dc/da, partials with respect to P-wave velocity. these are not used in the surface-wave of joint inversion; the changes were made for the design of a new inversion of P, S first-arrivals and dispersion for shallow exploration

• hrftn96 in VOLVI changed the sign of the Z and R synthetics for SV incident to agree with expected particle motion.

• In VOLVII, wvfmt96, wvfmtd96  sets USER5 to the percentage of fit, e.g., 100 * correlation coefficient squared

• In VOLVIII/src  changed output format for the file elocate.sum in elocate

• In VOLVII/src introducted new program saciterdc which is C version of the FORTRAN saciterd for speed, and easy transliteration of routines to java. The -Q flag suppresses output at each iteration

• in SUBS  put in an error check for the routine bwsac in the file sacsubc.c

#### October 24, 2010

• Fixed bugs in the interpolate function of gsac which affected the convolve comand. The current version of gsac is GSAC - Computer Programs in Seismology [V1.1.39 24 OCT 2010].

• The current release of Computer Programs in Seismology is NP330.Oct-24-2010.tgz.

#### October 12, 2010

• The current release of Computer Programs in Seismology is NP330.Oct-12-2010.tgz

• Corrected Makefiles for OSX in VOLVIII/gsac.src/  so that arrow keys work for gsac on OSX

• Modified the Setup script so that the help reads:

• 	Usage: Setup SOL WIN32 SOL-EGCS LINUX LINUX64 CYGWIN OSX
SOL        SUN Solaris Compilers
SOL-GNU    SUN Solaris with gcc/g77 compilers
CYGWIN     CYGWIN 98/NT/2K/XP gcc/g77 Compilers
CYGWIN40   CYGWIN 98/NT/2K/XP gcc/gfortran Compilers
LINUX40    Linux with gcc/gfortran compilers
LINUX6440  Linux 64 bit with gcc/gfortran compilers
OSX40      Apple  with gcc/gfortran compilers
OSX40-32   Apple  with 32-bit gcc/gfortran compilers
Older systems - not modern
LINUX      Linux with gcc/g77 compilers
LINUX64    Linux 64 bit with gcc/g77 compilers
OSX        Apple  with gcc/g77 compilers

The purpose of the OSX40-32 flag is to force a 32 bit compile. This is useful on some Mac's where the gfortran and gcc compilers such that one is 32 bit and the other 64 bit

• Fixed hspec96 by adding a -ZSRC zsrc and -ZREC zerc flag that will define the maximum depth of the source/receiver location models to patch onto the teleseismic model.  Before this change the halfspace was replaced by a 1 km thick layer, which messed up the resulting model. Creating a merged model is not simple scientifically irrespective of the programming error.

#### October 3, 2010

• Upgraded the distributed IRIS programs evalresp from version 3.2.40 to 3.3.3 and rdseed from  4.8 to 5.0.  In order to support the CPS architectures, e.g., CYGWIN (Windows), MacOS-X, LINUX and Solaris, changes were made to these codes. The changes are given in the links in this section. This upgrade was required in order to continue to use the most recent SEED volumes.

The two programs are installed in PROGRAMS.330/bin with the executable names evalresp and rdseed. The IRIS distribution found in PROGRAMS.330/IRIS/rdseedv5.0   provides the executables rdseed.linux (64 bit), rdseed.MAC.x86 (64 bit), rdseed.powerpc, rdseed.sparc (64 bit), and rdseed.sun.x86 (64 bit). If you wish to use these programs, you must install them in your system and adjust the path to find them.

• evalresp  - The SLU version has the fixes suggested by Eric J. Haug of SLU  in May, 2008 to handle continuation for FIRST with a large number of coefficients. Basically we need to get past the tendency of some SEED writers to put the following phase in the RESP file "No Abbreviation Referenced" for Response in units lookup and Response out units lookup. This is a proble with the CMG-5TD accelerometer which has a very, very FIR stage.

• rdseed - Several changes were made.

• Another change relates to the strupr routine. Since this  is actually provided by some libraries, we must avoid a double definition.  By making this global change the rdseed version of this is preserved and there is no conflict with other definitions. Thus all instances of strupr are changed to Strupr and a definition is placed in Include/rdseed.h. Changes were made in Decoders/process_data.c, Decoders/process_event_requests.c, Main/ah_resp.c, Main/output_css.c, Main/rdseed.c and Main/summary.c.

• Utilities/alloc_linklist_element.c -   malloc wants a size_t cast. This is changed from an int

• Utilities/get_date.c - Change #include <time.h> add #include <sys/time.h> which is required to define gettimeofday() [ LINUX, Solaris. CYGWIN]

• Ah/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Decoders/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Parsers/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Printers/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Snoop/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Utilities/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Parsers/parse_50.c - get rid of non-printing characters in the output

• Main/output_sac.c - Place comments into the pole-zero file giving information about the station, network, date, etc. rdseedv.5.0 incorporates this. The code changes add some additional information about the instrument and also test the output strings.

• Main/output_resp.c - Send information about the network to stderr, e.g., station, network, component, location, on_date, off_date, latitude, longitude, elevation, response file : rdseed -f seed_volume -R > out 2>&1   which using the sh or bash shells.
(Sample stderr output)

• Updated the manual entry in Main/rdseed.1 to make the command line flags more complete on the first page. The pdf version of this document is in PROGRAMS.330/DOC/IRIS.pdf where you will find evalresp.pdf, rdseed.pdf and SEEDManual_V2.4.pdf (no change from November 17, 2007)

• Compiler warnings: : Significant warnings appear using gcc -c -g -Wall instead of the normal gcc -c in the makefifles. Most of the warnings are due to the fact that the current rdseed is not written in ANSI C but rather in the original Kernighan and Ritchie style. Someone else should clean up this code.

• The modified rdseed was successfully compiled on both 32 and 64 bit Intel LINUX machines using the GNU compilers, on a SPARC Solaris machine using SUN compilers, on an Intel Windows machine under the CYGWIN environment, and under MacOX-X using the GNU compilers. The only compiler warnings were the following:

• cc -I../Include -c -g -c -o decode_steim.o decode_steim.c

decode_steim.c: In function decode_steim':
decode_steim.c:231: warning: integer overflow in expression
decode_steim.c:231: warning: integer overflow in expression
• on just the 64-bit LINUX compiler (note that the -w32 compiler flag was not used since the change to Include/sac.h makes this unnecessary) there are a number of errors associated with printf formats (Output of gcc -c).

#### September 30, 2010

• Modified gsac to implement the convolve command.  This required significant effort because positioning and time delays are a problem. The syntax of the new command differs from that of Sac and  is:


Con Z [B|O|A|T0|T1|T2|T3|T4|T5|T6|T7|T8|T9] File external_sac_file [Suffix suffix]

One marker must be set to define the zero time of the pulse in order to prevent incorrect delays. The following example creates an impulse, sets "O" to the start of the impulse, convolves with the internal function "Triangle" and then uses the "Convolve" command to do the same:

gsac << EOF
#####
#       define the filter
#####
fg imp delta 0.1 npts 256
w imp.sac
r imp.sac
# set the o marker
ch o 0
w
triangle w 1.0
w tr.sac
#####
#       define the impulse with a different sample rate
#####
fg imp delta 0.05 npts 512
w nimp.sac
r nimp.sac
conv z o f tr.sac s ".con"
w
q
EOF

A comparison of the files tr.sac and nimp.sac.con  shows that they are the same and have the same spectral amplitudes, except at high frequency because of the different sample rates.

• The current version of gsac is now V1.1.38 30 SEP 2010.

#### September 20, 2010

• Fixed some declaration errors in sacpsd

#### September 17, 2010

• Introduced the simple plotting program genplt

• Added a new program wvfdly96 for use in the source inversion. This program interprets the time shift used for regional moment tensor inversion to a change in origin time and location. This was necessary when analyzing events whose locations are poorly constrained

#### August 30, 2010

• Surface-wave programs in VOLIII/src were rewritten:  sregn96, spulse96, slegn96, sdpegn96, sdpdsp96, slat2d96, sdisp96, sdpder96 and sdpsrf96. The purpose for the rewrite was to align the code with the development in the class notes

• Added new programs to VOLIII/src for transverse-isotropic (TI) media. The new programs provide TI synthetics and phase-velocity partials with respect to the four vertical and horizontal P- and S-wave velocitys and the parameter η = F/(A - 2L). These new programs are called tcomb96tdpder96,  tdpsrf96,  tprep96,     tsregn96,  tdisp96,  tdpegn96,    tlegn96,  tpulse96  and tregn96.  The programs work together in the same manner as the original codes for isotropic models:

	tprep96 -M model.TI -L -R -NMOD 500 -d dfile -HS 10 -HR 0
tdisp96
tlegn96
tregn96
tpulse96 -V -p -l 4 -d dfile > tifile96
fprof96 -A < tifile96 

Of course isotropic models can be used with these programs, also though the computations will take longer than for the corresponding isotropic codes.

#### August 20, 2010

• Modified gsac in the following manner. The prs (plotrecordsection) command supports the use of the MAG field as an axis. The sort command now permits a sort on MAG.

• The current version of gsac is now V1.1.37 20 AUG 2010.

#### July 15, 2010

• The mt (moment tensor) command in gsac was modified to permit point forces with the command syntax FN fn FE fe FD fd where the force(units of dynes) is described in the north, east and down directions, respectively.  The output will be in meters.  The current version of gsac is now V1.1.36 15 JUL 2010

• For the purposes of documenting the menus, gsac accepts a ~ interactive input under the ppk  and refr commands. This will create a unique CALPLOT graphics file with name DUMPxxx.PLT which will show the menu buttons and the waveforms.

#### June 17, 2010

• The mt (moment tensor) command in gsac was modified to permit an explosion source, mt iso mw 5.0 FILE fileproto. In addition the generation of separate Z, R, T, N or E components can also the new syntax mt TO UZ  or in general one of UR, UZ, UT, UN or UE in place  of the previous R, Z, T, N or E, respectively. The current version of gsac is now V1.1.33 17 JUN 2010.

#### June 14, 2010

• Fixed login error in the gsac command sync which failed when there were P picks in the file. The code is a bit cleaner now. The current version of gsac is now V1.1.32 13 JUN 2010

#### June 6, 2010

• Added a new command to gsac.   The shift command moves the position of a trace in time, adjusting the A T0 ... T9 markers.  The syntax is shift Fixed Amount. If amount is positive the trace is plotted later in absolute time.  This command was introduced to fupport finite fault modeling, which requires that subfaults fire a later times as the rupture front excites the subfault.   I previously manipulated the O (origin time marker) and did some fancy cutting that was not intuitively obvious. The current version of gsac is now 1.30 06 JUN 2010.

• Modified output of elocate to provide a command line that can be used with gsac or sac2000. The line appears as  ch o gmt YEAR DOY HR MN SEC MSEC

#### May 20, 2010

• Modified iputmod.f in PROGRAMS.330/SUBS. Created igetmod.c, iputmod.c, imodel.h to implement 'getmod', 'putmod' in C. The igetmod.c source also includes a routine  getdval  which search the model to get the model parameters at a given depth within the model.

• Added a conversion to ZNE in the gsac mt command in order to make Z N and E components. Previously we could make Z R and T components. I needed this for short distance finite fault simulation for which the back azimuth to each sub-fault can change. This makes the current version number of gsac  [V1.1.30 19 MAY 2010]

•  Added rot3 to zne to the gsac rot3 command.

•  At the suggestion of Carol Bryan, fixed a bug in the deviatoric constraint in the program wvfmtd96  in PROGRAMS.330/VOLVII/src/wvfmtd96.f

•  pltsac - added -TMIN fmin -TMAX tmax command line arguments to make better plots.

•  shwmod96 - added dashed line option - but be careful if used with automatic -K -1 command line flag.

• sdpegn96 - added -NOBLACK command line flat to prevent black outline around the observed dispersion for compatibility with sdpdsp96.

#### April 7, 2010

• Slight modifications to plotnps. When I implemented the background fill to implement the -I flag, I also had a default forced white backgound when the -I flat was not invoked. While this was the desired effect for converting an EPS to graphics format such ad PNG using ImageMagick  so that there would be no transparency, this actually broke the inclusion of EPS files using groff since the background fill erased the page.  As an example of the use of the current flags, this is how I now work with the ImageMagick command convert and the GraphicsMagick command gm convert:

	plotnps -F7 -W10 -EPS -K -BGFILL < SACPSD.PLT > 2.eps
convert -trim 2.eps im2.png
gm convert -trim 2.eps gm2.png
plotnps -F7 -W10 -EPS -K -I < SACPSD.PLT > 3.eps
convert -trim 3.eps im3.png
gm convert -trim 3.eps gm3.png
• The current version of Computer Programs in Seismology is NP330.Apr-07-2010.tgz

#### March 18, 2010

• Modified the generation of EPS by  plotnps to a) paint the background to avoid transparency in ImageMagrick, and b) to permit an inverted background. The inverted background is useful for presenting some colored lines, e.g., yellow, which are difficult to see against a white background. The following images show the syntax of the command abd also the result of converting the EPS to a PNG file using GraphicsMagick. The ImageMagick command trims to the smallest size
 plotnps  -F7 -W10 -EPS -K < GRAYSC.PLT > t.eps gm convert -trim t.eps graysc.png plotnps  -F7 -W10 -EPS -K -I < GRAYSC.PLT > t.eps gm convert -trim t.eps grayscinv.png

#### February 28, 2010

• Introduced a new program sacpsd to compute power spectral density of ground acceleration noise for use in seismic network quality control. This program is in PROGRAMS.330/VOLVIII/src. The program uses long data segments in Sac file and the acceleration sensitivity output of the IRIS evalresp from the RESP files created by the IRIS rdseed. Documentation describing the use of this program is given in sacpsd.pdf    and the shell scripts used in the document are given in the file  SACPSD.tgz    (which is unpacked using the command  sequence gunzip -c SACPSD.tgz | tar xvf -   .   Note that some errors were corrected on February 28, 2010.

#### February 16, 2010

• Made slight changes to two routines in the current version of rdseed:

• PROGRAMS.330/IRIS/rdseedv4.8.slu/Main/output_seed.c

• PROGRAMS.330/IRIS/rdseedv4.8.slu/Parsers/parse_11.c

to move two declarations to the top of  update_type74 in output_seed.c and in parse_type11 in parse_11.c.  The Solaris cc compiler
did not like the declarations to be intermixed with code farther down.

• Changed the Main/makefile, Printers/makefile, Snoop/makefile  and Utilities/makefile to remove a CC=gcc line that I had entered previously

• Made changes to elocate in PROGRAMS.330/VOLVIII/src to facilitate the exchange of data for use in other location programs and also for easier change of source parameters in gsac. The extra output (inlooks like

	STA   IWT       ARRIVAL TIME   PhIDQL PHASE  FM CHAN    LAT       LON      ELEV
OK001   0 20100213053058.250    2 0 i S        X  Z   35.5611  -97.2895    350.00
OK001   0 20100213053057.219    1 0 i P        C  Z   35.5611  -97.2895    350.00
OK002   0 20100213053100.168    2 0 i S        X  Z   35.5493  -97.1966    354.00
...    ...   ...   ...   ...
HYPO71 Quality   :                  BA
Gap              :                  84              deg

ch evla    35.5331 evlo   -97.2982 evdp       5.52
ch ocal   2010 02 13 05 30 55 968

The last entries are used by cut here and pasting into the gsac command line.

#### January 11, 2010

• Modified the lp, bp, br and bp commands of gsac to implement a Chebyshev Type - I filter using the C1 command argument.  Following sac/sac2000 the lowpass filter prototype is defined by the number of poles, the corner frequency, fc, a transition bandwidth, transbw, and a desired attenuation factor at the frequency (1+transbw)*fc.  References on this filter are given in wikipedia, but the reference providing all information, including the gain normalization constant is Rorabaugh (1997).  Implementation of the Type-II filter, also known as the Inverse Chebyshev filter (C2 command argument in sac/sac2000 usage) was not implemented since data structures in the file gsac_filt.c would have to be completely redesigned.

Digital Filter Designers Handbook with C++ algorithms
C.  Britton Rorabaugh 2nd Edition
McGraw Hill, New York, 479 pp 1997  Chapter 5

• While testing the Chebyshev implementation, an error in the implementation of the bandreject filter for a single-pole stage  changing

} else if (filt_lhp == BANDREJECT){
wc = fwarp(Wc,dt);
wl = fwarp(Wl/ffac,dt);
wu = fwarp(Wu/ffac,dt);
num[0] = wl*wu*dt*dt;
num[1] = 0.0;
to

} else if (filt_lhp == BANDREJECT){
wc = fwarp(Wc,dt);
wl = fwarp(Wl/ffac,dt);
wu = fwarp(Wu/ffac,dt);
num[0] = wc*wc*dt*dt;
num[1] = 0.0;

#### December 4, 2009

• Upgraded the distributed IRIS programs evalresp from version 3.2.38 to 3.2.40 and rdseed from  4.7.5 to 4.8.  In order to support the CPS architectures, e.g., CYGWIN (Windows), MacOS-X, LINUX and Solaris, changes were made to these codes. The changes are given in the links in this section. This upgrade was required in order to continue to use the most recent SEED volumes.

The two programs are installed in PROGRAMS.330/bin with the executable nanames evalresp and rdseed. The IRIS distribution found in PROGRAMS.330/IRIS/rdseedv4.8   provides the executables rdseed.linux, rdseed.ppcMac, rdseed.intelMac, rdseed.Sun86. rdseed.SunSparc. If you wish to use these programs, you must install them in your system and adjust the path to find them.

• evalresp  - The distributed code has an #include <malloc.h>. This is deprecated in MacOS-X to <malloc/malloc.h> which means that the code will not compile.  Since modern compilers provide the prototypes for malloc in <stdlib.h>. I changed alloc_fctns.c, parse_fctns.c, regexp.c and string_fctns.c by commenting out the #include <malloc.h> and ensuring the presence of a #include <stdlib.h>.
Detailed changes.

NOTE THAT evalresp-3.2.40 will configure properly for the MacOS.

• rdseed - Several changes were made.

• Include/sac.h: The most important was to change the use of long in the definition of the integer SAC header values to int. Becore computers were truly 32 bit, short meant a 2 byte integer, long a 4 byte integer, and int could mean a  byte integer. Currently, int means 4 bytes.  A problem arises on 64 bit machines, for which a long means an 8 byte integer.  One way around this is to use the -m32 flag on the compiler.  CPS changes the long to int.

• Another change relates to the strupr routine. Since this  is actually provided by some libraries, we must avoid a double definition.  By making this global change the rdseed version of this is preserved and there is no conflict with other definitions. Thus all instances of strupr are changed to Strupr and a definition is placed in Include/rdseed.h. Changes were made in Decoders/process_data.c, Decoders/process_event_requests.c, Main/ah_resp.c, Main/output_css.c, Main/rdseed.c and Main/summary.c.

• Utilities/alloc_linklist_element.c -   malloc wants a size_t cast. This is changed from an int

• Utilities/get_date.c - Change #include <time.h> add #include <sys/time.h> which is required to define gettimeofday()

• Ah/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Decoders/makefile - remove CC=gcc  since not all compilers use gcc

• Parsers/makefile - remove CC=gcc  since not all compilers use gcc

• Parsers/parse_50.c - get rid of non-printing characters in the output

• Main/output_sac.c - Place comments into the pole-zero file giving information about the station, network, date, etc. The SAC pole-zero file is created using the command  rdseed -f seed_volume -p
(Sample SAC pole-zero file)

• Main/output_resp.c - Send information about the network to stderr, e.g., station, network, component, location, on_date, off_date, latitude, longitude, elevation, response file : rdseed -f seed_volume -R > out 2>&1   which using the sh or bash shells.
(Sample stderr output)

• Updated the manual entry in Main/rdseed.1 to make the command line flags more complete on the first page. The pdf version of this document is in PROGRAMS.330/DOC/IRIS.pdf where you will find evalresp.pdf, rdseed.pdf and SEEDManual_V2.4.pdf. (no change from November 17, 2007)

• Compiler warnings: : Significant warnings appear using gcc -c -g -Wall instead of the normal gcc -c in the makefifles. Most of the warnings are due to the fact that the current rdseed is not written in ANSI C but rather in the original Kernighan and Ritchie style. Someone else should clean up this code.

• The modified rdseed was successfully compiled on both 32 and 64 bit Intel LINUX machines using the GNU compilers, on a SPARC Solaris machine using SUN compilers, on an Intel Windows machine under the CYGWIN environment, and under MacOX-X using the GNU compilers. The only compiler warnings were the following:

• cc -I../Include -c -g   -c -o decode_steim.o decode_steim.c
decode_steim.c: In function decode_steim':
decode_steim.c:231: warning: integer overflow in expression
decode_steim.c:231: warning: integer overflow in expression

for all compiles, and

cc -I../Include -c -g   -c -o parse_30.o parse_30.c
parse_30.c: In function parse_type30':
parse_30.c:96: warning: cast to pointer from integer of different size

on just the 64-bit LINUX compiler (note that the -w32 compiler flag was not used since the change to Include/sac.h makes this unnecessary).

#### October 8, 2009

Summary of changes made since July 22, 2009

• gsac is now at version 1.1.28 30 SEP 2009. The following changes were made:

• The command plotpk (ppk) regional or plotpk (ppk) teleseism now had an Undo button which removes the last phase pick made. This was introduced to make phase picks for location easier to use

• The command merge had an improper index at line 214 which caused memory corruption.

• The subroutine gsac_sac.c, which is used for low level reading and writing of the SAC files, had an improper test of the file type. It now correctly reads the files with an IXY format, e.g., the spectral amplitude files.

• fmtp has a flag -X  to output nodal plane strike dip and rakes to the standard output with a line that reads

SDRPLANES stk1 dip1 rake1 stk2 dip2 rake2
The reason is that I wanted to examine the dip variation in the L'Aquila aftershock sequence

• elocate was changed so that the location time is also written in a format for use by the gsac command ch ocal, e.g., ch ocal 2009 07 29 10 00 37 123. The new output is indicated in the color red:


RMS Error        :               0.150              sec
Travel_Time_Table:          nnCIA
Latitude         :             35.9288 +-    0.0130 N         1.4501 km
Longitude        :             43.2796 +-    0.0173 E         1.5470 km
Depth            :               13.65 +-      0.88 km
Epoch Time       :      1247949147.952 +-      0.31 sec
Event Time       :  20090718203227.952 +-      0.31 sec
Event (OCAL)     :  2009 07 18 20 32 27 952
HYPO71 Quality   :                  CD
Gap              :                 161              deg

In addition two more prototype velocity models were added. The KOREA model corresponding to my pre-computed Green's functions for Korea, and the nnCIA model corresponding to my pre-computed Green's functions for waveform of earthquakes in the Central Italian Appenines.

#### July 22, 2009

• Corrected an inconvenience in the gsac write command.  The program was in the specification of the APPEND suffix and PREPEND prefix option. Until this modification the shortcuts AP suffix and PR prefix were permitted. The string matching in gsac_arg.c only considered the first two characters. This an attempt to execute the command  GSAC> write APZNZ10HHR APZNZ10HHT APZNZ10HHZ was doomed since the first AP was interpreted as APPEND rather than a station name.  This is fied now, but NEVER use APPEND or PREPEND as a file name in a write command.  Implementing this feature of sac force the use of a particulatly bad command grammar. [Thanks to Harley Benz]

• Documentation for the Inversion of Surface-Wave Dispersion and Receiver Functions for Structure was corrected to show the correct use of the -ZMAX flag for shwmod96 [thanks to Wen-Jie Jiao]

#### June 19, 2009

• Jessie Bonner noted some compile problems:

• Change fortran code that had type declaration after a data statement: VOLIV/src/jsamat.f (line 45) used in joint96 and gprep96 (line 92)

• Modify subroutine four. Evidently a g77 compiler, which is no longer supported, did not like the fact the there was a call four(complex_array …) while internally the code was subroutine four(real_array … ). This was valid in FORTRAN77 since the memory storage would be the same if one had a real array x(20), a two dimensional array x(2,10) or a complex array x(10) since the memory organization was:

         MEMORY          |     1    |     2    |     3    |     4    | ...
1-D real array  |   x(1)   |   x(2)   |   x(2)   |   x(4)   | ...
2-D real array  |  x(1,1)  |  x(2,1)  |  x(1,2)  |  x(2,2)  | ...
Complex array   | Re[x(1)] | Im[x(1)] | Re[x(2)] | Im[x(2)] | …

This was carefully fixed by defining a revised routine zfour(z,npts,isign,dt,df) in all code. This is a step in defining every variable in a manner compatible with modern computer languages.

The program files  updated were as follow:
VOLII/src/sacmat96.f
VOLII/src/sacmft96.f
VOLII/src/sacpom96.f
VOLIII/src/spulse96.f
VOLIV/src/rftndr96.f
VOLIV/src/twofft.f
VOLV/src/ffilt96.f
VOLV/src/fspec96.f
VOLVI/src/hpulse96.f
VOLVI/src/hrftn96.f
VOLVI/src/hudson96.f
VOLVI/src/trftn96.f
VOLVIII/src/sacdecon.f
VOLVIII/src/sacevalr.f
VOLVIII/src/sacfilt.f
VOLVIII/src/saciterd.f
VOLVIII/src/sacspc96.f
VOLIX/src/cpulse96.f

The following documentation was updated:
DOC/OVERVIEW.pdf/cps330o.pdf
DOC/GSAC.pdf/cps330g.pdf
DOC/STRUCT.pdf/cps330c.pdf
DOC/SOURCE.pdf/cps330s.pdf
CALPLOT.pdf/cps330p.pdf

#### June 11, 2009

• Hermann Zeyen noticed a problem in the FORTRAN  implementation of the subroutine arsac in the file SUBS/sacsubf.f.  "if the number of data is multipleof five, the program tries to read one line too much."  This is not a problem in the C version in SUBS/sacsubc.c

• Modified the SLU version of rdseed to improve comments in the sac pole-zero file.

#### May 23, 2009

• Corrected moment magnitude definition in fmech96 to replace the 16.05 by 16.1 in log Mo (dyne-cm) = 16.1 + 1.5 Mw

• Added a new command to gsac. The mt command uses the precompute Green's functions (e.g., f96tosac -B -G or f96tosac -B -T ) to compute a Z, R, T  seismogram or the 3 ZRT seismograms for a doubel couple source. The output will be in units of m/s for a step dislocation.  The purpose of this command is to be able to quickly create forward synthetics for use in finite-fault inversion.

• Changed the part sacmft96 that selects the period limits.  The previous code looked like if(p .ge. pmin .and. p .le. pmax.and. p.ge.permin .and. p.le.permax)then which had roundoff problems as noted by Renata Schivardi at INGV Bologna. The code  was changed to be slightly more generous by changing the test to           if(p .ge. 0.99*pmin .and. p .le. 1.01*pmax .and. p.ge. 0.99*permin .and. p.le.1.01*permax)then

• Changed the part of sacmat96 that selects the period limits. The previous code was if(p .ge. pmin .and. p .le. pmax)then This was changed to if(p .ge. 0.99*pmin .and. p .le. 1.01*pmax)then to overcome roundoff error. This problem was noted by Renata Schivardi at INGV Bologna.

• Modified gsac to permit a title to be placed on the figure generated by the plotsp command.

#### March 30, 2009

• Corrected xlim command of gsac so that the plot axes are correctly reset after the command xlim off.

• Added new options to the gsac rotate3 command:  rotate3 to uvwsts2 and rotate3 to uvwtril.

• The reason for this is that a field system had a broadband sensor with internal UVW sensors and also an accelerometer. As part of seismic network quality control, the recordings of each were corrected for instrument response and the derived ground motions were compared for that part of the signal with high signal-to-noise. These should be identical since the ground motion is the same.  For our field system, this did  not occur, and we returned the broadband and accelerometer to the university and compared the recordings to a reference system.  A similar comparison showed that the problem was with the broadband sensor. The rotate3 to uvw was added to gsac to isolate the internal UVW sensor that was the cause of the problem. To do this the ZNE motion of the test and reference systems were rotated to the UVW coordinate system of the test system and the individual UVW traces were compared. Because we compared two different 3-component sensors, recording of the UVW channels of each would not have worked since each was different.

A description of the coordinate system transformation inplemented now is gsac given in the PDF document rotate3.pdf

#### March 19,2009

• At the request of Chuck Ammon, hudson96 now has a flag that permits one to override the model computed T*. This is invoked with the command option -TSTAR tstar .  A valid range of the tstar is >= 0.  If one enters a value of 0.0 the program sets this internally to 0.0001 as a quick way to implement this without a major code rewrite.

#### March 18, 2009

• Incorporated a fix in subroutines PROGRAMS.330/SUBS/sacsubf.f and PROGRAMS.330/SUBS/sacsubc.c to use the Enumerated Types for IZTYPE of IUNKN, IB, IDAY, IO, IA, IT0, IT1, ..., IT9.

• Modified PROGRAMS.330/VOLVI/src/hudson96.f to fix logical errors in the time delay correction for the T* operator, and more importantly to overcome numerical proplems for the S-wave incident at the receiver. This was done by replacing an explicit matrix sub-determinant by a compound matrix value.

• Modified PROGRAMS.330/VOLV/src/time96.f to compute the travel time and ray parameter for the sS arrival.

• The distribution now includes the source code for the menu program for terminals, dialog (dialog-1.1-20080819)

#### March 3, 2009

• Modified the file gsac_sac.c in PROGRAMS.330/VOLVIII/gsac.src  to improve the detection of a non-sac file.  Hopefully this will prevent problems when trying to read sac files with the incorrect byte order.

#### February 21, 2009

• Renamed the utility program CAL in PROGRAMS.330/CALPLOT/Utility to calplt.  The reason for this is that MacOS-X does not distinguish between upper and lowercase executable names, and CAL conflicts with the calendar program cal.  Since CAL was used in shell scripts to annotate plots, instances of CAL will be changed to calplt until all are corrected.  The graphics output of the program is renamed CALPLT.PLT to be consistent with the naming convention of other programs. Correct documentation in PROGRAMS.330/DOC and MOMENT_TENSOR tutorial and others tutorials to reflect this change.

• Corrected the MOMENT_TENSOR tutorial scripts to provide more information on copying files. Also corrected the DOBYTEORDER in the GREEN directory of precomputed Green's functions.

• #### February 18, 2009

• Fixed a but in the programs spulse96, hspec96, hudson96, time96, timmod96 and refmod96  which was due to the fact that the parameter refdep was not defined in subroutine frstar.  This was discovered because of memory dump problems on some computers.

#### January 5, 2009

• Introduced a new program fmlpr which is based on the program lprmech written by Charles J. Ammon.  This program is controlled from the command line and creates two files. The first is a printer-plot tabulation of a double-couple focal mechanism solution, and the second provides those results in a global CMT 'ndk' format. A example run if given by the line

fmlpr "TITLE" "DATETIME" "GEOG" EVDP MW strike dip rake "FILE1" "FILE2" "FILE3"

TITLE DATETIME and GEOG are enclosed in quotes
FILE1 - name of file with stations used
FILE2 - file with data summary: NSTA NCMP MINPER
FILE3 - file with filtering commands

fmlpr "USGS/SLU Moment Tensor Solution" "${YEAR}${MO} ${DY}${HR} ${MN}${SEC} ${MSEC}${LAT} ${LON}${DEPTH} ${MAG}" "${STATE}" ${HS}${MW} ${STK}${DIP} ${RAKE} sumstaused.txt sumswavendk.txt sumwfilter.txt Here the first string provides a title for the printer plot, the second string provides information about the earthquake location and the third string provides a geographical description of the earthquake sourcce zone. These nest items, e.g., HS, MW, STK, DIP and RAKE are the source depth, moment magnitude and double-couple source parameters. The last three items are text files for turther annotation. The sumstaused.txt lists the stations used for the source inversion, sumswavendk.txt consists of a single line with three numbers - the number of stations used, the number of components and the minimum period used in the inversion. The last file, sumwfilter.txt, lists the gsac signal preparation commands so that the filter band is known. When invoked as fmlpr "USGS/SLU Moment Tensor Solution" "2008 10 14 03 07 28 000 35.76 -100.704 5.0 3.80" "Texas" 11.0 3.65 285 70 -80 sumstaused.txt sumswavendk.txt sumwfilter.txt with sumstaused.txt containing IU.ANMO TA.MSTX US.CBKS US.JCT US.KSU1 US.OGNE US.SDCO US.WMOK and sumswavendk.txt contining 8 19 10 and sumwfilter.txt containing hp c 0.05 n 3 lp c 0.10 n 3 br c 0.12 0.25 n 4 p 2 two files are created: 20081014030728.msg, which is the printer plot, and 20081014030728.ndk, which is the solution information in the Global CMT ndk format (http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/allorder.ndk_explained ). The contents of the file 20081014030728.ndk are ENS 2008/10/14 03:07:28:0 35.76 -100.70 5.0 3.8 0.0 EL SALVADOR C20081014030728A B: 0 0 0 S: 8 19 10 M: 0 0 0 CMT: 1 TRIHD: 0.2 CENTROID: 0.0 0.0 35.76 0.00 -100.70 0.00 11.0 0.0 FREE S-20090105234025 21 -2.379 0.000 2.526 0.000 -0.147 0.000 2.681 0.000 -0.949 0.000 -0.064 0.000 V10 3.758 24 7 0.000 9 102 -3.758 64 211 3.758 285 70 -80 78 22 -116 #### December 22, 2008 • Modified two routines in gsac: gsac_rot.c and gsac_rot3.c. Basically I had problems with statements like if( az[0] != az[1] ) and if ( baz[0] != baz[1] ) because of roundoff error. So I repalced this with if(ABS(az[0] - az[1]) > 0.01) and the same for the baz test #### November 12. 2008 • Modified plotxvig to permit the use of a -C flag. When this is set, the arrow cursor for display is replaced by a crosshair. The pupose of this is to be able to read axis numbers for a point within a plot. • Modified do_mft to permit the use of a cross-hair cursor. Note that this is not perfectly implemented. The purpose of both of these changes is to be able to compare the dispersion shown in a regular plot to that of the multiple filter analysis when plotting both side-by-side. I was comparing observed dispersion to an independent tomographic estimate. #### September 3, 2008 • Modified udelaz so that one can also use -GCARC to get the epicentral distance in degrees and -DIST to get the epicentral distance in kilometers. The previoous -DELDEG and -DELKM still work. This was done to make the command syntax closer to that of gsac #### August 4, 2008 • Modified the CALPLOT graphics programs plotxvig, plot4014, plotnps, plotgif and plotpng to support -X0x_offset and -Y0y_offset flags. These were always part of the reframe command for repositioning a plot. I added this feature for use mostly with plotxvig so that I can display plots that are too high on the page. • Changes #define NCHAR 1000 in VOLVIII/gsac.src/gsac_parse_command.c to #define NCHAR 2000 to permit more file names on the command line. I niticed this when trying to read in a lot of files, but I dod not have this problem with using a wild card. #### June 13, 2008 • Corrected the gsac_refr.c (refraction analysis routine) of gsac, so that the refraction velocities are given correctly for a reverse spread, e.g., negative distances which would yield negative slopes and velocities. • Modified refmod96 so that the line width parameter (e.g., command line -W width) applies to the reflected arrivals. For multilayered models, the line width or pattern is a function of the reflection coefficient of the deepest reflecting interface to be able to focus on the important reflection. The reflection times are still the times to the deepest layer, e.g., if the deepest layer boundary is 3, then the rays computed are 1 1, 1 2 2 1 for the first reflection, and 1 1 1 1 and 1 2 2 1 1 2 2 1 for the first multiple; however the mixed multiple 1 2 2 1 1 1 is not computed. A plot showing the refraction and reflection predictions is given in refmod96 examples #### May 25-29, 2008 • Corrected hudson96 so that the SH arrival times are computed and ultimately placed in the file96 and sac file trace headers (May 25, 2008). To properly set the arrival time, the start of the teleseismic attenuation pulse is determined empirically using two FFT's. This must be done using the same DT as for the synthetic because of the frequency dispersion inherent in the causal T* operator (May 28, 2008). Corrected indexing when merging source and receiver structure with teleseism structure. The halfspace was not inserted correctly (May 29, 2008). #### May 7, 2008 • Modified evalresp in PROGRAMS.330/IRIS/evalresp-3.2.38.slu to handle FIR stages with a large number of coefficients. Specifically, Change Eric J. Haug SLU May, 2008 to handle continuation for FIRST with a large number of coefficients. Basically we need to get past the tendency of some SEED writers to put the following phase in the RESP file "No Abbreviation Referenced" for Response in units lookup and Response out units lookup. The modified code provides the same output as the current version of the IRIS jevalresp #### April 9, 2008 • Corrected error in gsac (gsac_sort.c) which prevented a sort on an integer variable. Also correct a typo so that a sort up ihdr11 works. #### March 31, 2008 • Created a new program hudson96 which makes approximate teleseismic synthetics following the work of Hudson(1969). This approach focuses on the direct teleseismic P or S arrival through the use of a stationary phase approximation for the ray parameter connecting the direct arrival from the source to the receiver. The interesting feature of this approach is the use of propagator matrices for the source and reciever crusts, so that the depth pahses are automatically included as well as the details fo the source and receiver crustal models. The program is self contained in that the travel time, T* and geometrical spreading are comptued within the program. A tutorial on the use of the program is on the tutorial page. Hudson, J. A. (1969). A quantitative evaluation of seismic signals at teleseismic distances - II. Body waves and surface waves from an extended source, Geophys. J. Roy. astr. Soc. 18, 353-370. • The programs time96, timmod86, refmod96, spulse96, and hspec96 were upgraded with the same travel time routine used in hudson96, to now compute the geometrical spreading for teleseismic arrivals. #### March 14, 2008 • Added a new flag to hprep96. The -ML length_mul is typically used in increase the length for wavenumber integration,L , estimated using the Bouchon criteria, which is that dk = 2 pi / L . There are circumstances when the default value of L is a little too small, and significant noise signals occur at later time. This option is easieer than trying to guess the correct value of L for the command line override, e.g., -XL xleng #### February 8, 2008 • Corrected an error in the prs command of gsac. Because of an incorrect clipping command in the graphics, the trace annotation was not done correctly when using the prs ann sta command. • Tested the Earth flattening used in hspec96, tspec96, surf96, joint96, sdisp96, sregn96 and slegn96 #### January 28, 2008 • Created a new program ttime96 • Installed identical travel time routines for TI models in the following programs: ttime96, ttimmod96, tspec96 and hprep96. • The travel time validation is given at the links Isotropic Model Test and Transverse Isotropic Model Test • Fixes errors in tspec96 that affected synthetics with Q. For an isotropic model, the output of tspec96 agrees with that of hspec96 #### January 22, 2008 • Installed identical travel time routines for isotropic models in the following programs: time96, timmod86, refmod96, spulse96, hspec96. The new routines provide the correct results for flat and spherical velocity models for travel times, ray parameter, T* of the P, S, pP and sP phases. They also correctly handle the artifice of representing the atmosphere by a stack of layers with negative layer thickness so that a depth of 0 corresponds to the surface of the solid. #### January 17, 2008 • Modified gsac. Under ppk / plotpk the 'l' or 'L' command returns the epoch time to the right of the amplitude. The purpose of this is to facility estimating time differences from two successive 'L' commands on different parts of the trace #### January 10, 2008 • Corrected a memory allocation error in PROGRAMS.330/VOLII/src/do_pom3.c #### January 3, 2008 • There were a number of changes made to the wavenumber integration and surface-wave codes to permit a better Earth-flattening transformation to permit the computation of seismigrams for a spherical Earth by using the plane-layered code. These transformations were tested using mineos. A description of the current status of the Earth flattening studies is given at the following link. http://www.eas.slu.edu/TUTORIAL/SPHERICITY A detailed tutorial on teleseismic synthetics and source inversion must be prepared. This update changes the following programs: sdisp96, sregn96 and slegn96 in PROGRAMS.330/VOLIII/src, srfdis96, srfdrr96 and srfdrl96 in PROGRAMS.330/VOLIV/src, and hspec96 in PROGRAMS.330/VOLVI/src. Note hspec96 was modified on January 14, 2008 to correctly handle the source and receiver depth insertion for the flattening. The use fo the reference depth still needs work to handle atmospheric layers about the surface. #### January 2, 2008 • #### Created a new program: redodateThe purpose of redodateis to change a data/time specification through an offset, e.g., redodate 2009 01 01 00 01 02 456 -3660 returns 2008 12 31 23 00 02 456 I require this in order to specify window segments for traces that I wish to download. Previously i started at the origin time, but I needed more signal for inversion, so I would like to start some offset (a negative number) before the specified origin time. This code segment indicates how this would be used in a shell script, where I would like to start 60 seconds before the origin time: YEAR=2009 MONTH=01 DAY=01 HOUR=00 MINUTE=01 SECOND=02 MSEC=456 OFFSET=-60 redodate${YEAR} ${MONTH}${DAY} ${HOUR}${MINUTE} ${SECOND}${MSEC} ${OFFSET} > tmp read YEAR MONTH DAY HOUR MINUTE SECOND MSEC < tmp echo "YEAR :"${YEAR}
echo "MONTH    :" ${MONTH} echo "DAY :"${DAY}
echo "HOUR     :" ${HOUR} echo "MINUTE :"${MINUTE}
echo "SECOND   :" ${SECOND} echo "MSEC :"${MSEC}
rm -f tmp

#### November 17, 2007

• Upgraded the distributed IRIS programs evalresp from version 3.2.25 to 3.2.38 and rdseed from 4.5 to 4.7.5.  In order to support the CPS architectures, e.g., CYGWIN (Windows), MacOS-X, LINUX and Solaris, changes were made to these codes. The changes are given in the links in this section. This upgrade was required in order to continue to use the most recent SEED volumes.

The two programs are installed in PROGRAMS.330/bin with the executable nanames evalresp and rdseed. The IRIS distribution of rdseed provides the executables rdseed.linux, rdseed.ppc-mac, rdseed.intelmac and rdseed.solaris. These binaries are available in the PROGRAMS.330/IRIS/rdseed4.7.5 directory of CPS. If you wish to use these programs, you must install them in your system and adjust the path to find them.

• evalresp  - The distributed code has an #include <malloc.h>. This is deprecated in MacOS-X to <malloc/malloc.h> which means that the code will not compile.  Since modern compilers provide the prototypes for malloc in <stdlib.h>. I changed alloc_fctns.c, parse_fctns.c, regexp.c and string_fctns.c by commenting out the #include <malloc.h> and ensuring the presence of a #include <stdlib.h>.
Detailed changes.

NOTE THAT evalresp-3.2.40 will configure properly for the MacOS. When this is officially released, it will be included in CPS!

• rdseed - Several changes were made.

• The most important was to change the use of long in the definition of the integer SAC header values to int. Becore computers were truly 32 bit, short meant a 2 byte integer, long a 4 byte integer, and int could mean a  byte integer. Currently, int means 4 bytes.  A problem arises on 64 bit machines, for which a long means an 8 byte integer.  One way around this is to use the -m32 flag on the compiler.  CPS changes the long to int.

• Another change relates to the strupr routine. This is actually provided by some libraries.  Thus all instances of strupr are changed to Strupr and a definition is placed in Include/rdseed.h

• Utilities/alloc_linklist_element.c - add #include <stdlib.h>  to define malloc()

• Ah/makefile - remove CC=gcc  since not all compilers use gcc - note that  CC is defined in the main Makefile

• Decoders/makefile - remove CC=gcc  since not all compilers use gcc

• Parsers/makefile - remove CC=gcc  since not all compilers use gcc

• Utilities/log_errors.c - change the syntax of stderr_to_file from stderr_to_file(char *fmt, char*ap) to void stderr_to_file(char *fmt, va_list ap) for comparibility with <stdio.h>

• Main/output_sac.c - Place comments into the pole-zero file giving information about the station, network, date, etc.
(Sample SAC pole-zero file)

• Main/output_resp.c - Send information about the network to stderr, e.g., station, network, component, location, on_date, off_date, latitude, longitude, elevation, response file name, component azimuth and dip and sample rate in print_resp.
(Sample stderr output)

• Updated the manual entry in Main/rdseed.1 to make the command line flags more complete on the first page. The pdf version of this document is in PROGRAMS.330/DOC/IRIS.pdf where you will find evalresp.pdf, rdseed.pdf and SEEDManual_V2.4.pdf.
Detailed changes
The modified rdseed was successfully compiled on both 32 and 64 bit Intel LINUX machines using the GNU compilers, on a SPARC Solaris machine using SUN compilers, on an Intel Windows machine under the CYGWIN environment, and under MacOX-X using the GNU compilers. The only compiler warnings were the following:

cc -I../Include -c -g   -c -o decode_steim.o decode_steim.c
decode_steim.c: In function decode_steim':
decode_steim.c:231: warning: integer overflow in expression
decode_steim.c:231: warning: integer overflow in expression

for all compiles, and

cc -I../Include -c -g   -c -o parse_30.o parse_30.c
parse_30.c: In function parse_type30':
parse_30.c:96: warning: cast to pointer from integer of different size

on just the 64-bit LINUX compiler (note that the -w32 compiler flag was not used).

#### November 9, 2007

• Modified the code of wvfmch96 so that there is no internal limit on the number of waveforms for which forward synthetics are required. In the present distribution, the maximum number of waveforms permitted for use by wvfmch96, wvfgrd96, wvfmt96 and wvfmtd96 is 2000, with each trace permitted to consist of 4096 points.  Because wvfmch96 creates a forward synthetic for a given moment tensor description, this program can be used as part of a QC process to evaluate all recorded waveforms in a data center.  Since there can easily be more than 2000 such waveforms today, hard dimension limits are a major inconvenience.  The logic of wvfmch96 was modified so that the forward synthetic is created separately for each input file, rather than reading in all files and corresponding Green's functions into memory.

#### October 21. 2007

• Hacked a fix for the distance computation used in gsac (gsac_subs.c) and udelaz.  The problem was that if the event - station latitude/longitude pair differed in coordinates by 0.0001 degree or less, then the computed distance in kilometers was incorrect because of roundoff error of the angles used for the elliptic integrals. The artifice used for distances less than 1.0 km, is to use the spherical relation  delkm = 6371.003 * delradian  for the short distances.

#### October 19. 2007

• Error in the DIF routine found by Harley Benz in gsac_dif.c. The erros was that the first point was incorrectly defined.

October 04, 2007

• Modified sacmft96 and sacpom96 to take ABS(DIST) where DIST is from the SAC header.  This is because I sometimes require negative distances when using shallow refraction data.  For a split spread, the  gsac refr command will display the split spread.  The refl command will also work now.  The only bad thing is that the refraction velocities may now appear to be negative.  These changes were made as part of an effort to develop a procedure for determining shallow S-wave velocity structure.

A problem with the distance/azimuth computations in gsac  was noticed in the process of computing epicentral distances less than 1 km.  This problem will also affect udelaz. This is currently under investigation.

#### August 30, 2007

• gsac  - implemented two new commands in prs (plotrecordsection): ScaleRelative (SR) and ScaleAbsolute power (SA Power) to control the amplitude of the traces in the record section.  ScaleRelative is the old behavior which just scales each trace to the same maximum amplitude (AMP value), except that the scaling now applies only to the actual trace window and not the whole trace.  ScaleAbsolute 0.0 adjusts all traces to plot true amplitude. If the non-time axis is distance, then ScaleAbsolute 1.0  adjusts for geometrical spreading by multiplying all traces by distpower . For a true refracted arrival (e.g., with no velocity gradients), power=1.5 would show the refractions at the same level. The current  gsac version is  now[V1.1.21 30 AUG 2007]

• wvfgrd96 is now modified so that the output traces, e.g., TRACE.obs and TRACE.pre have the time shift applied. This means that if pltsac is used with the flag -USER9 the time shift will be displayed even though the traces align. The objective here is to present the observed and predicted in such a way to focus on amplitude and phase differences and not time shifts.

#### August 19, 2007

• gsac  - modified the write command to support 'append suffix' and 'prepend prefix'.  To illustrate the differences between sac2000, we will work in a subdirectory with the files located in the level above:

 > sac2000  SEISMIC ANALYSIS CODE [8/8/2001 (Version 00.59.44)]  Copyright 1995 Regents of the University of California SAC> r ../*.sac ../imp.sac ...nimp.sac SAC> w append .s ../imp.sac.s ...nimp.sac.s SAC> ls SAC> r ../imp.sac ...nimp.sac SAC> w prepend p. p.../imp.sac ...nimp.sac  ERROR  101: opening file SAC> GSAC - Computer Programs in Seismology [V1.1.19 01 AUG 2007]        Copyright 2004, 2005, 2006, 2007 R. B. Herrmann GSAC> r ../*.sac ../imp.sac ../nimp.sac GSAC> w append .s imp.sac.s nimp.sac.s   [Note write in local directory not original ] GSAC> r ../imp.sac ../nimp.sac GSAC> w prepend p.    [Note write in local directory not original - note Sac2000 ] p.imp.sac p.nimp.sac      [cannot prepend in parent directory ] GSAC> w append .S prepend P.   [Note write in local directory not original ] P.imp.sac.S P.nimp.sac.S GSAC> ls                      [show what was created in the current directory] imp.sac.s  nimp.sac.s  p.imp.sac  P.imp.sac.S  p.nimp.sac  P.nimp.sac.S

#### August 16, 2007

• gsac modified to permit the use of the SAC commands plot1 (p1) and plot2 (p2) to force multi-trace plot (which is the default of the gsac plot ) and overlay (respectively). Note that the SAC plot shows a single trace at a time.  This is inconvenient for normal use and gsac will reserve the plot (p) for a multitrace view.

#### August 14, 2007

• Initial development of new source model for surface-wave and waveform inversion - the crack. Following Nakano, M., and H. Kumagi (2005). Waveform inversion of volcano-seismic signals assuming possible source geometries, Geophys. Res. Letters 32, L12302, doi:10.1029/2005GL022666, a new command line option -CRACK is permitted for the programs srfgrd96, sdprad96, sdpspc96 and wvfmch96. The crack is specified in terms of the parameters -DIP dip -STK -RAKE rake -MW mw or -M0 mo where dip is the dip of the crack,  stk is the dip-direction of the crack (dips down to the right when looking in the direction of strike).  The rake is used to indicate an expanding crack (positive) or closing crack (negative).  The  moment magnitude or moment are positive numbers.

At present the program fmplot is not modified, nor is the output of sdprad96 to indicate that the crack model was used instead of the earthquake double couple.

The reason for this introduction of a new source was to study the Utah mine collapse of 2007 08 06 08 48 40!

#### August 1, 2007

• gsac - corrected map command so that it now correctly handles the station off and epicenter off commands. Previously the epicenters were always plotted.

• gsac is at version V1.1.19  01 AUG 2007

#### July 16, 2007

• gsac - corrected error in the remove trend rtr command. Roundoff error messed up the determination of the regression equation a + bx . Previously x was the time index, e.g., i=0, i < npts. For long time seeries this let to large errors in the determinant. I replaced this by defining x as i/npts with varies between 0 and 1 for the time series.

#### July 1, 2007

• gsac- implemented a new option to ppk (plotpk), e.g.,

ppk pq    or      ppk pquality

This is functionally similar to the previously introduced  'ppk quality'  except that the purpose is to simultaneously define a useful trace and to reset the P-time. For my source inversion, I set the theoretical P-time using the AK135 predicted times.  This is OK for source inversion since the inversion code permits a simple waveform time shift. However for broadband depth or energy determination, it is necessary to start at the P time.  To speed up the review, I added this new option. So now move the crosshair to the P-wave first arrival, and then click any mouse button. Simultaneously the P-arrival time is reset and the trace is marked as one to be used (IHDR20 is set to +1).

#### April 23, 2007

• saciterd - repaired an error with the use of the -2 flag. This was not done correctly in the earlier code. A change was made in late 2006 - early 2007 on line 490. This was not correctly done. This was repaired today. Running the code with the tutorial data set yields the same figures as in the tutorial.  Matt Sibol of Multi max.com noted the discrepancy.

#### March 31, 2007

• saclhdr - modified the format statement for printing reals from %f to %g so that a band amplitude because of dividing by zero gives +1e+38 and -1e38 for DEPMAX and DEPMIN

#### March 19, 2007

• gsac - Added an option to ppk (plotpk), e.g.,
ppk q
For source inversion, one must examine many traces to decide which to use since the inversion results are strongly influenced by bad data.  For teleseismic and regional data, one may review hundreds of traces. So gsac was modified to permit a rapid assessment of the trace quality. The assessment is placed in the integer header value IHDR20 (gsac notation) and a subsequent shell script examines the value there to decide if the trace should be moved to the final directory for source inversion.

The purpose is to use the plotpk display to quickly define whether  to use a trace for a source inversion. For teleseismic inversion one must examine the quality of the P-wave arrivals, so one reads the traces, bandpasses, sorts with distance, applies a suitable xlim and executes the command
GSAC> ppk   relative perplot 4 q
The result is the image

The default behavior is that a mouse click says accept - of course one can use the 'a' key to accept and the 'r' key to reject. The default behavior of the mouse click can be changed using the menu at the bottom. When a trace is accepted or rejected, a red '+' or a blue 'X' is drawn at the lower left corner of the trace.  The effect of these operations is to set the integer header position 20, which I call IHDR20 to +1 for accept and 0 for reject. The default value of this header field is -12345  This changes added about 300 lines of trivial code to gsac_plotpk.c. (In the image above, the P times are from the velocity model. The slight offset is no problem for the inversion since the inversion programs implement a time shift). This process of performing QC for P, SH and SV is shown the shell script IDOQC

The value of this approach to assess trace quality was immediately apparent by the resulting improved source inversions. The current version of gsac is V1.1.18 31 MAR 2007

#### March 15, 2007

• time96 - corrected the on-line help to reflect the -P -SV and -SH flags

• wvfmch96, wvfmtd96, wvfmt96 and wvfgrd96 - added error checking on the input files - reject a data set if the number of sample points is <= 4 or if DEPMAX=DEPMIN, which indicates a faulty data set. Current version of wvfmtd96  and wvfmt96  were modified to use both the external a priori  and internal variance weighting schemes.

#### February 22. 2007

• time96 – now supports the -P -SV and -SH flags to give first arrival time and corresponding ray parameter for the first arrival. Previously only supported P-wave. Why? I wanted to automate source inversion by having this program define first arrival times rather than interactively reading them. This is possible because I already know the location. The section of the shell script that does this is

for i in *sac
do
GCARC=saclhdr -GCARC $i EVDP=saclhdr -EVDP$i
A=time96 -M ${MODEL} -GCARC${GCARC} -T -EVDP ${EVDP} T0=time96 -M${MODEL} -GCARC ${GCARC} -T -EVDP${EVDP} -SV
echo ${GCARC}${EVDP} ${A}${T0}
gsac << EOF
rh $i synchronize O ch A${A} T0 ${T0} T1${T0}
wh
quit
EOF
done

The 'synchronize O' is used to reset the time stamp so that 'O' , the offset of the origin time with respect to the reference time, is zero. In this case then the 'A' will be the travel time of the P first arrival.

#### February 20. 2007

• The fft command of gsac was augmented to permit the syntax fft length 1 or fft length 2 or fft length 4 which adds additional zeros to increase the frequency domain resolution. gsac is now at version V1.1.17 22 Feb 2007

#### February 8, 2007

• The computation of epicentral distance in kilometers was changed in gsac, udelaz and edabac. Prior to this it was estimated from the great circle arc, thus assuming a spherical Earth. It is now computed numerically by direct evaluation of an incomplete elliptic integral of the second kind. Sac uses Rudoe's reverse formula (Geodesy, G. Bomford, Carlendon Press, 4th ed, 1980, Section 2.15(b), page 121.  Rather than getting confused about inverse tangents, the idea of Rudoe was used, e.g., the plane through the center of the Earth, the station and epicenter intersects the spheroid as an ellipse, and the linear distance is computed along the ellipse. However, this derivation uses vector analysis more that it does trigonometric relations.  The computation of the SAC  header value dist agrees with the values of sac2000 (based on a comparison of 19554 different combinations of EVLA, EVLO, STLA and STLO.

• Output formats of sacmft96 were changed to avoid  format overflow (***** in FORTRAN) when  AZ, BAZ and GCARC  had the undefined values of -12345.

#### January 24, 2007

• timmod96 has two new out  options:

timmod96 -TIMEP pfile
timmod96 -TIMES sfile

where pfile or sfile are ASCII files with a one distance – time pair per line, e.g.,

124.454140 20.464310
238.482346 35.048897
145.274033 24.015045
350.404083 48.721909

which can be created from the DIST and A (for P) sac header values through a loop

for i in *Z
do
saclhdr -NL -DIST -A >> ptime
done

This option permits plotting measured first arrivals on the same plot as the model predicted values

• gsac – the merge command now has the command line option GAP gap_value. Before this change, a merge of SAC files having gaps resulted in a zero fill of the gaps. Now the user value is inserted. This may permit a user program to detect the gap. For example if the SAC trace consists of integers from the ray 24 bit digitizer, the expected extreme values are ±8388608. A “MERGE GAP 1.0e+9” would fill the trace with an impossible number which can then be flagged.

• gsac – new command outcsv This command causes all traces in memory to be output in ASCII as comma separated variables (CSV) in the file name f001.csv. Each row of output consists of the values at one time sample with entries “Time,amp_trace_1,amp_trace2 ...” The reason for this command was to create a simple output format that could be read by Excel or OpenOffice spreadsheet programs so that students can perform simple time series analysis and plots. Note that OpenOffice does not permit more than 65536 to be read in

#### January 13, 2007

• Modifications to gsac: Introduced new commands grid, xgrid, ygrid, background. The grid options permit setting a grid on plots and also too define the color and type of the grid lines. The background command permits setting a background color in the plotting frames.

• GSAC>
GSAC - Computer Programs in Seismology [V1.1.15 17 JAN 2007]
Copyright 2004, 2005, 2006, 2007 R. B. Herrmann

Also corrected appearance of the help screen and got the grids to work on the plotsp command. My excuse is that I only had battery power in a house without electricity for two days because of the ices storm

#### January 3, 2007

• Fixed fmech96, which did not set the back azimuth correctly. Now one can create Z R T and Z N E three-component synthetics using fmech96, convert them to SAC format using f96tosac and used gsac to read in the Z N E to make the Z R T and everything agrees.

#### January 3, 2007

• Correct error in surf96, joint96. Love only dispersion data were not processed correctly. The minor change was in surf96.f and joint96.f

#### January 3, 2007

• Changed all dimensions for layering to now have 200 layers instead of 100. Note this does not apply to the Cerveny programs cprep96, cpulse96, cseis96.

• changed time96 so that it now can output the ray parameter. This gets around the problem of using udtdd which is an implementation of the Jeffreys-Bullen curves. There is now a consistency between the model96 synthetics and the travel time prediction. time96 is used to predict ray parameter and first arrival for teleseismic P-wave synthetics. See the tutorial.

#### December 28, 2006

• Reformatted help entries in gsac to be easier to read.

• Implemented new commands in gsac: version, triangle, boxcar, trapezoid. version lists the version number. triangle, boxcar and trapezoid are time domain filters that convolve the traces in memory with an isosceles triangle defined by its halfwidth, a boxcar defined by its width, and a trapezoid defined by the length of its three segments.  All of these pulses are designed to have unit area and function as a lowpass filter. There is some granularity in the implementation that means that the sequence of commands

triangle width 1.0

is not the same as two commands

boxcar width 1.0
boxcar width 1.0

even though this should be true by virtue of the convolution theorem.

#### December 15, 2006

• Changed SUBS/f96subf.f and recompiled everything using file96 format because of incorrect format for writes. This now fully implements the change of December 12, 2006.

#### December 12, 2006

• Changed VOLIII/src/spulse96.f VOLV/src/genray96.f VOLVI/src/hpulse96.f VOLIX/src/cpulse96.f to set the file96 header values for evlat evlon stlat and stlon to -12345 in order to set up the ultimate SAC header correctly.

#### November 8, 2006

• sdpdsp96: Corrected error in test of xmin/xmax/ymin/ymax that precluded a plot

#### October 18, 2006

• Changed map.sh created by gsac map command to remove references to RESIN which is already set by the GMTXXX/share/dbase/grdraster.info. The latest version of gsac is [V1.1.12 19 OCT 2006]

#### October 14, 2006

• fixed obtuse error in sregn96 slegn96 sdisp96 for an all fluid velocity model. this did not affect solid or fluid/solid model computations. Also fixed the depth axis notation for sdpder96.:

#### October 1, 2006

• Fixed allocation error in do_mft3.c of do_mft

#### September 16, 2006

• gsac Version 1.1.11 – Commands ignore a leading #. Also the pole-zero and frequency-amplitude-phase files ignore leading blanks. Finally the ABS macro in gsac.h was correctly written, meaning that the x -x works properly in plotpk whether the first x is the earliest or the latest time. The ignore leading blanks in the response files was required to overcome some laxness in some pole-zero files. Finally the ignore # serves the same purpose as in the sh or bash shells – e.g., to have an inline comment. This means that gsac commands in a shell script can now have comments. For example, we can now

     #!/bin/sh
#####
#    generate a pulse and low pass filter it
#####
gsac << EOF
#####
#        generate a pulse
#####
fg triangle delta 0.05 npts 1024
            #####
                  #        lowpass at 1 Hz
#####
lp c 1 n 2
#####
#         save as the file LP.sac
#####
w LP .sac  # a comment can be placed here too
quit
EOF  

#### August 20, 2006

• Modified the CALPLOT routines plotnps and plotdbg to work in a 64-bit environment by removing and reference to the old MS-DOS. The problem was that negative numbers were incorrectly written with the "%ld" format specification. These were replaced globally by "%d".

#### August 8, 2006

• Corrected a major error in the first arrival computations in gsac in the interpolation routine. This was caused by a time series with more than 10,000,000 and involved floating point roundoff problems in the time array used for interpolation.

#### August 8, 2006

• Began adapting codes to the gcc/gfortran (4.2) compiler suite. For the most part this just required adhering to the strict 72 character line for the FORTRAN 77 standard.

• The gfortran 4.2 uses buffered IO, thus overcoming slow IO of the gfortran 4.0 suite

• The new makefiles force a 4 byte header for the unformatted sequential IO instead of the default 8 byte.

#### August 4, 2006

• Corrected a major error in the first arrival computations in time96, timmod96, refmod96, hprep96, hspec96, tspec96, spulse96. Faulty logic gave too early an arrival time if the first arrival was the direct wave. This was fixed by computing the distance at which the critical ray appears; if this is greater than the desired distance, then use the direct arrival time. In addition we found that one of the refraction times was not comptued correctly which velocity decreases with depth. The error means that the pole-zero option was not correctly implemented. (Thanks to C. J. Ammon).

#### June 16, 2006

• Correct a major error in the transfer command in gsac_trans.c in gsac. The error means that the pole-zero option was not correctly implemented. (Thanks to C. J. Ammon).

#### June 10, 2006

• Added map command to gsac. This command creates a well commented GMT script that plots stations, epicenters, raypaths and topography. The shell script can be easily edited to change the default plotting parameters.

#### May 27, 2006

• Corrected login in gsac_in.c for interpolation. The existing routine lost points.

• Added decimation command to gsac

#### April 26, 2006

• Corrected gsac_lh.c and gsac_ch.c to properly handle the IEVREG header variable. sac2000 actually has internal list of the Flinn-Engdahl regions that print out a nice string for the integer variable. gsac now just outputs the integer value of this field.

• prfmod96 is available. This was written to display regional velocity models from the joint inversion of receiver functions and surface waves in Korea. The user defines a profile in terms of the latitude and longitude at the ends of the profile. The coordinates of each model (station) and station model file are read in. This program plots a contoured velocity map from stations within +-W of the profile. The contours are labeled. The use has control over the contour intervals and the labeling intervals.

• Major changes to the VOLVII waveform inversion programs. wvfmtd96 and wvfmt96 are now multi-pass programs that perform and inversion, determine an alignment time shift and weighting, and then perform final inversion. The logic of the code follows that developed by Harley Benz at USGS for teleseismic inversion. As part of the code modifications, wvfgrd96 and wvfmch96 were also modified.

#### March 19, 2006

• Error in gsac/ gsac_prs.c that tried to flush X11 when in plot mode and not X11 mode. This caused the program to crash. Also prs relative works. I noticed this when I wants to plot observed versus predicted teleseismic P-wave traces, windowed about the P arrival as a function of distance. This was done to support the newly modified wvfmtd96 which permits time shifting

#### February 27, 2006

• Fixed major error in PROGRAMS.330/CALPLOT/src/XVIG/src/window.c The execl( ...,init_height, 0) was WRONG!. It should have been execl(...,init_height, NULL) where NUL is a (char *)0. This error appeared using a 64 bit LINUX. Thanks to Domenico Di Giacomo, INGV - Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Milano.

#### February 27, 2006

• Added a backup or return to previous display option to the plot command of gsac. This is useful when previewing many traces using the perplot option to skip back to re-examing missed features.

• Changed internal routines to prevent a mis-setting of USER1 and USER2. USER1 and USER2 are set to be the minimum and maximum periods of the filtered trace. When picking P- and S-arrival times, one may wish to filter the trace with a zero-phase bandpass filter to see the arrival. The pick is then saved using the wh command. Until this fix the USER1 and USER2 were also set. Then then led to an inconsistency that the header information would not reflect the trace content. Now the minimum and maximum periods are preserved, but only written out when the w command rewrites the trace. This fix was introduced in the source code for all bandpass filters, the transfer and whiten commands.

#### January 3, 2006

• Began adapting codes to the gcc/gfortran (4.0) compiler suite. For the most part this just required adhereing to the strict 72 character line for the FORTRAN 77 standard. There are still problems with the surf96/rftn96/joint96 suites

• Introduced an O or o option to gsac to permit synchronize to set the reference time as the origin time.

#### December 14, 2005

• Corrected the PDF documentation for the synthetic seismogram overview in PROGRAMS.330/DOC/OVERVIEW.pdf/cps330.pdf. This error was introduced in June, 2005 with the introduction of the new commands time96 and refmod96. Sorry

#### September 25, 2005

• Extended gsac to have a GUI based refraction/reflection analysis tool with the refr command

• This required modifications to the CALPLOT library by introducing a gcontrol command and hyperbolic cursor.

#### August 10, 2005

• fixed but in gsac interal time routine to correctly set milliseconds for epoch to human time

#### August 3, 2005

• Modified gsac transfer routine to handle GSE frequency-amplitude-phase FAP response file.

#### June 30, 2005

• Fixed cut option in gsac so that cut CAL YEAR MONTH DAY HOUR MINUTE SECOND MILLISECOND works correctly.

• time96 created to predict P first arrival times for flat and spherical earth velocity models. This is used to define the starting time window for teleseismic synthetics to avoid unnecessary white space before the first arrival. The use of this program in a shell script to compute the teleseismic P-wave signal but to apply a phase-velocity filter to focus on teleseismic P for speed

     for GCARC in \
20 21 22 23 24 25 26 27 28 29 \
30 31 32 33 34 35 36 37 38 39 \
40 41 42 43 44 45 46 47 48 49 \
50 51 52 53 54 55 56 57 58 59 \
60 61 62 63 64 65 66 67 68 69 \
70 71 72 73 74 75 76 77 78 79 \
80 81 82 83 84 85 86 87 88 89 \
90 91 92 93 94 95 96 97 98 99 \
100 101 102 103 104 105
do

#####
#  convert epicentral distance in degrees to kilometers
#  for the Earth
#####

dist=echo ${GCARC} | awk '{print 111.195*$1}'
#####
#       for this distance and source depth get the ray parameter
#       and define the phase velocity limits
#####
RAYP=udtdd -EVDP ${HS} -GCARC${GCARC}
CVAL=echo ${RAYP} | awk '{print 1.0/$1}'
CMAX=echo ${RAYP} | awk '{print 4.0/$1}'
C1=echo ${RAYP} | awk '{print 3.2/$1}'
C2=echo ${RAYP} | awk '{print 0.49/$1}'
CMIN=echo ${RAYP} | awk '{print 0.46/$1}'
#####
#       get P travel time == arrival time
##### A=time96 -M AK135sph.mod -GCARC ${GCARC} -T -EVDP${HS}
T0=echo $A | awk '{printf "%f",$1 - 100}' `
cat > dfile << EOF ${dist} 1.0 1024${T0} 0.0 EOF
hprep96 -M AK135sph.mod -HS ${HS} -HR 0.0 -EQEX -TF -BH -d dfile -CMAX${CMAX} -C1 ${C1} -C2${C2} -CMIN \${CMIN} -NDEC 1
hspec96 -H
hpulse96 -p -V -l 1 >> file96v
hpulse96 -p -D -l 1 >> file96d
• Modified sacmft96 and do_mft to create a temporary file called MFT96.disp instead of MFT96.DSP to avoid problems when using a caseless file system, e.g., Cygwin under Windows

• Modified sacmft96 so that color contouring ande the spectral amplitude plots reflect the envelope amplitudes within in the chosen group velocity window. Prior to this, colors and amplitudes were set by the amplitudes of the entire filtered trace.

#### June 16, 2005

• Fixed typos in gsac which had Pp instead of pP in the Teleseismic phase menu for plotpk. Also ensured that the reset in prs actually reset all parameters

• hspec96 now has an earth flattening routine for spherical earth models. Note hprep96 and tspec96 must be upgraded to handle a spherical model.

#### March 23, 2005

• Implemented B option to plotpk command of gsac to permit moving back in a multipage display. (requested by Risheng Chu of SLU)

• Implemented fileid command of gsac to permit trace labeling. The syntax was extended to permit placement at upper (UC) and lower (LC) center of each plot. (Requested by Luis Rivera and  Jeroen Ritsema)

February 11, 2005

• Added boxcar  to gsac funcgen  command - PROGRAMS.330/XOLVIII/gsac.src/gsac_fg.c

• Added debug statements to gsac transfer command to catch types in response input - PROGRAMS.330/XOLVIII/gsac.src/gsac_trans.c

• Modified PROGRAMS.330/CALPLOT/src/cmd/plotxvig.c to introduce a new window behavior. If the window is resized before an erase, the next frame will then use the window. This means that a non-default window size does nto ahve to be changed. The gframe calplot command has been changed to perform an automatic erase but to prohibit a resize.

• Corrected gsac xlim commmand to correctly handle GMT and CAL arguments and to forbid xlim o -10 o 10 if o is not set withing the header

#### February 5, 2005

• Modified PROGRAMS.330/CALPLOT/src/XVIG/src/xvig.c to unlock the keyboard for the terminal window. Successful, except that the first character is lost because of the event catching mechanism. More low level X11 work must be done. On Solaris with the CDE window manager, the trace window can only be moved use the Move under the menu option and not by dragging the window - I will not worry about this in the future

• Modified VOLIII/src/spulse96.f  VOLVI/src/hpulse96.f VOLIX/src/cpulse96.f   VOLV/src/gpulse96.f to add a -Z flag to indicate that the internal parabolic or triangular pulses are to be zero phase - this is required if we attempt to model the interstation Green's functions obtained by cross-correlation of background noise, otherwise there will be a delay causing a mismatch of synthetic and observed responses. This also necessitated a change in the documentation DOC/OVERVIEW.pdf/cps330o.pdf

#### January 31, 2005

• SUBS/SUBS

• sacsubc.c, sacsubc.h, sacsubf.c Added brsach and arsach only header information. Previously brsac was used with the overhead of reading all the data. This fixed a problem in saclhdr with very large files

• grphsubc.c fixed logic in clipping routine (baker@usgs.gov)

• VOLI/src

• stereo.f new program to plot stereonets

• VOLIII/src

• sio.f corrected minor typos

• sdprad96.f corrected minor typos

• srfgrd96.f corrected minor typos

• sdpspc96.f corrected minor typos

• VOLV/src

• genray96.f subroutine chkdep, hm1 not defined changed to dm1,  baker@usgs.gov

• gpulse96.f LER not defined in subroutine  source - baker@usgs.gov

• ttimmod96.f subroutine fstarr - vlmm undefined - changed to vlmn, baker@usgs.gov

• fmdfit.f extra long lines split

• fmplot changed default so that no nodal planes are plotted unless requested by strike,. dip, rake, force or moment tensor flags. In this simnple mode the command fmplot -F fmplot.tmp  will just plot first motion data. It is then is then up to the analyst to use the stereonet producted by stereo to determine the nodal planes and then to rerun fmplot for the final plot

• VOLVI/src

• tspec96.f vlmm not defined in subroutine fstarr relace with vlmn, baker@usgs.gov

• hprep96.f vlmm not defined in subroutine fstarr relace with vlmn, baker@usgs.gov

• VOLVII/src

• mtd.f slip not defined in subroutine trans1, baker@usgs.gov

• VOLVIII/src

• saclhdr.c use brsach instead of brsac - do not need data but want header of potentially large trace file

• saciterd.f nbumps undefined in subroutine tdomain, changed to nshift. tshift undefined in subroutine tdomain changed the theshift, baker@usgs.gov

• elocate.f emares not defined in subroutine sumup - now set to -9.99
-VELMOD flag cretes the file VEL.MOD. It no longer loists the sample velocity model fileN

• sacfilt.f LER not defined in subroutine gcmdln baker@usgs.gov

• sacdecon.f thshift undefined in main routine - change dto the shift baker@usgs.gov

• VOLVIII/gsac.src

• Many changes

• Improvements:

• plot, plotpk - low level plot routines changed to let qdp decimation be controlled by screen resolution

• plot, plotpk - low level plot routines changed to reduce number of non-drawing penup moves

• prs (plotrecordsection) - low level plot routines changed to reduce number of non-drawing penup moves

• cut, xlim - now also take GMT 2004 002 01 02 03 445  and CAL 2004 01 02 01 02 03 455 to define the time window

• GNU readline implemented - the user may use the 4 arrow keys to move around the command line and to view history of commands

• New commands

• reverse - reverse a time series

• correlate - cross-correlate time series

• sign - one bit sgn() operator

• stack - stack traces in memory

• whiten - force flat amplitude spectrum

• markt - put up velocity tics at 1, 2, 3, 4, 5, 6, 7 and 8 km/sec (note that these values are pre-assigned. The user specification is not yet implemented)

• taper - trace taper implemented

• DOC/GSAC.pdf/cps330g.pdf

• Documentation was modified to reflect new gsac commands, change to elocate, and description of fmplot (duplicated in the DOC/SOURCE.pdf/cps330s.pdf manual) and stereo

#### TO BE DONE

• add convolve to gsac

• test, test, test

• document, document, document