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

recfunk_estack.f -- an event-stacking algorithm for RF with long delay times

recfunk_estack.f -- an event-stacking algorithm for RF with long delay times

Click the title to access the Fortran code recfunk_estack.f in

recfunk_estack is a version of the multiple taper correlation receiver-function estimator that allows the time-domain RF to extend to times longer than (roughly) T/10, where T is the duration of the analysis window. Recfunk_estack does NOT apply a moveout or migration correction, so the ability to resolve Ps from deeper interfaces is limited. Other codes to be posted soon will have these properties.

The major shortcoming of the MTC RF estimator of Park & Levin (2000) is its difficulty retrieving Ps conversions at long delay times. Receiver function studies of upper mantle structure at the base of the lithosphere and the transition zone involve Ps converted phases with delay times between 20 and 70 seconds (e.g., Bostock, 1998; Sheehan et al 2000), but the frequency averaging properties of the Slepian data tapers with time-bandwidth product p\ge 2 tends to limit typical RFs to delay times \tau<20 br="" s.="">
Frequency averaging affords the MTC algorithm desirable statistical properties, especially for single-event RF estimates. We can achieve a narrower frequency average for the retrieval of Ps from interfaces deeper than 100 km if we sum single spectral estimates over K seismic events, rather than K eigenspectra from one event, in the RF cross-correlation formulas (2) and(3). Using only one spectrum estimate per component for each P wave, we estimate the spectra with a Slepian taper of time-bandwidth p=1, that is, the "one-\pi-prolate" taper. This taper averages local frequency information over a central bandwidth similar to the simple boxcar taper, but is optimized for spectral leakage resistance.

c  program recfunk_estack

c  5/11/01 JJP -->  REV2 5/16/01

c  updated 08/17/07 JJP, updated for GFORTRAN compiler 07/28/15

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

c  xf77 -o /Users/jjpark/bin/recfunk_estack recfunk_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  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

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

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c

recfunk_pick.f -- version of multiple-taper RF estimation

recfunk_pick.f -- version of multiple-taper RF estimation

click title to obtain code from Google Drive directory JParkCodes

This code is a simple variant of the multiple-taper RF estimator recfunk.f, which has been public for years on my Yale website. This version reads the start time of the P-wave analysis window from the "AMARKER" entry of the SAC header. It also rotates the vertical and radial seismogram components into LQT coordinates using a formula based on epicentral distance. This formula for LQ rotation angle is very crude, but simple measures of the actual rotation angle scatter quite a bit, and every station has unique behavior. This variability suggests that careful derivation of a precise LQ rotation angle is a fool's game.  In this code there is a simple formula used, but most of my other RF codes use the slowness of the incoming P wave to estimate the rotation angle.  This choice works very well for synthetic seismograms, but real data can be diverse in behavior.

c  program recfunk_pick -- UPDATED
c  updated for GFORTRAN compiler 07/28/15
c
c  original code written for the paper
c  Park, J., and V. Levin, Receiver functions from multiple-taper
c  spectral correlation estimates, Bull. Seismo. Soc Amer., v90,
c  pp1507-1520, 2000.
c
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  the code reads filenames from a file that defaults to "in_recpick"
c  the user inputs the length of the data window to be processed after the
c  time marker in the SAC header.
c
c  08/16/07 JJP
c  multitaper program to compute receiver functions
c  SINGLE STATION WITH SPECTRAL WEIGHTING
c
c  uses spectral estimates of pre-event noise
c  to weight the inversion in the frequency domain, thereby avoiding the
c  usual problem with microseism noise in 0.1-0.5 Hz
c  IN PRACTICE, however, the regularization of what might be a spectral division is achieved
c  via multiple taper cross-correlation
c
c  Program has extra code to perform chi-squared tests, discussed in Park and Levin paper
c
c  xf77 -o /park/backup/bin/recfunk_pick recfunk_pick.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c  xf77 -o /Users/jjpark/bin/recfunk_pick recfunk_pick.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c

recfunk.f -- original recfunk version without moveout corrections

recfunk.f -- original recfunk version without moveout corrections

Click the title of the post to access the code recfunk.f on Google Drive directory JParkCodes.

c  program recfunk

c  5/25/99 JJP updated for GFORTRAN compiler 07/27/15

c

c  original code written for the paper 

c  Park, J., and V. Levin, Receiver functions from multiple-taper 

c  spectral correlation estimates, Bull. Seismo. Soc Amer., v90, 

c  pp1507-1520, 2000.

c

c  this version of the RF code reads datafile names from standard input

c  then the time intervals in the filename file

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

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                <-- 399="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">

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  multitaper program to invert for receiver functions

c  SINGLE STATION WITH SPECTRAL WEIGHTING

c  

c  uses spectral estimates of pre-event noise 

c  to weight the inversion in the frequency domain, thereby avoiding the 

c  usual problem with microseism noise in 0.1-0.5 Hz

c

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

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

rfmigrate.f -- recfunk version with migration

rfmigrate.f -- recfunk version with migration

Click the title of the post to access the code rfmigrate.f on Google Drive directory JParkCodes.

c  program rfmigrate

c  9/20/00 JJP -- adapted from recfunk  --- UPDATED 08/19/07 JJP

c   updated again 07/22/15 JJP

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, because anirec_synth was updated 

c  in 2015 to have a proper-length pre-event time window.

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  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 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 oum[rt]_epi.grid oum[rt]_baz.grid

c

c  out --> no moveout correction;  oum --> moveout correction

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

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

c