Enhancements Bug Fixes330

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 code

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 two links.   http://www.eas.slu.edu/Earthquake_Center/NORMAL_MODE/ and http://www.eas.slu.edu/Earthquake_Center/NORMAL_MODE/SW

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.

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 dist^power . 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:

One major difference is that although a  'w'  under both sac2000 and gsac will overwrite the original file in the upper directory, the gsac append and prepend will create files in the current directory.  Because gsac strips off the leading directory path prior to implementing the prepend command, this operation will work.

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

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.

March 19, 2006

•      #####
               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

Last changed 22 February 2007