Wednesday, July 29, 2015

rfmig_cboot.f -- code to compute harmonic expansion of RFs with back-azimuth

click title to see Google Drive directory JParkCodes
c  program rfmig_cboot
c  10/12/04 JJP -- adapted from rfmigrate   -- UPDATED 12/01/07 JJP and 07/29/15
c  xf77 -o /Users/jjpark/bin/rfmig_cboot rfmig_cboot.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c  code computes frequency-domain stacks of receiver functions that 
c  follow a harmonic expansion in baz for both radial and transverse RFs 
c  there are 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, and bootstrap-resamples 
c    the data to estimate the uncertainty of the harmonic terms.
c   output files are out[rt]_cexp.grid  -- harmonic-expansions of the RFs, in time domain
c     out[rt]1_cexp.grid  -- harmonic-expansion RFs plus bootstrap uncertainty
c     out[rt]2_cexp.grid  -- harmonic-expansion RFs minus bootstrap uncertainty
c     out[rt]_bbaz.grid  -- harmonic-expansion RFs computed for ordered baz values
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 the output of this program is a least-square regression of RFs in the freq
c domain, with complex-valued coefficients for the 
c constant, cos(baz)R+sin(baz)T, cos(baz)R-sin(baz)T,
c cos(2*baz)R+sin(2*baz)T, cos(2*baz)R-sin(2*baz)T variations
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  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  for start times in the file of filenames
c  the file defaults to "in_recfunk" and has lines of the form:
c  <-- code="" div="" r="" replaces="" t="" with="" z="">
c  57 52       <-- analysis="" div="" duration="" of="" sec="" start="" time="" window="">
c  62 62 
c  ...
c  ...
c  ...
c  stop                <-- 799="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">
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  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  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   many intermediate quantities are plotted with PLOTIT as the code proceeds.
c  other intermediate quantities can be plotted by uncommenting calls to PLOTIT

No comments: