## Friday, August 1, 2008

### Rayleigh Wave Dispersion Code chenray.f

click title to see Google Drive directory JParkCodes

This a a computer code that performs the 1D layered-media reflectivity method to compute Rayleigh-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them. The code only plots dispersion curves for phase velocity CVEL, but does not compute GVEL, the group velocity. The values of CVEL are *not* written out for later use. You can modify the code yourself to write out cvel values.

c chenray - to calculate propagating Rayleigh modes of layered structure
c revised to calculate dispersion curves 7/7/95, debugged 7/27/08
c 9/2/94 - cannibalized from aniprop to test its structure
c xf77 -o /home/jjpark/bin/chenray chenray.f /home/jjpark/Plotxy/plotlib.a /home/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/chenray chenray.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c
c will calculate Rayleigh wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c ee, e1 and e2 are 4x4 matrices that get partitioned into 2x2 submatricies to
c separate out the reflection and transmission of waves at an interface.
c
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, but does not write phase velocity (cvel) to a file
c or compute the group velocity (gvel). The code has not worked in 1994, and I chose to
c push on to the general anisotropic surface wave code aniprop.f, rather than to debug this
c Rayleigh-wave code. This code was revived and debugged in 2008(!) by JPark
c

### Love-Wave Dispersion Code: chenlove.f

click title to see Google Drive directory JParkCodes

chenlove.f

This a a computer code that performs the 1D layered-media reflectivity method to compute Love-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them.

c will calculate love wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c the variables:
c c(j,1) are C_d; c(j,2) are the C_u from Chen (1993) for the jth layer
c c(nlayer+1,*) are coefs in the halfspace.
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, and also writes phase velocity (cvel)
c group velocity (gvel) and three integration scalars that are used in the variational
c computation of group velocity from phase velocity
c do the sum
c ...
c ci1=ci1+rho(il)*(dci1+dci2)
c ci2=ci2+xmu(il)*(dci1+dci2)
c ci4=ci4+xmu(il)*(dci1-dci2)*xnu(il)**2
c end do
c gvel=ci2/(cvel*ci1)
c
c

## Tuesday, July 1, 2008

### agu.f Viewing the output of the sacw_corr_many code

click title to see Google Drive directory JParkCodes

The code agu.f (guess where I wrote it!) reads the cvel_YYYY.DDD.HH files and generates some plots of the phase-velocity as fct of period using a collection of algorithms: weighted sum, trimmed weighted sum, median. Version posted can process up to 100 cvel files. You can expand the arrays for more.

input:

cvel_2006.175.02 --> input file (one per event , just "ls cvel_20* >! c_agu"
cvel_2006.225.10
cvel_2006.225.19
cvel_2006.240.22
stop --------------------> end of input files
3 --------------------> make a plot
4 --------------------> make a plot and write to PS file
cvel_eastwest.ps ---------> PS plotfile

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

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

### Code for surface-wave phase-velocity estimation sacw_corr_many.f

click title to see Google Drive directory JParkCodes

c program sacw_corr_many
c JJP 04/13/07 updated 7/24/07 updated 12/11/07, again 6/25/08
c compute phase delays as function of frequency using the multiwavelet transform
c input two fil,es, typically BHZ components at different stations
c cross-correlate the wavelet transforms
c output is a colormap of time and frequency dependent time delay
c plotted with the wavelet transform of the two time series
c
c f77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c f77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr.f /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr_many.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c to read SAC-format data, filename.bh? <-- br="" e="" for="" n="" routine="" subs="" z="">c
c generates wavelets with eigenvector decomposition + interpolation
c calculate padded FFT of data for convolution with wavelet FFTs
c general use - calculates set of wavelets for a specified number of points
c interpolates and fft other lengths on the fly
c max=12 real-wavelets max = 6 complex-valued wavelets
c
c code is based on the algorithms described in
c
c Lilly, J., and J. Park, Multiwavelet spectral and polarization analysis
c of seismic records, Geophysical J. International, v122, 1001-1021, 1995.
c
c Park, J. and M. E. Mann, Interannual temperature events and shifts in
c global temperature: A multiwavelet correlation approach, Earth Interactions,
c v. 4, online paper 4-001, 2000.
c
c Bear, L. K., and G. L. Pavlis, Estimation of slowness vectors and
c their uncertainties using multi-wavelet seismic array processing,
c Bull. Seismol. Soc. Am., v87, 755--769, 1997.
c
c this code generates 2-D arrays of spectra, polarization and statistical stuff
c it writes GMT scripts to display these arrays in dazzling color
c
c for an input filename sacfile.bh? the GMT scripts are written to
c GMD_sacfile.bh
c
c
c the ASCII files to plot are named
c cpt_nilcut.sh spw_nilcut.sh1 spw_nilcut.sh2
c dat_nilcut.sh
c cpd_nilcut.sh del_nilcut.sh
c
c dat_* is data to plot
c cpt_* and cohere.cpt are colorfiles for GMT
c del_* dep_* and dem_* are the time delay arrays
c spw_* wavelet spectra
c
c you can fill a lot of disk space if you are not tidy with old files!
c
c user must make scripts executable ("chmod +x GM*_sacfile.bh"),
c execute the scripts and display the postscript files

## Tuesday, May 27, 2008

### changes in the inverse FFT in recfunk codes

click title to see Google Drive directory JParkCodes

After some comparison runs, Vadim Levin convinced me that the "improvement" I made to the inverse FFT in the recfunk*, rfmig* codes was a bad innovation.  The original RF-estimation codes applied a cos^2 factor to the frequency-domain RFs H_R(f) and H_T(f) with weighting unity at the f=0 point, decreasing to zero at f_max while following a cos-squared dependence.  I changed all codes to have a sharper bandpass (I was unhappy that the half-amplitude point occurred at f_max/2!) with unit factor from f=0 to f=0.5*f_max, then a cos-squared rolloff within (0.5*f_max , f_max).  This sharper bandpass led to small negative sidelobes at the flanks of sharp Ps pulses.  Although one can live with this, the temptation to interpret a wiggle could be strong, so it is safer to use the original weighting of H(f) before inverse FFT to obtain the time-domain RF.  The sidelobes do not disappear completely, but are far smaller.

recfunk_pick.f
recfunk_estack.f
rfmigrate_estack.f
rfmigrate.f
recfunk_svd.f
recfunk_mwm.f
recfunk_mwmj.f
recfunk_mwms.f
rfmig_boot.f
rfmig_mboot.f
rfmig_cboot.f
rfmig_mcboot.f

### Update of anirec_synth.f and anirec_synth_circle.f

click title to see Google Drive directory JParkCodes

A small amount of random noise was added to synthetics to prevent zero-divide errors (in cases where the recfunk codes tries to compute an RF for an interval of x=0 data). In the update, this bit of noise was made smaller, and explicitly zero-mean. There were some low-frequency artifacts cropping up in synthetic RFs. The updated versions seem to avoid this problem:

anirec_synth.f

anirec_synth_circle.f

## Monday, March 24, 2008

### The erstwhile version of aniprop.f has had a bugfix

click title to see Google Drive directory JParkCodes

The erstwhile version of aniprop.f that was (is?) available on my website has been superceded by an all-in-one version that computes dispersion curves and synthetic seismograms all in one program. Adapted and debugged on an Intel Mac in gnu fortran.