Click the link to access the code rfmigrate_estack.f
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).
For GMT scripts to plot these output files, see blog posts from July.
Saturday, August 18, 2007
ttimes_PKP.dat -- accessory file for RF migration and moveout corrections
click title to see Google Drive directory JParkCodes
This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for PKP/PKIKP. Numbers massaged from ttimes, the public code available from IRIS. Click the title to access the file.
This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for PKP/PKIKP. Numbers massaged from ttimes, the public code available from IRIS. Click the title to access the file.
ttimes_P.dat -- accessory file for receiver function moveout correction
click title to see Google Drive directory JParkCodes
This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for P and Pdiff. Numbers massaged from ttimes, the public code available from IRIS.
This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for P and Pdiff. Numbers massaged from ttimes, the public code available from IRIS.
Friday, August 17, 2007
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.
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.
Thursday, August 16, 2007
Modified recfunk_pick_through.f
click title to see Google Drive directory JParkCodes
More options for specifying the start time of the data record are recognised.
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
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Ritz/jlib.a
c
......implicit real*4(a-h,o-z)
......real*8 ar,ai,el(12),aaa
......complex*8 zc,zero
......complex*8 afft,rf,rfs
......character*80 name,subtext,string
......character*28 title(2)
......character*10 cmps(3)
......character comp(3)
......common/nnaamm/iwflag
......common/npiprol/anpi
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......dimension iah(158),chead(158),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data comp/'r','t','z'/
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
......open(10,file='in_recpick_raw',form='formatted')
......open(11,file='in_recpick',form='formatted')
10 print *,'enter filename e.g. 1998.361.bhz'
......read(10,101) name
101 format(a)
......if(name(1:4).eq.'stop') go to 111
......subtext=name
c we assume that the last character is either 'r','t',or 'z'
......do i=80,1,-1
......if(name(i:i).ne.' ') then
......namlen=i
......go to 98
......endif
......end do
......print *,name
......print *,'no filename?'
......stop
98 continue
......do kcmp=3,3
......name(namlen:namlen)=comp(kcmp)
......print *,name
......call sacread(name,a(1,kcmp),ierr)
......if(ierr.eq.1) stop
......dt=ahead(1)
......tt(1)=0.
......tt(2)=dt
......tt(3)=dt
......nscan=iah(80)
......if(nscan.gt.99000) then
......print *,'Careful! data length greater than 99000:', nscan
......pause
......endif
......end do
c end kluge
c if(baz.lt.70) then
c call plotit_axes(0.,0.,0.,0.)
c call manyplot(99000,1,nscan,3,tt,a,1.1,name,'time(sec)',' ')
c endif
c plot kluge
c KLUGE: CHANGES FOR READING TIME INTERVAL FROM SAC HEADER
c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk
c ahead(12,13,14) are SAC's T1,T2,T3 parameters, for arrival time. Set by pressing "T" and "N" in plotpk
......if(ahead(9).gt.1.0.or.ahead(12).gt.1.0.or.ahead(13).gt.1.0.
......x or.ahead(14).gt.1.0) write(11,101) subtext
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart,
......x (ahead(n),n=12,14)
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end
More options for specifying the start time of the data record are recognised.
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
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Ritz/jlib.a
c
......implicit real*4(a-h,o-z)
......real*8 ar,ai,el(12),aaa
......complex*8 zc,zero
......complex*8 afft,rf,rfs
......character*80 name,subtext,string
......character*28 title(2)
......character*10 cmps(3)
......character comp(3)
......common/nnaamm/iwflag
......common/npiprol/anpi
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......dimension iah(158),chead(158),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data comp/'r','t','z'/
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
......open(10,file='in_recpick_raw',form='formatted')
......open(11,file='in_recpick',form='formatted')
10 print *,'enter filename e.g. 1998.361.bhz'
......read(10,101) name
101 format(a)
......if(name(1:4).eq.'stop') go to 111
......subtext=name
c we assume that the last character is either 'r','t',or 'z'
......do i=80,1,-1
......if(name(i:i).ne.' ') then
......namlen=i
......go to 98
......endif
......end do
......print *,name
......print *,'no filename?'
......stop
98 continue
......do kcmp=3,3
......name(namlen:namlen)=comp(kcmp)
......print *,name
......call sacread(name,a(1,kcmp),ierr)
......if(ierr.eq.1) stop
......dt=ahead(1)
......tt(1)=0.
......tt(2)=dt
......tt(3)=dt
......nscan=iah(80)
......if(nscan.gt.99000) then
......print *,'Careful! data length greater than 99000:', nscan
......pause
......endif
......end do
c end kluge
c if(baz.lt.70) then
c call plotit_axes(0.,0.,0.,0.)
c call manyplot(99000,1,nscan,3,tt,a,1.1,name,'time(sec)',' ')
c endif
c plot kluge
c KLUGE: CHANGES FOR READING TIME INTERVAL FROM SAC HEADER
c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk
c ahead(12,13,14) are SAC's T1,T2,T3 parameters, for arrival time. Set by pressing "T" and "N" in plotpk
......if(ahead(9).gt.1.0.or.ahead(12).gt.1.0.or.ahead(13).gt.1.0.
......x or.ahead(14).gt.1.0) write(11,101) subtext
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart,
......x (ahead(n),n=12,14)
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end
recfunk_pick.f -- update to allow the identification of P phase in the SAC header
click title to see Google Drive directory JParkCodes
Update suggested by Vadim Levin
recfunk_pick.f
c program recfunk_pick -- UPDATED
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 08/16/07 JJP
This code also asks the user whether he/she wishes to apply the LQT coordinates (not everyone does!)
When you prepare P-wave files using the previously blogged procedures, or with your own algorithm, you can select "T" & "1" at the start of a P wave or Pdiff, "T" & "2" at the start of a PKP/PKIKP window, and "T" & "3" at the start of a PP window. The code change for recfunk_pick creates an epicentral distance variable to determine the LQ-rotation angle and (in other codes) the moveout/migration parameters. GCARC --> GCARC/2 for the PP phase, because PP has two equal P-phase legs.
Update suggested by Vadim Levin
recfunk_pick.f
c program recfunk_pick -- UPDATED
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 08/16/07 JJP
This code also asks the user whether he/she wishes to apply the LQT coordinates (not everyone does!)
When you prepare P-wave files using the previously blogged procedures, or with your own algorithm, you can select "T" & "1" at the start of a P wave or Pdiff, "T" & "2" at the start of a PKP/PKIKP window, and "T" & "3" at the start of a PP window. The code change for recfunk_pick creates an epicentral distance variable to determine the LQ-rotation angle and (in other codes) the moveout/migration parameters. GCARC --> GCARC/2 for the PP phase, because PP has two equal P-phase legs.
Friday, August 3, 2007
Vadim Levin links to RF codes
From my good friend and colleague Vadim Levin:
hey, cool blog :)
why don't you put a plug for my version of the RF software there,
it's at
http://www.ldeo.columbia.edu/~vadim/RF/RF-manual.html
hey, cool blog :)
why don't you put a plug for my version of the RF software there,
it's at
http://www.ldeo.columbia.edu/~vadim/RF/RF-manual.html
Subscribe to:
Posts (Atom)