recfunk_mwms.f -- moving-window migration that splices together a long RF from a series of target depths

recfunk_mwms.f -- moving-window migration that splices together a long RF from a series of target depths

This code can be accessed by clicking the title link

recfunk_mwms.f generates a succession of RF estimates for a succession of migration target depths. It splices together the results into a long RF. I have used this code to reconstruct BAZ-sweeps of RFs from 0 to 75-s delay, encompassing the crust to the 660-km discontinuity.

This code will warm your laptop. For a splice from 0 to 680 depth, targeted every 40 km with 524 events, my Macbook Pro required 15 minutes of run time.

INPUT 
'velocity model for stacking (blank = stack_model)'
'depth range and inc for the RFs (km) (dep1,dep2,ddep)'
'enter fmax (Hz)'
'time intervals in file (0) or in SAC-header (1)'
'file of seismic records to read? (space-ret: in_recpick)'
'duration of data windows for analysis'
'minimum number of events in a binsum (eg 1 or 2)'
'rotate to LQT coordinates to isolate P and SV? (0=no)'

--->> then the first pass through the seismic data -- typically for target depth = 0 km.

'change baz-interval or baz-spacing? (1=yes)'
'enter back-azimuth range: bazmax, bazmin, bazinc'

--->> compute the BAZ-interval binsums

'compute RFs binned with epicentral distance: enter back-azimuth limits ib1,ib2 (integers!)'
' do you want to change 'epicentral spacing? (1=yes)'
'enter delta range and spacing (ep1,ep2,dep)'


Then Computer code takes off, remembering the binsum parameters for successive target depths.

OUTPUT

grid files for short interval of time around each target depth NNN-km

out[rt]_baz_NNN.grid, out[rt]_epi_NNN.grid

grid file for spliced baz- and epicen-sweeps

ous[rt]_baz.grid, ous[rt]_epi.grid

You can plot these files with the posted GMT scripts, after adjusting the BOX variable to change plot limits.

c  program recfunk_mwms

c

c  08/21/07 JJP  --- adapted from recfunk_mwm.f

c  updated 07/27/15 JJP

c  MIGRATION/MOVEOUT variant to compute cross-correlation in a moving time window

c  you tell the program the depth range in which you want to focus the cross-correlation

c  and the depth spacing

c  and it computes the moving-window delay from an input model

c  using the epicentral distance to look up the horizontal slowness

c  the program loops through the data set N times for N depth levels and

c  splices the binsummed RFs

c  

c  multitaper program to invert for receiver functions

c  SINGLE STATION WITH SPECTRAL WEIGHTING

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

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

c

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_NNN.grid out[rt]_epi_NNN.grid 

c  where NNN corresponds to the target depth in kilometers

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c

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

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

c