Tuesday, July 28, 2015

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

No comments:

 
Link