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.
Wednesday, August 22, 2007
Monday, August 20, 2007
recfunk_mwmj.f -- moving-window moveout RF estimation with jackknife uncertainty on the bin-sums
c program recfunk_mwmj
c 6/3/04 JJP --- UPDATED 08/20/07 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 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 use the GMT scripts in the previous few blog posts for plotting these files.
c 6/3/04 JJP --- UPDATED 08/20/07 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 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 use the GMT scripts in the previous few blog posts for plotting these files.
GMT shellscript syntax that is not showing up in the blog posts CAVEAT EMPTOR
NOTE -- there is a bit of shellscript syntax that is not showing up in the blog post. There is a shell command
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
GMT script for the moving-window with jackknife uncertainty -- epicentral bin-sums
NOTE -- there is a bit of shellscript syntax that is not showing up in the blog post. There is a shell command
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 180
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 180
END
echo gv is JPark's alias for ghostscript
gv plot.ps
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
7.5 -10
7.5 180
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <
7.5 -10
7.5 180
END
echo gv is JPark's alias for ghostscript
gv plot.ps
GMT script for the moving-window with jackknife uncertainty -- back-azimuth bin-sums
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
gmtset ANOT_FONT Times-Roman HEADER_FONT Times-Roman LABEL_FONT Times-Roman
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END
echo gv is JPark's alias for ghostview
gv plot.ps $
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
gmtset ANOT_FONT Times-Roman HEADER_FONT Times-Roman LABEL_FONT Times-Roman
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
7.5 -10
7.5 380
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <
7.5 -10
7.5 380
END
echo gv is JPark's alias for ghostview
gv plot.ps $
stacking model example
stacking model example (used for GSN Station RAYN) 4 layers over a halfspace
4
50 145
5000 5500 -0.06 0 3070 -0.06 2700
0 0
8000 6350 0 0 3494 0 2700
0 0
21000 6420 0.02 0 3510 0.02 2700
0 0
40000 6930 0.03 0 3804 0.03 2700
0 0
220000 7900 -0.15 0 4300 -0.15 3350
4
50 145
5000 5500 -0.06 0 3070 -0.06 2700
0 0
8000 6350 0 0 3494 0 2700
0 0
21000 6420 0.02 0 3510 0.02 2700
0 0
40000 6930 0.03 0 3804 0.03 2700
0 0
220000 7900 -0.15 0 4300 -0.15 3350
GMT script for plotting epicentral RF sweep from moving-window-moveout grid files
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces
/bin/rm wigtab? tmp? plot.ps
echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-10/10/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )
echo working on $i
awk '$2 == A && $1 >= -10' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -10' A=$i $2 >> wigtabT
#in case you decide to demean...
#awk '$2 == A && $1 >= -2' A=$i $1 > tmp1
#set mean = `awk '{a+=$3} END {print a/NR}' tmp1 `
#echo baz $i mean $mean
#awk '{print $1,$2,$3-mn}' mn=$mean tmp1 >> wigtab
echo ">" >> wigtabR
echo ">" >> wigtabT
end
set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
5 -10
5 180
END
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/255/0/0 -O <> plot.ps
5 -10
5 180
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces
/bin/rm wigtab? tmp? plot.ps
echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-10/10/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )
echo working on $i
awk '$2 == A && $1 >= -10' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -10' A=$i $2 >> wigtabT
#in case you decide to demean...
#awk '$2 == A && $1 >= -2' A=$i $1 > tmp1
#set mean = `awk '{a+=$3} END {print a/NR}' tmp1 `
#echo baz $i mean $mean
#awk '{print $1,$2,$3-mn}' mn=$mean tmp1 >> wigtab
echo ">" >> wigtabR
echo ">" >> wigtabT
end
set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
5 -10
5 180
END
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/255/0/0 -O <
5 -10
5 180
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
GMT script for plotting backazimuthal RF sweep from moving-window-moveout grid files
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_baz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
# set BOX = -R25/65/-5/365
set BOX = -R-10/10/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $3
set title = $4
pswiggle $1 $FRAME $BOX -Z$scale -M -G0/0/255 -P -K > plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SWen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END
pswiggle $2 $FRAME $BOX -Z$scale -M -G0/0/255 -P -O -K -X3.5 >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SwEn -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $2 scale $3 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -1.25 2.7 24 0 5 CT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_baz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
# set BOX = -R25/65/-5/365
set BOX = -R-10/10/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $3
set title = $4
pswiggle $1 $FRAME $BOX -Z$scale -M -G0/0/255 -P -K > plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SWen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
7.5 -10
7.5 380
END
pswiggle $2 $FRAME $BOX -Z$scale -M -G0/0/255 -P -O -K -X3.5 >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SwEn -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $2 scale $3 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -1.25 2.7 24 0 5 CT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <
7.5 -10
7.5 380
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
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
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
Sunday, August 19, 2007
rfmigrate.f -- recfunk version with migration
Click the title of the post to access the code rfmigrate.f on Google Drive directory JParkCodes.
Version of MTC RF estimation that computes single-station RFs with 3 2.5-pi Slepian eigentapers, similar to recfunk_pick.f. This code uses a migration scheme to align Ps pulses to the time expected for Ps from a vertical-incidence wave. See comments for the rfmigrate_estack code, posted yesterday (August 18, 2007).
Version of MTC RF estimation that computes single-station RFs with 3 2.5-pi Slepian eigentapers, similar to recfunk_pick.f. This code uses a migration scheme to align Ps pulses to the time expected for Ps from a vertical-incidence wave. See comments for the rfmigrate_estack code, posted yesterday (August 18, 2007).
recfunk_svd.f version of recfunk with SVD in the inner loop
click the title of post to obtain the code recfunk_svd.f in
From unpublished manuscript on RF estimation:
Receiver Functions from Multiple-Taper Spectral Correlation:
Statistics, Single-Station Migration, and Jackknife Uncertainty
J. Park and V. Levin, submitted to GJI
RF-computation suggests that the coherence between P and SH-converted energy in the P coda is comparable to the coherence between P and SV-converted energy. Because aggregate coherence C(f)< 0.5, one can question whether the SV and SH converted phases are mutually correlated. The alternative hypothesis is that SV and SH scattering arise from statistically distinct portions of the incoming P wave. Poor correlation between SV and SH scattering can be tested by simultaneous cross-correlation of all three particle-motion components, computed via the singular-value decomposition (SVD). Define a K ×3 matrix A(f) of eigenspectral estimates in the P-SV -SH coordinate system (in our application K = 3):
For noise-free, perfectly correlated data, the matrix A(f) = u (f) v(f), a single dyad in which the 3-component vector v(f) contains the correlation coefficients between components and the K-component vector u(f) describes the spectral signature of the correlated signal in a narrow frequency band about f. The receiver functions HSV (f),HSH(f) correspond to the component ratios vSV /vP and vSH/vP , where v(f) = [vP , vSV , vSH]. For real noisy data, we associate the SVD receiver function with the singular vectors of A that correspond to the largest singular value, similar to the ”pure state” filtering algorithm (Samson, 1983). We SVD decompose A(f). The uncertainty of the SVD receiver function can be estimated by analogy to uncertainties of RF computed from complex cross-correlation and coherence.
In tests with real data, the added sophistication of the SVD receiver function algorithm leads only to small changes in the estimated RFs. We infer that P-SV conversion and P-SH conversion in the P coda are co-dependent, in a statistical sense. It should therefore be feasible to estimate HSV (f) and HSH(f) jointly. On the other hand, our experiments suggest that the ordinary pairwise MTC RF estimator is adequate in practice. The SVD-based RF estimator may have utility in complex situations, such as the estimation of S receiver functions for S-to-P converted energy.
From unpublished manuscript on RF estimation:
Receiver Functions from Multiple-Taper Spectral Correlation:
Statistics, Single-Station Migration, and Jackknife Uncertainty
J. Park and V. Levin, submitted to GJI
RF-computation suggests that the coherence between P and SH-converted energy in the P coda is comparable to the coherence between P and SV-converted energy. Because aggregate coherence C(f)< 0.5, one can question whether the SV and SH converted phases are mutually correlated. The alternative hypothesis is that SV and SH scattering arise from statistically distinct portions of the incoming P wave. Poor correlation between SV and SH scattering can be tested by simultaneous cross-correlation of all three particle-motion components, computed via the singular-value decomposition (SVD). Define a K ×3 matrix A(f) of eigenspectral estimates in the P-SV -SH coordinate system (in our application K = 3):
For noise-free, perfectly correlated data, the matrix A(f) = u (f)
In tests with real data, the added sophistication of the SVD receiver function algorithm leads only to small changes in the estimated RFs. We infer that P-SV conversion and P-SH conversion in the P coda are co-dependent, in a statistical sense. It should therefore be feasible to estimate HSV (f) and HSH(f) jointly. On the other hand, our experiments suggest that the ordinary pairwise MTC RF estimator is adequate in practice. The SVD-based RF estimator may have utility in complex situations, such as the estimation of S receiver functions for S-to-P converted energy.
Saturday, August 18, 2007
rfmigrate_estack.f -- code to migrate the RFs in the frequency domain
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.
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.
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)