Sunday, October 16, 2016

New Versions of rfmig_mcboot.f, rfmig_boot.f, rfmig_mboot.f, rfmig_cboot.f, rfmig_mcbp.f, rfmig_mcbf.f, all with improved handling of larger data sets (500+ events).


Avoided some NaNs in the harmonic-regression calculations.  Earlier versions of the codes included some single-precision computations that led the algorithm into darkness.  Tested with a 1000-record dataset.  The link takes you to the Google Drive folder where the codes reside.

Wednesday, June 15, 2016

Picking P waves, PKP waves, PP waves, etc for Recfunk codes

c program recfunk_pick_through
c version of recfunk_pick that reads seismic record start times from
c the sac header, and writes to in_recpick the filenames that have been picked
c using EITHER the A-marker or the T1, T2, or T3 markers in the SAC header, placed by running plotpk
c within SAC. Set by pressing "T" and 1,2 or 3 in plotpk sequentially
c 10/13/05 JJP
c input file: in_recpick_raw
c output file: in_recpick

Tuesday, December 8, 2015

Converting JWeed event files to julian days (new for 2015)


JWeed event files identify earthquakes in month-day format, but seismic data filenames follow Julian-day format. Unless you have a magic patch of brain cells that converts the calendar to julian days, this is a big hassle when matching files with events. I adapted a little fortran code to convert, and some awk lines that convert a JWeed event file into a format that is used by other codes to pre-process body wave data. Here are the awk lines:

sed 's/,/ /g' events*.events | sed 's[-[ [g' >! blob.dat
awk '{print "echo " $3 " " $4 " " $2 " >! blobb; echo " $0 " >>blobb ; Julian1 < blobb >> ev_julian " }' blob.dat >! c_events

IF THE EVENT INFO IS IN A JWEED SUMMARY FILE, start with extracting the EVENT lines, then remove the dashes, but only the first two (NOT the minus signs on latitude and longitude!), then sort the events from earliest to latest, then generate a command file to invoke the Julian1 utility.  Remember to change the first line of command file c_events (from >> to >!) to create the file ev_julian.

grep EVENT RAR_events.summary >! RAR.events
sed 's/,/ /g' RAR.events | sed 's[-[ [' | sed 's[-[ [' >! blobb.dat
sort blobb.dat >! blob.dat
awk '{print "echo " $3 " " $4 " " $2 " >! blobb; echo " $0 " >>blobb ; Julian1 < blobb >> ev_julian " }' blob.dat >! c_events

optional to remove extra "words" from blob.dat:
sed 's/xxxx yyyyy//g' c_events >! c_ev

cd RAR (from the Teleseismic directory)
ls *.SAC >! infiles_raw
cp ../ev_julian events_rawlist
cp infiles_raw infiles1
events_prune
te events &
te infiles1 &

the JWeed event filenaming convention (for me) is events_YYYY.events, with variants. The first step strips the syntax separating numbers in the events file, so that awk can parse the numbers. The second line creates a command file c_events that uses Julian1 to read the output of the awk step. You must change the ">> ev_julian" in the first line to ">! ev_julian" so that the event file ev_julian is created to write into when you "source c_events"    The optional step takes out some superfluous matter from c_events, then you "source c_ev" to make ev_julian. I rename this file events.julian as a template for the file "events" that each station needs in preprocessing (see previous post).

Oh yes, there is a Fortran code

c program Julian1
c f77 -o /Users/jjpark/bin/Julian1 Julian1.f
character string*80
dimension mmonth(12)
data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
c print *,'enter month, day, year'
read(5,*) jmo,jd,jy
leap=1-(mod(jy,4)+3)/4
jd_julian=mmonth(jmo)+jd
if(jmo.gt.2) jd_julian=jd_julian+leap
c print *,'julian day = ',jd_julian
read(5,101) string
print 102,jy,jd_julian,string
101 format(a)
102 format(2i5,4x,a)
stop
end

Wednesday, July 29, 2015


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

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

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

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

 
Link