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 For each year, the changes are (listed in in order from newest to oldest)


May 1, 2017


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
(Th)fixedX=(Th)fixedp-(Tp)fixedh(Xp)fixedh-1(Xh)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




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

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



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.


December 15, 2012

December 10, 2012

October 22, 2012

September 20, 2012

September 9, 2012

August 6, 2012

July 26, 2012

April 9, 2012

March 26, 2012

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



    November 23, 2011

    November 1, 2011

    October 20, 2011

    October 13, 2011

    September 23, 2011



    October 24, 2010

    October 12, 2010

    October 3, 2010

    September 30, 2010

    September 20, 2010

    September 17, 2010

    August 30, 2010

    	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

    July 15, 2010

    June 17, 2010

    June 14, 2010

    June 6, 2010

    May 20, 2010

    April 7, 2010

    March 18, 2010

    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;



    December 4, 2009

    October 8, 2009

           Summary of changes made since July 22, 2009

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

    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

  • 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

    January 5, 2009

    January 2, 2008



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



    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



    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