rfmig_boot.f -- code to estimate harmonic expansion of RFs in back azimuth

rfmig_boot.f -- code to estimate harmonic expansion of RFs in back azimuth

click title to see Google Drive directory JParkCodes

c  program rfmig_boot

c  10/12/04 JJP -- adapted from rfmigrate  -- UPDATED 12/01/07 JJP, and 07/27/15

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

c

c  code computes frequency-domain stacks of receiver functions that follow 

c  a harmonic expansion in baz for both radial and transverse RFs there are 

c  constant terms and sin/cos terms for 2- and 4-lobed amplitude dependence.  

c  The constant term should be zero for the transverse RF.  

c  The 2-lobed terms govern dipping interface effects and tilted symmetry-axis ansotropy. 

c  The 4-lobed term is anisotropy with a horizontal axis

c  The code regresses for the harmonic expansion, 

c  and bootstrap-resamples the data to estimate the uncertainty of the harmonic terms.

c  

c   output files are out[rt]_bexp.grid  -- harmonic-expansions of the RFs, in time domain

c     out[rt]1_bexp.grid  -- harmonic-expansion RFs plus bootstrap uncertainty

c     out[rt]2_bexp.grid  -- harmonic-expansion RFs minus bootstrap uncertainty

c     out[rt]_bbaz.grid  -- harmonic-expansion RFs computed for ordered baz values

c

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

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 of filenames

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 computes a 

c  the code writes the BAZ-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  these files are overwritten everytime you run the program

c  so rename them to save them

c