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 and provided to users.

Enhancements to Version 3.30 codes



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



2024

March 7, 2024

February 14, 2024

February 12, 2024





2023



November 29, 2023

October 26, 2023

August 28, 2023

May 11, 2023

May 2, 2023

April 24, 2023

March 23, 2023

March 6, 2023

February 28, 2023

February 6, 2023

February 6, 2023



2022



November 7, 2022

October 12, 2022

July 29,  2022

June 9, 2022

June 6, 2022

March 20, 2022

From
float tr, ti;
...
temp = sqrt(tr*tr + ti*ti);

To
double tr, ti;
...
temp = (float)sqrt(tr*tr + ti*ti);

The problem was that some spectra that should have been non-zero, e.g., tr and ti = 1.0e-20, were set to zero because of the multiplication inside the square root, even thought the square root is a double.

March 9, 2022

March 6, 2022

February 7, 2022



2021



December 11, 2021

September 20, 2021

September 19, 2021

WVFMT96    1.0  288.   74.  -91.   4.15     1.000 0.186E-12     1.000     1.000 0.154E-12  54.5 
where the columns are
 1 WVFMT96 if created by wvf96.f
 WVFMTD96 if created by wvfmtd96
 The next four columns describe the best (major) double couple derived from the moment tensor
 2  depth - depth of the Greens functions
 3 stk - strike of best double couple
 4 dip - dip of best double couple
 5  rake - rake of best double couple
 6  Mw  - moment magnitude of best double couple
 7 Reduction of Variance: 1.0-total_err/total_sum_sq - raw data 0 zero is best
 8 Std Error of fit: sqrt(total_err*tmpsig) 
 9 xy/sqrt(xx*yy) - cosine of angle between observed and predicted
10  Weighted reduction of variance: 1.0-total_wtd/total_sum_sq_wt 
11  Weighted Std error of fit: sqrt(total_wtd*tmpsig) 
12  Percent CLVD - 0 means pure double couple



The new output is
WVFMT961    1.0  288.   74.  -91.   4.15     1.000 0.186E-12     1.000     1.000 0.154E-12  54.5 -0.1010000E+23 -0.1060000E+23  0.3000003E+21  0.8799998E+22  0.2799999E+22 -0.2190001E+23
The program name in the first column changes and there are 6 additional columns


        1 WVFMT961  if
        created by the updated wvf96.f

          WVFMTD961  if
        created by the updated wvfmtd96

        Columns 2 - 12 are the 
        13  Mxx dyne-cm)
        14  Myy (dyne-cm)
        15  Mxy (dyne-cm)
        16  Mxz (dyne-cm)
        17  Myz (dyne-cm)
        18  Mzz (dyne-cm)

      

September 15, 2021

September 4, 2021

The reason for this option is that users may wish to use another program to plot the model parameters as a function of depth.

Examples of the plots created are given in SHWMOD96/index.html

August 31, 2021

June 2, 2021

February 1, 2021

January 31, 2021



2020



October 29, 2020

October 1, 2020

'eml_H.xy' 0 0.01 'CI' 0.05  '1 sig'
where eml_H.xy has entries   
   2.5000000       3.3192308       0.0000000      0.40308103   
   2.7000000       3.3750000       0.0000000      0.29140610   
   2.9000001       3.4560604       0.0000000      0.28769836   

The result is shown in the next figure:

genplt

This figure was created by concatenating the CALPLOT files of three invocations of genplt:
#####
#   plot the local earthquake Mw vs ML
#####cat > lcmdfil << EOF
'local.xy' 4 0.01 'CI' 0.05 'Local'
EOF
genplt -XLEN 5 -YLEN 5 -X0 2 -Y0 1 -XMIN 1.0 -XMAX 7.5 -YMIN 1.0 -YMAX 7.5 -TX "ML (H)" -TY "Mw" -L lcmdfi
mv GENPLT.PLT LSLU.PLT
#####
# plot the SLU Mw vs ML
#####
cat > ecmdfil << EOF
'eml_H.xy' 0 0.01 'CI' 0.05 '1 sig'
EOF
genplt -XLEN 5 -YLEN 5 -X0 2 -Y0 1 -XMIN 1.0 -XMAX 7.5 -YMIN 1.0 -YMAX 7.5 -TX "ML (H)" -TY "Mw" -E ecmdfil
mv GENPLT.PLT ESLU.PLT
#####
# plot the SMSIM bounds as a polynomial. The first line says use the file SMSIMML.xy, plot in red, and fill
#    The second line says use the file SMSIMML.xy, use the color black, and only plot the outline of the polygon
#    The result is a shaded red area outlined in black
#####
cat > pcmdfil << EOF
'SMSIMML.xy' 2 1
'SMSIMML.xy' 1 0
EOF
#####
#    note that since there was no -LPOS "TR" on the command line, no legend is plotted.
genplt -XLEN 5 -YLEN 5 -X0 2 -Y0 1 -XMIN 1.0 -XMAX 7.5 -YMIN 1.0 -YMAX 7.5 -TX "ML (H)" -TY "Mw" -P pcmdfil
mv GENPLT.PLT SMSIM.PLT
#####
#      concatenate placing red on bottom, then results with error bars and then the local
#####
cat SMSIM.PLT ESLU.PLT LSLU.PLT > ALL.PLT



September 16, 2020

May 20, 2020

February 24, 2020



2019



December 30, 2019

hprep96 -M model -d dfile 
hspec96
hpulse96 -V -p -l 1
hprep96 -M model -d dfile
rspec96
hpulse96 -V -p -l 1

February 19, 2019



2018



November 30, 2018


With the -XDOWN flag


Without the -XDOWN flag

June 5, 2018

 -LDOT ldot (default none ) dotted line with length ldot
 -LDOTDASH ldotdash (default none ) dash dotted line with length ldotdash
 -LDASH ldash (default none ) dashed line with length ldash
For the default axes lengths, a value of ldot, ldotdash or ldash = 0.05 is OK, e.g., -LDOT 0.05

May 1, 2018

April 15, 2018

March 1, 2018

February 22, 2018



2017



December 31, 2017

October 16, 2017

September 12, 2017

#ifdef DEBUG
fprintf(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, Argentina

July 7, 2017

May 27, 2017

May 1, 2017



2016



May 8, 2016

April 28, 2016

March 17, 2016

February 17, 2016

 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
#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.592
1969 01 01 (001) 18 27 23 408
#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.592
1969 01 01 18 27 23.408  1969001 Jan  1,1969 18:27:23.408
T(p) = ∑iNHiVi(1-p2Vi2)12∑_i^N \frac{H_i}{V_i (1 - p^2V_i^2)^\frac{1}{2}}    
X(p)=∑iHipVi(1-p2Vi2)12X(p) = ∑_i \frac{H_i p V_i}{ ( 1 - p^2 V_i^2) ^  \frac{1}{2} }    
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(p+Δp,h+Δh)=T(p,h)+(∂T/∂p)Δp+(∂T/∂h)Δh
X(p+Δp,h+Δh)=X(p,Δ)+(∂X/∂p)Δp+(∂X/∂h)ΔhX(p+Δp, h+Δh) = X(p,Δ) + (∂X/∂p) Δp + (∂X/∂h) Δ
Since we do not want the X to change, the change in p is related to the change in h, and thus
(∂T∂h)fixedX=(∂T∂h)fixedp-(∂T∂p)fixedh(∂X∂p)fixedh-1(∂X∂h)fixedp(\frac{∂T}{∂h})_{fixed X} = (\frac{∂T}{∂h})_{fixed_p} - (\frac{∂T}{∂p})_{fixed_h } ( \frac{∂X}{∂p})_{fixed_h}^{-1} (\frac{∂X}{∂h})_{fixed_p}

      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.

example of cut vel command

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

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

example of xlim command

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



2015



September 11, 2015

The help psppk give the following:  
SUMMARY:
     Interactive spectra pick
PlotSPPK
     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 marking
DESCRIPTION:
     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.WSTK
An example of the use of the Z-Z option of pkppk is in the following image
example of t*
        determination

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 40
r *Z
fileid name
background on color 1
color rainbow
fileid name
sort up dist
#p relative
rtr
taper w 0.25
#####
#    differentiate ground velocity to get
#    acceleration. For large earthquake
#    acceleration spectrum is flat, any trend
#    is due to tstar
dif
fft
psppk xlin fmin 0 fmax 4 overlay on  smooth on
q
EOF
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 information 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






2014



November 11, 2014

June 10,2014

May 24, 2014

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.

May 6, 2014

April 28, 2014

                at line516

March 28, 2014

March 13, 2014

February 7, 2014

February 6, 2014



2013



November 8, 2013

August 22, 2013

August 8, 2013

July 26, 2013

        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

May 22, 2013

                    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.


2012



December 15, 2012

December 10, 2012

October 22, 2012

Original code

Corrected code

September 20, 2012

September 9, 2012

August 6, 2012

July 26, 2012

April 9, 2012

March 26, 2012

February 8, 2012

January 14, 2012



2011



November 23, 2011

November 1, 2011

October 20, 2011

October 13, 2011

September 23, 2011



2010



October 24, 2010

September 20, 2010

September 17, 2010

August 30, 2010

August 20, 2010

July 15, 2010

June 17, 2010

June 14, 2010

June 6, 2010

May 20, 2010

April 7, 2010

March 18, 2010

graysc.png

grayscinv.png

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

February 16, 2010

January 11, 2010

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

                                } 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;



2009



December 4, 2009

October 8, 2009

July 22, 2009

 June 19, 2009

June 11, 2009

May 23, 2009

March 30, 2009

March 19 ,2009

March 18, 2009

March 3, 2009

February 21, 2009

February 18, 2009

January 5, 2009



2008



December 22, 2008

November 12. 2008

September 3, 2008

August 4, 2008

June 13, 2008

May 25-29, 2008

May 7, 2008

April 9, 2008

March 31, 2008

March 14, 2008

February 8, 2008

January 28, 2008

January 22, 2008

January 17, 2008

January 10, 2008

January 3, 2008

January 2, 2008



2007



November 17, 2007

November 9, 2007

October 21. 2007

October 19. 2007

October 04, 2007

August 30, 2007

August 19, 2007

August 16, 2007

August 14, 2007

August 1, 2007

July 16, 2007

July 1, 2007

April 23, 2007

March 31, 2007

March 19, 2007

March 15, 2007

February 22. 2007

February 20. 2007

February 8, 2007

January 24, 2007

January 13, 2007

January 3, 2007

January 3, 2007

January 3, 2007



2006



December 28, 2006

December 15, 2006

December 12, 2006

November 8, 2006

October 18, 2006

October 14, 2006

October 1, 2006

September 16, 2006

August 20, 2006

August 8, 2006

August 8, 2006

August 4, 2006

June 16, 2006

June 10, 2006

May 27, 2006

April 26, 2006

March 19, 2006

February 27, 2006

February 27, 2006

January 3, 2006



2005



December 14, 2005

September 25, 2005

August 10, 2005

August 3, 2005

June 30, 2005

                   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

June 16, 2005

March 23, 2005

February 11, 2005

February 5, 2005

January 31, 2005