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

click title to see Google Drive directory JParkCodes

c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. 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, and bootstrap-resamples the data to estimate the
c uncertainty of the harmonic terms. The posted version only
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 (see previous posts on codes that apply
c a frequency-domain moveout correction for the variation of Ps delay time with epicentral distance.)
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:
c the file is "in_recfunk" and has lines of the form:
c
c 1997.070.19.33.bh? <-- br="" code="" r="" replaces="" t="" with="" z="">c 57 52 <-- analysis="" br="" duration="" of="" sec="" start="" time="" window="">c 1997.076.08.15.bh?
c 62 62
c ...
c ...
c ...
c stop <-- 799="" br="" code="" data="" events="" finished="" is="" max="" tells="" that="">c
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