rfmig_mcboot.f -- code to compute complex-valued RF harmonic expansion in back-azimuth, moving-window moveout correction

rfmig_mcboot.f -- code to compute complex-valued RF harmonic expansion in back-azimuth, moving-window moveout correction

click title to see Google Drive directory JParkCodes

c  program rfmig_mcboot

c  10/12/04 JJP -- adapted from rfmigrate 

c  07/27/15 JJP

c  does radial/transverse phase-shift stack

c

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

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

c

c  computes moving-window moveout correction for MTC receiver functions 

c  applied 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  code computes frequency-domain stacks of receiver functions that follow 

c  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

c  amplitude dependence.  

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

c  The 2-lobed terms govern dipping interface effects and tilted 

c  symmetry-axis ansotropy. 

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

c    The code regresses for the harmonic expansion, using combined 

c   radial/transverse stack bootstrap-resamples the data to estimate the 

c    uncertainty of the harmonic terms.

c  

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

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

rfmig_mboot.f -- code to compute RF harmonic expansion in back-azimuth, moving-window moveout correction

rfmig_mboot.f -- code to compute RF harmonic expansion in back-azimuth, moving-window moveout correction

click title to see Google Drive directory JParkCodes

c  program rfmig_mboot

c  10/12/04 JJP -- adapted from rfmigrate

c

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

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

c

c  computes moving-window moveout correction for MTC receiver functions 

c  applied 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  code computes frequency-domain stacks of receiver functions that follow a 

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

c  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 the data 

c    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  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

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

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

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

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  

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

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

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

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

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

recfunk_mwmj.f -- moving-window moveout RF estimation with jackknife uncertainty on the bin-sums

recfunk_mwmj.f -- moving-window moveout RF estimation with jackknife uncertainty on the bin-sums

click title to access the code on Google Drive

c  program recfunk_mwmj

c  6/3/04 JJP  --- UPDATED 08/20/07 JJP, and 07/27/15 JJP

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

c  This version performs a jackknife uncertainty estimate of bin-summed RFs

c

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

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

c  using the epicentral distance to look up the horizontal slowness

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:

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  jackknife uncertainty envelope (lower) out[rt]1_baz_NNN.grid out[rt]1_epi_NNN.grid

c  jackknife uncertainty envelope (upper) out[rt]2_baz_NNN.grid out[rt]2_epi_NNN.grid

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_mwmj recfunk_mwmj.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a

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

c

recfunk_mwm.f -- moving-window moveout

recfunk_mwm.f -- moving-window moveout

click the title of the post to access the code.

This code applies a moveout correction by shifting the time window of the horizontal-component time series before DFT and multiple-taper spectral correlation. The code prompts the user to specify a stacking model in the anirec format, and a target depth. The code uses the phase type (P,PP or PKP) and epicentral distance to compute a predicted time that a Ps phase from the target depth should arrive. The horizontal time series is shifted by this time, so that the RF at t=0 is the Ps pulse at the target depth. The zero-time pulse on the radial RF associated with the direct P wave is shifted to negative times, and will often lose focus in back-azimuthal stacks. The GMT script for plotting is similar to that used for other recfunk codes, but has a shifted X-axis to focus the eye on the signals at Ps delay-times both before and after the target Ps arrival.

The code writes output files that flag the target depth:

outr_baz_040.grid,outt_baz_040.grid,outr_epi_040.grid,outt_epi_040.grid

where the integer 040 stands for a target depth of 40 km in the stacking model. GMT scripts for plotting are executed this way:

./Plot_baz outr_baz_040.grid outt_baz_040.grid 1. RAYN_mwm

./Plot_epi outr_epi_040.grid outt_epi_040.grid 1. RAYN_mwm

c  program recfunk_mwm

c  6/3/04 JJP  --- UPDATED 08/19/07 JJP

c   updated again 07/22/15 JJP

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

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

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

c  using the epicentral distance to look up the horizontal slowness

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:

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

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_mwm recfunk_mwm.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a

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

c

rfmigrate_estack.f -- code to migrate the RFs in the frequency domain

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