recfunk_estack.f -- an event-stacking algorithm for RF with long delay times

recfunk_estack.f -- an event-stacking algorithm for RF with long delay times

Click the title to access the Fortran code recfunk_estack.f in

recfunk_estack is a version of the multiple taper correlation receiver-function estimator that allows the time-domain RF to extend to times longer than (roughly) T/10, where T is the duration of the analysis window. Recfunk_estack does NOT apply a moveout or migration correction, so the ability to resolve Ps from deeper interfaces is limited. Other codes to be posted soon will have these properties.

The major shortcoming of the MTC RF estimator of Park & Levin (2000) is its difficulty retrieving Ps conversions at long delay times. Receiver function studies of upper mantle structure at the base of the lithosphere and the transition zone involve Ps converted phases with delay times between 20 and 70 seconds (e.g., Bostock, 1998; Sheehan et al 2000), but the frequency averaging properties of the Slepian data tapers with time-bandwidth product p\ge 2 tends to limit typical RFs to delay times \tau<20 br="" s.="">
Frequency averaging affords the MTC algorithm desirable statistical properties, especially for single-event RF estimates. We can achieve a narrower frequency average for the retrieval of Ps from interfaces deeper than 100 km if we sum single spectral estimates over K seismic events, rather than K eigenspectra from one event, in the RF cross-correlation formulas (2) and(3). Using only one spectrum estimate per component for each P wave, we estimate the spectra with a Slepian taper of time-bandwidth p=1, that is, the "one-\pi-prolate" taper. This taper averages local frequency information over a central bandwidth similar to the simple boxcar taper, but is optimized for spectral leakage resistance.

c  program recfunk_estack

c  5/11/01 JJP -->  REV2 5/16/01

c  updated 08/17/07 JJP, updated for GFORTRAN compiler 07/28/15

c  xf77 -o /park/backup/bin/recfunk_estack recfunk_estack.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a

c  xf77 -o /Users/jjpark/bin/recfunk_estack recfunk_estack.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a

c

c  will compute RFs using single eigenspectrum estimate per record

c  instead of cross-correlating eigenspectra in one event

c  --> cross-correlate across events

c  steps:

c  LOOP1

c  read the filenames to analyse, compute 1-pi prolate spectra 

c and coherence estimates

c  save baz/epicenter info

c  LOOP2

c  choose criteria for estack

c  for each group of events: loop over freq:

c assemble the spectral vectors, cross-correlate

c   with coherence weighting

c  compute the H(f), display

c

c

c  BASED ON

c  program recfunk_noplot

c  5/25/99 JJP

c  multitaper program to invert for receiver functions

c  SINGLE STATION WITH SPECTRAL WEIGHTING

c  

c  "Multiple Taper Correlation" (MTC) Receiver functions are described in 

c

c  Park, J., and V. Levin, Receiver functions from multiple-taper spectral 

c  correlation estimates, Bull. Seismo. Soc Amer., v90, pp1507-1520,  2000.

c  

c  spectral estimates of pre-event noise 

c  weight the inversion in the frequency domain, thereby avoiding the 

c  usual problem with microseism noise in 0.1-0.5 Hz

c

c  the main advantage of MTC is not the noise damping, however -- the

c  cross-correlation of the horizontal and vertical components is used

c  in place of spectral division when estimating the RF.

c  this approach appears to be more stable and gives RF-uncertainties

c  in the frequency domain that are useful for stacking.

c

c  this version of the RF code reads a file of data filenames 

c  you have two choices: either read the time intervals in the filename file

c  or read them in the sac header

c  the data must be binary SAC format

c  horizontals must be rotated to radial and transverse

c

c  for start times in the file:

c  the file is "in_recfunk" and has lines of the form:

c

c  1997.070.19.33.bh?  <-- code="" div="" r="" replaces="" t="" with="" z="">

c  57 52       <-- analysis="" div="" duration="" of="" sec="" start="" time="" window="">

c  1997.076.08.15.bh?

c  62 62 

c  ...

c  ...

c  ...

c  stop                <-- 799="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">

c

c  for start times in the SAC header

c  reads seismic record start times from the sac header

c  will search the SAC header for specific markers of P phases

c  T1 - P, Pdiff    ahead(12)

c  T2 - PKP,PKIKP   ahead(13)

c  T3 - PP          ahead(14)

c  T1=T2=T3=0 ==> use original A-marker ahead(9)

c

c   code does NOT combine data with different sample rates

c  data files limited to 99K pnts. To increase, see common block /datastuff/

c

c   many intermediate quantities are plotted with PLOTIT as the code proceeds.

c  other intermediate quantities can be plotted by uncommenting calls to PLOTIT

c

c   the program computes RFs in the Freq domain for multiple records

c   then it computes composite RFs that are dependent on either 

c   Back-azimuth or epicentral distance.  The standard bin width is 10 degrees

c   with stacked RFs computed every 5 degrees.  You can change this in the code.

c  RF are not migrated in this process, so the user is prompted for parameters 

c  that narrow the focus of the stacks 

c  e.g. a BAZ-sector over which to compute an epicentral RF profile.

c

c  the PLOTIT graphs of RF sweeps can be improved upon, natch

c  the code writes the BAZ- and EPICEN-dependent RFs to files

c  in a format easily digested by GMT (traces are separated by ">" lines)

c

c  filenames: out[rt]_baz.grid out[rt]_epi.grid

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c