rfmigrate_estack.f -- code to migrate the RFs in the frequency domain

rfmigrate_estack.f -- code to migrate the RFs in the frequency domain

Click the link to access the code rfmigrate_estack.f  on Google Drive

Long-time window event-stacking of RFs. This code takes three swipes at the data. First it does a normal event-stack without migration in back-azimuth bins (output files out[rt]_baz.grid). Then a normal event-stack without migration in epicentral bins (output files out[rt]_epi.grid). Lastly it does an event-stack with migration in back-azimuth bins (output files oum[rt]_baz.grid).

c  program rfmigrate_estack

c  5/10/02

c  adapted from recfunk_estack

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

c  updated 08/17/07 JJP

c  updated 07/23/15 JJP

c

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

c  xf77 -o /Users/jjpark/bin/rfmigrate_estack rfmigrate_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

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 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 oum[rt]_baz.grid

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c