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

recfunk_pick.f -- version of multiple-taper RF estimation

recfunk_pick.f -- version of multiple-taper RF estimation

click title to obtain code from Google Drive directory JParkCodes

This code is a simple variant of the multiple-taper RF estimator recfunk.f, which has been public for years on my Yale website. This version reads the start time of the P-wave analysis window from the "AMARKER" entry of the SAC header. It also rotates the vertical and radial seismogram components into LQT coordinates using a formula based on epicentral distance. This formula for LQ rotation angle is very crude, but simple measures of the actual rotation angle scatter quite a bit, and every station has unique behavior. This variability suggests that careful derivation of a precise LQ rotation angle is a fool's game.  In this code there is a simple formula used, but most of my other RF codes use the slowness of the incoming P wave to estimate the rotation angle.  This choice works very well for synthetic seismograms, but real data can be diverse in behavior.

c  program recfunk_pick -- UPDATED
c  updated for GFORTRAN compiler 07/28/15
c
c  original code written for the paper
c  Park, J., and V. Levin, Receiver functions from multiple-taper
c  spectral correlation estimates, Bull. Seismo. Soc Amer., v90,
c  pp1507-1520, 2000.
c
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  the code reads filenames from a file that defaults to "in_recpick"
c  the user inputs the length of the data window to be processed after the
c  time marker in the SAC header.
c
c  08/16/07 JJP
c  multitaper program to compute receiver functions
c  SINGLE STATION WITH SPECTRAL WEIGHTING
c
c  uses spectral estimates of pre-event noise
c  to weight the inversion in the frequency domain, thereby avoiding the
c  usual problem with microseism noise in 0.1-0.5 Hz
c  IN PRACTICE, however, the regularization of what might be a spectral division is achieved
c  via multiple taper cross-correlation
c
c  Program has extra code to perform chi-squared tests, discussed in Park and Levin paper
c
c  xf77 -o /park/backup/bin/recfunk_pick recfunk_pick.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c  xf77 -o /Users/jjpark/bin/recfunk_pick recfunk_pick.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c

recfunk.f -- original recfunk version without moveout corrections

recfunk.f -- original recfunk version without moveout corrections

Click the title of the post to access the code recfunk.f on Google Drive directory JParkCodes.

c  program recfunk

c  5/25/99 JJP updated for GFORTRAN compiler 07/27/15

c

c  original code written for the paper 

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

c  spectral correlation estimates, Bull. Seismo. Soc Amer., v90, 

c  pp1507-1520, 2000.

c

c  this version of the RF code reads datafile names from standard input

c  then the time intervals in the filename file

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

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                <-- 399="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">

c

c  has kluge to cheat the pre-event noise for synthetic records  3/12/00 JJP

c  check to see if the kluge is commented out

c

c  multitaper program to invert for receiver functions

c  SINGLE STATION WITH SPECTRAL WEIGHTING

c  

c  uses spectral estimates of pre-event noise 

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

c  usual problem with microseism noise in 0.1-0.5 Hz

c

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

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

rfmigrate.f -- recfunk version with migration

rfmigrate.f -- recfunk version with migration

Click the title of the post to access the code rfmigrate.f on Google Drive directory JParkCodes.

c  program rfmigrate

c  9/20/00 JJP -- adapted from recfunk  --- UPDATED 08/19/07 JJP

c   updated again 07/22/15 JJP

c  migrates MTC receiver functions in the frequency domain.

c  requires a stacking model in the anirec format

c  such a model may have anisotropy parameters, 

c  but migration code only uses the isotropic velocities. 

c

c  has kluge to cheat the pre-event noise for synthetic records  3/12/00 JJP

c  check to see if the kluge is commented out, because anirec_synth was updated 

c  in 2015 to have a proper-length pre-event time window.

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 defaults to "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

c  for start times in the SAC header, the default file is "in_recpick"

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  legacy data has timing in the

c   A-marker (ahead(9)) which doesnt indicate which phase was marked, 

c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 div="" nbsp="">

c  PKP/PKIKP for DELTA>120

c

c  code does NOT combine data with different sample rates

c  it is possible to spline interpolate the data files for a uniform sample rate

c  see the code in rfmig_mcboot.f, which should be spliceable into this code.

c  data files are 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 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 oum[rt]_epi.grid oum[rt]_baz.grid

c

c  out --> no moveout correction;  oum --> moveout correction

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c

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

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

c

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

click title to see Google Drive directory JParkCodes

c  anirec_synth_circle - anirec to compute synthetics for a specified model and
c  a specified station location from one seismic record read from in_recfunk
c  the code will take the station location and generate a ring of 120 events
c  that are 90 degrees away, one for every 3 degrees back-azimuth.
c  To do this we generate a great circle path
c  with the station at the pole of the great circle
c
c  To minimize aliasing in the time domain, we create synthetics that last 200 seconds
c  The inverse FFT of the frequency-domain reflectivity response is time shifted to
c  place the P wave just beyond the 100-sec mark, using a wraparound of the inverse FFT.
c  The code can use several pulse shapes, according to the flag ipulse
c   ipulse=1 -- one-sided pulse (squared sine wave) with user-chosen period
c   ipulse=2 -- two-sided sine pulse with user-chosen period
c   ipulse=3 -- two-sided dispersive pulse, 160-points
c   ipulse=4 -- two-sided 10-second dispersive pulse with random-walk phase
c   ipulse=5 -- random-walk pulse (integration of RAND) -- emergent 10 sec pulse with tapered end
c
c   This code is hardwired for incident P waves (ity=1), but the reflectivity calculation has info
c   to construct incident SV (ity=2) and SH (ity=3) synthetics as well.  The exact timing of these
c   S-coda are different from the P-coda, so you would need to check the wraparound timing to get it right.
c
c  code accepts as input SACfiles with start time info in the header
c  legacy data has timing in the
c   A-marker (ahead(9)) which doesnt indicate which phase was marked,
c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 nbsp="" p="">c  PKP/PKIKP for DELTA>120
c  the value of A-marker is not important.
c  In the output files of synthetics, the A-marker and/or the T1,T2,T3-marker
c  in the header are set to 95.0 seconds to allow for some pre-event noise
c  that is generated by intrinsic function RAND
c  for A.ne.0 by recfunk codes that use the header information for timing.
c  code divides between T1 (P) or T2 (PKP) or T3 (PP) markers based on epicentral distance,
c  subroutine call_ttimes reads files that give slowness values for 0c  Preactically, these wont be used for all DELTA, but consistency!
c  the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c  the code writes in_recpick_synth, a list of the filenames for recfunk_pick
c  and other recfunk codes that expect the start time in the SAC header
c
c  written originally 12/3/03 JJP
c  3/3/04 -- Added Kluge -- I have added code that forces the phase velocity to
c  fluctuate with BZ to mimic the effect of a dipping layer,
c  changing its incidence angle with a cos(\phi) dependence.
c  This kluge mimics the P-SV conversion variation caused by a dipping interface
c  but not the P-SH conversion behavior.
c  This kluge was used for the last check on RF interpretation in
c  Park et al 2004 in JGR, the paper on the Cascadia subduction zone.
c  Look for the word "kluge" in the code and make certain that the code
c  between begin-kluge and end-kluge comment-lines is commented
c
c  xf77 -o /Users/jjpark/bin/anirec_synth_circle anirec_synth_circle.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  anirec - program to calculate receiver-function response of
c  a stack of anisotropic layers to a plane wave incident from below
c  cannibalized from aniprop.f  11/18/95
c  for hexagonally symmetric media
c  reads fast axis orientation, constants A,B,C,D,E from file animodel
c  calculate quadratic eigenvalue problem based on the Christoffel matrix
c  see appendix of P. Shearer's thesis
c
c  read model, phase velocity of incident wave, P, SV, or SH
c
c  calc the eigenvector decomps for the layers
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
c  find distortion of reference wavelet
c

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

click title to see Google Drive directory JParkCodes

c  anirec_synth - anirec to compute synthetics for a specified model and
c  a specified collection of seismic records
c  (read from in_recpick-format file of filenames)
c the goal is to create P-coda synthetics that match the event distribution of your dataset.
c
c  To minimize aliasing in the time domain, we create synthetics that last 200 seconds
c  The inverse FFT of the frequency-domain reflectivity response is time shifted to
c  place the P wave just beyond the 100-sec mark, using a wraparound of the inverse FFT.
c  The code can use several pulse shapes, according to the flag ipulse
c   ipulse=1 -- one-sided pulse (squared sine wave) with user-chosen period
c   ipulse=2 -- two-sided sine pulse with user-chosen period
c   ipulse=3 -- two-sided dispersive pulse, 160-points
c   ipulse=4 -- two-sided 10-second dispersive pulse with random-walk phase
c   ipulse=5 -- random-walk pulse (integration of RAND) -- emergent 10 sec pulse with tapered end
c
c   This code is hardwired for incident P waves (ity=1), but the reflectivity calculation has info
c   to construct incident SV (ity=2) and SH (ity=3) synthetics as well.  The exact timing of these
c   S-coda are different from the P-coda, so you would need to check the wraparound timing to get it right.
c
c  code accepts as input SACfiles with start time info in the header
c  legacy data has timing in the
c   A-marker (ahead(9)) which doesnt indicate which phase was marked,
c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 font="" nbsp="">
c  PKP/PKIKP for DELTA>120
c  the value of A-marker is not important.
c  In the output files of synthetics, the A-marker and/or the T1,T2,T3-marker
c  in the header are set to 95.0 seconds to allow for some pre-event noise
c  that is generated by intrinsic function RAND
c  for A.ne.0 by recfunk codes that use the header information for timing.
c  code divides between T1 (P) or T2 (PKP) or T3 (PP) markers based on epicentral distance,
c  subroutine call_ttimes reads files that give slowness values for 0
c  Preactically, these wont be used for all DELTA, but consistency!
c  the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c  the code writes in_recpick_synth, a list of the filenames for recfunk_pick
c  and other recfunk codes that expect the start time in the SAC header
c
c  modified 11/30/07 JJP
c  5/12/01 JJP, 7/23/2015
c  xf77 -o /Users/jjpark/bin/anirec_synth anirec_synth.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  anirec - program to calculate receiver-function response of
c  a stack of anisotropic layers to a plane wave incident from below
c  cannibalized from aniprop.f  11/18/95
c  for hexagonally symmetric media
c  reads fast axis orientation, constants A,B,C,D,E from file animodel
c  calculate quadratic eigenvalue problem based on the Christoffel matrix
c  see appendix of P. Shearer's thesis
cawk '{print "s_" $1}' in_recpick >! in_recpick_synth

c  read model, phase velocity of incident wave, P, SV, or SH
c
c  calc the eigenvector decomps for the layers
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
c  find distortion of reference wavelet
c

Updated Fortran Programs for Receiver-Function Processing

Sportsfans,

After upgrading my Macs to MacOS X 10.10.4 Yosemite, I loaded the gcc5 compiler from Fink.org and started modifying my codes to compile properly with the gfortran compiler that comes with Rev5 of gcc.  Although the codes are still written mostly in the oldest flavors of Fortran, they will indeed compile and run in the new computing environment.  I will post descriptions of the updated codes over the next days, and load them into my Google Drive directory.

One update is that the RF codes that implement a moveout correction while stacking have been altered to use P-wave slowness files for the P, PKP-PKIKP and PP cases.  These files are newly computed for the iaspei91 model from the TauP codes that superseded the venerable ttimes program in the IRIS software library.  You will see that the subroutine read_ttimes opens three files in a subdirectory of an Applications folder where TauP sits.  I will upload my own versions of these slowness files to Google Drive. These files average over triplications and allow Pdiff to reach DELTA=180 and follow the PKiKP phase to its zero-slowness value at DELTA=0.

look for

taup_P.dat
taup_PP.dat
taup_PKP.dat