sacw_corr_many.f -- Using the wavelet-based surface wave code

click title to see Google Drive directory JParkCodes

sacw_corr_many.f does a 2-station estimation of phase velocity.

sacw --> sac-format data, analyzed with wavelets

corr --> cross-correlation is used to estimated phase shits and phase velocity

many --> code is set up to run a series of data-record pairs.

The posted code is designed for 10 sps data, and computes Slepian wavelets in 16 banks logarithmically spaced between 8 and 56 second center-period. Three complex-valued Slepian wavelets are used, in order to estimate phased cross-correlations with 2 and 4 degrees of freedom. The phase of the cross-correlation is scaled with the center-period of the wavelet to compute a phase-delay. This phase delay is computed at evenly spaced points in the seismogram. The station locations and the back-azimuth of the earthquake are used to estimate a phase velocity from the phase delay, but only at the time point where the maximum wavelet spectrum is found. The wavelet spectrum is computed only from the first Slepian wavelet, because experience demonstrates that this affords a more stable determination of the max-amplitude arrival time.

The input is at the terminal, so the command

sacw_corr_many < c_files_process

will be standard analysis method. This command file has the form (for stations PIIR and VOLR)

2003.325.04.12.54.PIIR.DHZ
2003.325.04.12.56.VOLR.DHZ
3 ---> plot to screen, could be "0" to skip plotting
4 ---> plot to screen and write output to PS file (also could skip plotting with "0")
2003.325.ps ---> name of the PS file for phase-velocity graph
2004.024.13.05.37.PIIR.DHZ
2004.024.13.05.40.VOLR.DHZ
3
4
2004.024.ps
stop ----> end of data input


The data files must have a common start time (despite the differing numbers in the filenames in this example) In practice, one can use SAC commands to achieve this easily. It is convenient to drop the minutes and seconds entries in the filenames, saving the year/julianday/hour filename format, because these parameters will match for different stations, and the minutes/seconds wont.. For instance, if you have a file "infiles" of the filename pairs to analyze, then

ls 20*VOLR*DHZ | sed 's/\./ /g' | awk '{print "read " $1 "." $2 "." $3 "*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off"}' >! c_synch

will generate the SAC macro to synchronize and cut to a common start point:

read 2003.325.04*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.009.19*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.024.13*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.038.21*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off

and so on . . . .

The code writes a small ASCII file of phase-vel information, cvel_YYYY.DDD.HH, which looks like this

cat cvel_2004.339.10
56.00 1.547 1.858 13.76 236.6 341.4 7.340
49.18 4.596 12.74 13.76 236.6 342.0 2.470
43.20 2.893 1.626 13.76 236.6 319.2 3.925
37.93 4.331 1.237 13.76 236.6 321.0 2.621
33.32 3.411 0.3311 13.76 236.6 322.8 3.328
29.27 3.258 0.9897E-01 13.76 236.6 334.8 3.484
25.70 2.934 0.2715 13.76 236.6 334.8 3.869
22.58 2.948 0.1586 13.76 236.6 347.4 3.850
19.83 2.939 0.2432 13.76 236.6 358.2 3.863
17.42 2.759 1.053 13.76 236.6 372.0 4.114
15.30 3.739 2.711 13.76 236.6 370.2 3.036
13.43 5.761 2.135 13.76 236.6 417.6 1.971
11.80 4.417 0.4976 13.76 236.6 456.0 2.570
10.37 1.556 0.1375 13.76 236.6 559.2 7.295
9.100 2.039 0.4498 13.76 236.6 725.4 5.568
8.000 -27.40 99.12 13.76 236.6 995.4 -0.4143

period cvel d(cvel) gcarc baz timeOcorrelation delay-time