Thursday, December 6, 2007

plotting script for rfmig_cboot.f and rfmig_mcboot.f output WITH UNCERTAINTIES

click title to see Google Drive directory JParkCodes



# ./Plotj_cexp outr_cexp.grid outr1_cexp.grid outr2_cexp.grid outt_cexp.grid outt1_cexp.grid outt2_cexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
# Input - $1 - mean R+T $2 lower R+T $3 upper R+T $4 mean R-T $5 lower R-T $6 upper R-T $7 scale $8 TITLE

plotting script for rfmig_boot.f and rfmig_mboot.f output WITH UNCERTAINTIES

click title to see Google Drive directory JParkCodes


# ../Plotj_bexp outr_bexp.grid outr1_bexp.grid outr2_bexp.grid outt_bexp.grid outt1_bexp.grid outt2_bexp.grid 0.3 title

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
# Input - $1 - mean R $2 lower R $3 upper R $4 mean T $5 lower T $6 upper T $7 scale $8 TITLE

plotting script for rfmig_cboot.f and rfmig_mcboot.f output

click title to see Google Drive directory JParkCodes

#../Plot_cexp outr_cexp.grid outt_cexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
#Input - $1 - file name $2 file name $3 - scale $4 TITLE

plotting script for rfmig_boot.f and rfmig_mboot.f output

click title to see Google Drive directory JParkCodes

#../Plot_bexp outr_bexp.grid outt_bexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_boot.f or rfmig_mboot.f
#Input - $1 - file name $2 file name $3 - scale $4 TITLE

Tuesday, December 4, 2007

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 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 harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. 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, using combined radial/transverse stack
c 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

nboot=0: compute simple regression, no bootstrap computation of variance.

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 nboot=0 -- only compute a single regression for RF harmonic expansion --> no bootstrap uncertainty estimate
c
c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. 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 to estimate the
c 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

Sunday, December 2, 2007

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

click title to see Google Drive directory JParkCodes

c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. 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 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 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 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

Saturday, December 1, 2007

rfmig_boot.f -- code to estimate harmonic expansion of RFs in back azimuth

click title to see Google Drive directory JParkCodes

c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. 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 to estimate the
c uncertainty of the harmonic terms. The posted version only
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 (see previous posts on codes that apply
c a frequency-domain moveout correction for the variation of Ps delay time with epicentral distance.)
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:
c the file is "in_recfunk" and has lines of the form:
c
c 1997.070.19.33.bh? <-- br="" code="" r="" replaces="" t="" with="" z="">c 57 52 <-- analysis="" br="" duration="" of="" sec="" start="" time="" window="">c 1997.076.08.15.bh?
c 62 62
c ...
c ...
c ...
c stop <-- 799="" br="" code="" data="" events="" finished="" is="" max="" tells="" that="">c
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

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

click title to see Google Drive directory JParkCodes

c  anirec_synth_circle - anirec to compute synthetics for a specified model and
c  a specified station location from one seismic record read from in_recfunk
c  the code will take the station location  and generate a ring of 120 events that are 90 degrees
c  away, one for every 3 degrees back-azimuth.  To do this we generate a great circle path
c  with the station at the pole of the great circle
c  modified to accept SACfiles with start time info in the header, in A-marker (ahead(9))
c  because the timing information is not used, the value of A-marker is not important.
c  However, in the synthetics, the A-marker and the T1-marker variables in the header are set to 5.0 seconds
c  for easy analysis by recfunk codes that use the header information for timing.
c  the synthetic files are s_NNN.bh[zrt], where NNN is the back azimuth.
c  the code writes in_recpick_circle, a list of the filenames for recfunk_pick and other recfunk codes that
c  expect the start time in the SAC header

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

click title to see Google Drive directory JParkCodes

c anirec_synth - anirec to compute synthetics for a specified model and
c a specified collection of seismic records (read from in_recpick-format file of filenames)
c modified to accept SACfiles with start time info in the header, in A-marker (ahead(9))
c because the timing information is not used, the value of A-marker is not important.
c However, in the synthetics, the A-marker and the T1-marker (or T2) variables in the header are set to 5.0 seconds
c for easy analysis by recfunk codes that use the header information for timing.
c code divides between T1 (P) or T2 (PKP) markers based on epicentral distance, with 120 degrees the divide
c the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c the code writes in_recpick_synth, a list of the filenames for recfunk_pick and other recfunk codes that
c expect the start time in the SAC header
c

The goal of this code is to generate P-coda form a specified model to replicate the earthquake distribution of a real data set. This allows the analyst to test whether the oddities in an RF sweep are likely caused by an imperfect data distribution, or by a shortcoming of his/her model.

Wednesday, August 22, 2007

recfunk_mwms.f -- moving-window migration that splices together a long RF from a series of target depths

This code can be accessed by clicking the title link

recfunk_mwms.f generates a succession of RF estimates for a succession of migration target depths. It splices together the results into a long RF. I have used this code to reconstruct BAZ-sweeps of RFs from 0 to 75-s delay, encompassing the crust to the 660-km discontinuity.

This code will warm your laptop. For a splice from 0 to 680 depth, targeted every 40 km with 524 events, my Macbook Pro required 15 minutes of run time.

INPUT
'velocity model for stacking (blank = stack_model)'
'depth range and inc for the RFs (km) (dep1,dep2,ddep)'
'enter fmax (Hz)'
'time intervals in file (0) or in SAC-header (1)'
'file of seismic records to read? (space-ret: in_recpick)'
'duration of data windows for analysis'
'minimum number of events in a binsum (eg 1 or 2)'
'rotate to LQT coordinates to isolate P and SV? (0=no)'

--->> then the first pass through the seismic data -- typically for target depth = 0 km.

'change baz-interval or baz-spacing? (1=yes)'
'enter back-azimuth range: bazmax, bazmin, bazinc'

--->> compute the BAZ-interval binsums

'compute RFs binned with epicentral distance: enter back-azimuth limits ib1,ib2 (integers!)'
' do you want to change 'epicentral spacing? (1=yes)'
'enter delta range and spacing (ep1,ep2,dep)'


Then Computer code takes off, remembering the binsum parameters for successive target depths.

OUTPUT

grid files for short interval of time around each target depth NNN-km

out[rt]_baz_NNN.grid, out[rt]_epi_NNN.grid

grid file for spliced baz- and epicen-sweeps

ous[rt]_baz.grid, ous[rt]_epi.grid

You can plot these files with the posted GMT scripts, after adjusting the BOX variable to change plot limits.

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.

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

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

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 $

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

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

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

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

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

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.

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.

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.

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.

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.

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

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.

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

Monday, July 23, 2007

aniprop_downrange.f -- version of aniprop to emulate a linear array of seismometers

click title to see Google Drive directory JParkCodes

c aniprop_downrange - program to calculate propagating modes of anisotropic layer
c writes a file of dispersion curves at evenly-spaced freq points in files out_cvel out_gvel
c also writes a set of SAC_Format data files SYN_nn.BH[ZNE] for downrange distances of
c 75km 150 225 300 ... 1200 1275 1350 1424 1500km downrange
c E-component is the x component (radial) and N-component is the y component (transverse)


http://earth.geology.yale.edu/~jjpark/Code/aniprop_downrange.f

aniprop.f -- a standalone program to generate surface waves in layered anisotropic structures

click title to see Google Drive directory JParkCodes

Sportsfans,

The erstwhile version of aniprop.f that was (is?) available on my website has been superceded by an all-in-one version that computes dispersion curves and synthetic seismograms all in one program. Adapted and debugged on an Intel Mac in gnu fortran.

 aniprop_072307.f  in Google Drive directory JParkCodes

(after a 24 March 2008 bugfix is now aniprop_032408.f)

The code computes surface wave modes at evenly-spaced frequency points between 0 and 0.5 Hz. You can change this by tinkering with the code. It computes synthetics for a source at a specified downrange distance and depth, and plots them using PLOTIT. The dispersion values are output to files out_cvel and out_gvel (phase and group velocities, natch), and the synthetics are written out to SAC-format files.

c aniprop - program to calculate propagating modes of anisotropic layer
c writes a file of dispersion curves at evenly-space freq points
c
c reads a layered model like the following file (ignore the "c "s)
c
c K&H SoCal model, deep crust horizontal anisotropy TITLE
c 3 # OF LAYERS OVER HSPACE
c 45 45 THETA,PHI ORIENTATION ANGLES
c 4000 5500 0.06 0.00 3175 0.03 2600 FOR SYMMETRY AXIS
c 0 0
c 27400 6300 0.00 0.00 3637 0.00 2800 DEPTH (M), VP (M/S), "B", "C",
c 90 45 VS (M/S), "E"
c 32400 6800 0.04 0.00 3925 0.02 2900 B,C,E ARE ANISOTROPIC PARAMETERS
c 0 0
c 60000 7800 0.00 0.00 4500 0.00 3200 NOTE: HSPACE MUST BE ISOTROPIC
c

One quirk of this code is misbehavior for isotropic models. The code solves for Rayleigh and Love dispersion simultaneously, because the two polarizations are hybridized in anisotropic models. Paradox: for a perfectly isotropic model, the rootfinder sometimes misses roots, typically at occasions where the Rayleigh and Love dispersion curves cross. The bookkeeping for solutions is not very rigorous, so loss of a mode can happen. With anisotropy, the hybridization of Rayleigh/Love motion seems to repel the roots and keep some minimal separation. As far as I can tell, the total effect on actual synthetic seismograms is negligible.

Wave propagation in a truly 1-D lossless model leads to some interesting reverbarations that one rarely sees in nature. With all the overtones included, you can see a succession of bouncing S waves that precede the main surface wave pulse. In Earth's crust, such phases would scatter away their high frequency energy from surface topography and sedimentary layers. There is a loop for mode-summing that can be altered to eliminate the overtones, leaving the fundamental Love and Rayleigh waves: (lines 681-3)

c loop over overtones at a certain frequency
do iov=1,maxbr
c do iov=1,2


Notes on the angle conventions for w-hat, the axis of symmetry:

In the anisotropic reflectivity code, subroutine zzget *assumes* a coordinate system in which z is down (anti-vertical), x is the radial direction, and y is anti-transverse. Therefore, the position angles theta,phi for w-hat are tilt relative to down, and azimuth defined as a rotation from x towards y. This rotation is CCW if viewed from below, and CW if viewed from above. Since w-hat and -(w-hat) define the same axis of symmetry, the position angles *also* can be defined as theta=(tilt from vertical) and phi=(rotation from anti-x (anti-radial) towards anti-y (transverse)). Viewed from above, this phi rotation is CW, and defines the strike of w-hat relative to the arrival azimuth of the wave.

In order to compute seismograms for a variety of back-azimuths, the synthetic code genrecd_az.f accepts a layered model with w-hat position angles defined as theta=(tilt from vertical) phi=(strike CW from N). For an event at back-azimuth psi (CW from N), the code rotates w-hat from geographic coordinates to ray-based coordinates before passing it to subroutine zzget. If a wave arrives at back-azimuth psi, the strike of the axis of symmetry w-hat relative to its arrival azimuth is phi'=phi-psi. The code performs this rotation with this code in subroutine matget, for w-hat azimuth "az":

caz=dcosd(az)

saz=-dsind(az) ! sin(-az)

do n=1,nlp

ww(3)=w(3,n)

ww(1)=w(1,n)*caz-w(2,n)*saz

ww(2)=w(1,n)*saz+w(2,n)*caz

...

In this manner, the axes of symmetry of the model, saved in array w(.,.), are never modified.

Journal Reference:

Park, J., Surface waves in layered anisotropic structures, Geophys. J. Int. v126, 173-184, 1996.

Saturday, July 21, 2007

GMT script for plotting epicentral RF sweep from recfunk *.grid files

GMT script for plotting epicentral RF sweep from recfunk *.grid files

Anyone planning to use this post should know what GMT commands are. If not, google "Generic Mapping Tools"
This shell script is based on an Ur-script written by Vadim Levin

#!/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-2/15/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 >= -2' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -2' A=$i $2 >> wigtabT
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

GMT script for plotting backazimuthal RF sweep from recfunk *.grid files

Anyone planning to use this post should know what GMT commands are. If not, google "Generic Mapping Tools"
This shell script is based on an Ur-script written by Vadim Levin

#!/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. ./Plotbaz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles

echo sample usage
echo "./Plotbaz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-2/15/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

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

open plot.ps &
# mv plot.ps plot.eps
# qlmanage -p plot.eps &
# convert -density 300 plot.ps plot.jpg &
echo convert -density 300 plot.ps plot.jpg


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

click title to see 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.

recfunk_pick_072107.f

(you didnt think that I would post a 1000-line Fortran code!)

The code needs the EISPACK library for matrix eigenvalue subroutines. Also Plotit and Sacread, which are the subject of prior posts.

Code is superceded on 8/16/07 to diversify the P-phase identification in the SAC header -- see later post

recfunk_pick_through.f -- read list of data files and retain only those which have a P-wave window marked for RF estimation

click title to see Google Drive directory JParkCodes

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 10/13/05 JJP
c input file: in_recpick_raw
c output file: outfile (to be renamed in_recpick)
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /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/stap2/tap(16400,16)
......common/taperz/ta(16400,8,3),tai(16400,8,3)
......common/stap0/ar(16400),ai(16400),bb(16400),time(99000)
......common/specs/afft(4100,8,3),rf(4100,2),sss(4100,6),crf(4100,2)
......common/saveit/rft(4100,799,2),crft(4100,2),rfs(2050,799,2),
x drfs(2050,799,2),bazs(799),epic(799),s2n(799),ip(799)
......common/saveit2/crfts(4100,799,2)
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......common/distrib/xi(72,36),yi(72,36),sumij(72,36)
......common/chisq/chi(160000),enchi(160000),bbaz(160000),
x chim(200,2),enchim(200,2)
......dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data comp/'r','t','z'/
......data cmps/'Radial ','Transverse','Vertical '/
......data title/'Radial Receiver Function ',
x 'Transverse Receiver Function'/
......con=180./3.14159265358979
......pi=3.14159265358979
......pih=pi/2.
......zero=cmplx(0.,0.)
......fmax=5.
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
c some default values (often superceded)
......fmax0=fmax
......fmin=0
......anpi=2.5
......nwin=3
......npad=16384
......nnf=npad/2
......do i=1,1000
...... crft(i,1)=0.
...... crft(i,2)=0.
......end do
......irec=0
......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 naturally, you are in deep poop if dt is not the same for all records
......npre=10./dt
......npost=30./dt
......ntot=npre+npost
......ttt(1)=-npre*dt
......ttt(2)=dt
......ttt(3)=dt
......do i=1,ntot
...... time(i)=-ttt(1)-(i-1)*ttt(2)
......end do
......irec=irec+1
......if(irec.gt.799) then
...... print *,'max records exceeded'
...... stop
......endif
......baz=ahead(53)
......bazs(irec)=ahead(53)
......epic(irec)=ahead(54)
......print *,baz,'= back azimuth'
c ddf is cycles/sec
......ddf=1./(dt*npad)
......nf1=1
......nf2=fmax/ddf+1
......nf=nf2-nf1+1
......if(nf.le.0) stop
......fr(1)=0.
......fr(2)=ddf
......fr(3)=ddf
c 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......print *,'P window start time, duration (sec)'
c......read(10,*) pstart,postwin
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart
......if(ahead(9).gt.0.1) write(11,101) subtext
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end

Fortran subroutine for X-window plotting of simple arrays

click title to see Google Drive directory JParkCodes

see the following files on Google Drive directory JParkCodes

plotit.html  and plotit.tar

examine_sac.f -- write a SAC macro to visually assess your body-wave data before receiver function estimation

c program examine_sac
c 7/8/99 JJP - adapted for WEED output, after running event_prune & rotate_sac
c 4/1/02 altered to order the input filenames by epicentral distance
c...... to facilitate quality control (P and S wave pattern varies smoothly)
c
c xf77 -o /Users/jjpark/bin/examine_sac examine_sac.f /Users/jjpark/Ritz/jlib.a /Users/jjpark/Plotxy/plotlib.a
c to read filenames of data files (bhr,t,z) for a single event
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85
c slat, slon are in radians, sdep is in km, slat=colatitude
c SAC headers are degrees and meters, respectively, and lat is lat
c ---> so must convert
......common/sour/slat,slon,sdep,sol(6),am0,strike,dip,slip
......common/ptitle/bigbanner,name,namo,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/allfiles/nami(3,6000),gcarc(6000),gcord(6000),ip(6000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12),xx(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data dumf/'ae','an','az'/
......con=180./3.14159265358979
......open(7,file='infiles2',form='formatted') ! sacfiles 19*.bh[rtz] 20*.bh[rtz]
......open(18,file='outfile1',form='formatted') ! to be in_recfunk
......open(9,file='outfile',form='formatted') ! to be sac_examine
......open(13,file='outfile3',form='formatted') ! to be infiles3 (for pruning SKS data sets)
......write(18,103) '0.4'
......j=1
10 print *,'starting again'
......do k=1,3
...... read(7,101,end=998) nami(k,j)
......end do
......name=nami(1,j)
......print *,name
......call sacread(name,a,ierr)
......gcarc(j)=ahead(54)
......j=j+1
......go to 10
998 nevents=j-1
......close(7)
......xx(1)=1.
......xx(2)=1.
c order the events by epicentral distance
......call plotit(xx,gcarc,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,21)
......do l=1,nevents
...... blob=200.
...... do m=1,nevents
...... if(blob.gt.gcarc(m)) then
............blob=gcarc(m)
............mm=m
...... endif
...... end do
...... ip(l)=mm
...... gcord(l)=gcarc(mm)
...... gcarc(mm)=200.
......end do
......call plotit(xx,gcord,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,22)
......do j=1,nevents
...... name=nami(1,ip(j))
...... write(18,103) name(1:17),'?'
...... write(18,103)'117 75 75'
......end do
......write(18,103) 'stop'
......close(18)
......write(9,101) 'qdp 2000'
......do j=1,nevents
...... name=nami(1,ip(j))
........write(13,103) nami(1,ip(j))
........write(13,103) nami(2,ip(j))
........write(13,103) nami(2,ip(j))
...... write(9,103) 'read ',name(1:17),'?; chnhdr a 115 '
...... write(9,103)'rtrend; lp c 0.15 np 8; write lowr lowt lowz'
...... write(9,103) 'read ',name(1:17),'?'
...... write(9,103)'rtrend; hp c .5 np 8; write hir hit hiz'
...... write(9,103) 'readhdr lowz;listhdr gcarc baz evla evlo evdp ',
x 'kzdate'
c to pick arrival times for pulses, change hi?;plot1;pause
c to hi?;plotpk markall on; write over
...... write(9,103) 'read ',name(1:17),'? low? hi?;plot1;pause'
......end do
999 close(9)
......close(13)
101 format(a)
103 format(7a)
102 format(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f4.1)
105 format(22x,f8.2,1x,f9.2,1x,f4.0)
104 format(a,f10.3,a,f8.3,a,f10.2)
......print *,'to pick arrival times, change '
......print *,'to '
......stop
......end

subroutines sacread and sacout -- SAC-binary I/O

These subroutines lack the neat feature of the new SAC release, which can detect and compensate for byte-swapped data. Beware the swapped byte! There are variants of the routines included that swap bytes, dogswork to be sure.

c......subroutine sacread_noswap(name,a,ierr)
......subroutine sacread(name,a,ierr)
c read a SAC-format file, dont swap the bytes
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
...... end do
......endif
c......print *,'rec=',nrec,' Im trying to read the last partial record'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
......subroutine sacread_swap(name,a,ierr)
c......subroutine sacread(name,a,ierr)
c read a SAC-format file -- swap bytes for Intel mac for all numbers, not characters
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhead(j,i)=chd(5-j:5-j)
........end do
......end do
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
c swap bytes
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a(i)=ahhd
......end do
c......print *,'rec=2'
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
c...... print *,'rec=',irec
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
c swap bytes
...... do i=1,128
............ahd=a((irec-3)*128+98+i)
...... .....do j=1,4
..............chhd(j:j)=chd(5-j:5-j)
........ end do
........ a((irec-3)*128+98+i)=ahhd
...... end do
...... end do
......endif
......print *,'byteswapping read'
c this here is a big fat kluge --- I cant read the last partial record
c so I skip the last fragment and adjust the header to discard the points.
c......print *,'NOT READING ',last,' points in last record -- SACBUGFIX'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
c swap bytes
......do i=1,last
...... ahd=a((nrec-3)*128+98+i)
...... .do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a((nrec-3)*128+98+i)=ahhd
......end do
c......nscan=nscan-last
c......iah(80)=nscan
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
c......subroutine sacout_noswap(name,a)
......subroutine sacout(name,a)
c write a SAC-format file
......real*4 ahead
......character*80 name
......character*4 chead(158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead)
......open(8,file=name,access='direct',recl=512)
c write the 158-word header
......write(8,rec=1)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
......write(8,rec=2)(ahead(i+128),i=1,30),(a(j),j=1,98)
......do irec=3,nrec
...... write(8,rec=irec)(a((irec-3)*128+98+i),i=1,128)
......end do
......print *,nrec,' records written'
......close(8)
......return
......end
c......subroutine sacout(name,a)
......subroutine sacout_swap(name,a)
c write a SAC-format file
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158),output(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......open(8,file=name,access='direct',recl=512)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
c write the 158-word header
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(i)=ahhd
......end do
......do i=121,128
...... output(i)=ahead(i)
......end do
......write(8,rec=1)(output(i),i=1,128)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
c swap bytes
......do i=1,30
...... output(i)=ahead(128+i)
......end do
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(30+i)=ahhd
......end do
......write(8,rec=2)(output(i),i=1,128)
......do irec=3,nrec
...... iii=(irec-3)*128+98
...... do i=1,128
...... ahd=a(iii+i)
........ do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........ end do
........ output(i)=ahhd
...... end do
...... write(8,rec=irec)(output(i),i=1,128)
......end do
......print *,nrec,'byteswapped records written'
......close(8)
......return
......end

rotate_sac.f -- Preprocessing of SAC data for receiver function analysis

c program rotate_sac_puteventinheader  (sac_rotate assumes the event was already inserted)
c 7/8/99 JJP -- adapted for JWEED output
c f77 -o /Users/jjpark/bin/rotate_sac rotate_sac.f /Users/jjpark/Ritz/jlib.a
c to read infiles1=filenames of data files (E,N,Z) for a single event
c cross-ref with file *events* of source locs and times
c writes a SAC macro to read three components of data, synchronize times, cut to equal length,
c put the event location in the header, and rotate the horizontals to radial and transverse.
c
c there must be three lines in infiles1 for every line in events, and the events must match the data
c take out the dots and replace with spaces before compiling
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85,chh*51
c slat, slon are in radians, sdep is in km, slat=colatitude
c SAC headers are degrees and meters, respectively, and lat is lat
c ---> so must convert
......common/sour/slat,slon,sdep,sol(6),am0,strike,dip,slip
......common/ptitle/bigbanner,nami(3),namo,name,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data dumf/'ae','an','az'/
......con=180./3.14159265358979
......print *,'open(7,file=infiles1) ! sacfilenames'
......print *,'open(18,file=events) ! source locations and times'
......open(7,file='infiles1',form='formatted') ! sacfilenames
......open(18,file='events',form='formatted') ! source locations and times
......open(9,file='outfile',form='formatted')
......write(9,101) 'qdp 2000'
10 print *,'starting again'
c......if(jjc.eq.1) go to 112
......read(18,101,end=110) name
......print *,name
c......go to 10
c 112 continue
......read(name,102) jy,jm,jd,jh,jmi,ssec
......print *, jy,jm,jd,jh,jmi,ssec
c......nami(1:38)=name(43:80)
c......read(name,105) slat,slon,sdep
......read(name(44:80),*) slat,slon,sdep
......print *, slat,slon,sdep
101 format(a)
103 format(7a)
102 format(19x,i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f6.3)
c here reads to column 45 -- if evdp can be read OK here, doesnt matter
c if the spacing of columns is perfect
105 format(22x,f8.2,f9.2,f4.0)
c 105 format(25x,3f8.2)
104 format(a,f10.3,a,f8.3,a,f10.2)
......leap=1-(mod(jy,4)+3)/4
......jd_julian=mmonth(jm)+jd
......if(jm.gt.2) jd_julian=jd_julian+leap
......write(9,103) 'cuterr fillz'
......sdep=sdep*1000.
......print *,slat,slon,sdep,' is depth in meters?'
......do k=1,3
...... read(7,101) nami(k)
......end do
......write(9,103) '/bin/rm ae an az'
......do k=1,3
...... name=nami(k)
...... ii=80
...... do while (name(ii:ii).eq.' ')
...... ii=ii-1
...... end do
c...... ii=ii-4 ! earlier SEED-output name convention
...... ii=ii-6 ! newer SEED-output name convention
...... print *,ii
...... if(name(ii:ii).eq.'Z') then
...... j=3
...... elseif(name(ii:ii).eq.'E'.or.name(ii:ii).eq.'2') then
...... j=1
...... elseif(name(ii:ii).eq.'N'.or.name(ii:ii).eq.'1') then
...... j=2
...... else
...... print *,'bad character: ',name(ii:ii)
...... endif
...... write(9,103) 'sc ln ',name(1:ii+7),dumf(j)
......end do
c LOVROK is the flag for allowing the user to overwrite a data file
......chh='synch;rmean;rtrend;p1;chnhdr LOVROK TRUE;write over'
......write(9,103)'read ae an az;',chh
......write(9,103) 'cut 0 600'
......write(9,103)'read;write over;cut off'
......do k=1,3
...... name=nami(k)
c read in the file to obtain its reference time
...... call sacread(name,a,ierr)
...... if(ierr.ne.0) go to 111 ! in case of mistyped filename
...... iy=iah(71)
...... id=iah(72)
...... ih=iah(73)
...... im=iah(74)
...... iss=iah(75)
...... ims=iah(76)
c compute origin time - year, month are assumed equal
...... torg=ssec-iss-float(ims)/1000.+
x......60.*(jmi-im+60*(jh-ih+24*(jd_julian-id)))
...... write(9,103) 'readhdr ',nami(k)
...... write(9,104) 'chnhdr evla ',slat,' evlo ',slon,' evdp ',sdep
...... write(9,104) 'chnhdr o ',torg
...... write(6,104) 'chnhdr o ',torg
...... write(9,103) 'writehdr'
......end do
......write(9,103) 'read ae an az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......write(9,103) 'read ae an'
......write(9,103) 'rotate to gcp reversed'
......write(9,103) 'read more az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......namo=name(1:15)
......write(9,103) 'write ',namo,'bhr ',namo,'bht ',namo,'bhz'
......go to 10
110 close(7)
......close(18)
......close(9)
......stop
111 print *,'error! ',name
......stop
......end

Synthetic Seismogram normalizations

I have received queries about constants in the synthetic seismogram code anirec.f

In order to avoid round-off error and numerical overflow in numerical computations, it is conventional to nondimensionalize physical units in synthetic seismogram calculations. My codes use the convetion that I inherited from Freeman Gilbert, my dissertation advisor, and the UCSD Free Oscillation crowd. The rules for normalizing Mass, Length and Time are derived from three constraints.

density is normalized by rho-bar = average Earth density = 5515 kg/m^3

length is normalized by the mean Earth radius R=6.371*10^6 m

The last constraint is that pi*G=1, where pi=3.14159265358979d0 and G=6.6723*10^{-11} m^3/(kg-s^2) is the Gravitational constant.

From the public free-oscillation code MINEOS

c the model has been normalised such that a density of 5515 mg/m**3 = 1 ;
c pi*g=1 where g is the grav constant (6.6723e-11) and the radius of the
c earth (rn=6371000 m) is 1. these normalizations result in
c acceleration normalisation = pi*g*rhobar*rn = 7.365081 m/s^2
c velocity normalisation = rn*sqrt(pi*g*rhobar)= vn = 6850.0396 m/s
c frequency normalization = vn/rn = sqrt(pi*g*rhobar) = 1.075190645e-3 s^1


This last constraint leads to a normalization constant for frequency of REN=1.075190645d-3. I have found a slightly different value for G quoted on the web, which is understandable considering that the seismological constants were set when G was not known as precisely as now. There can be danger in mistyping the numbers. One colleague of mine lost days and much sleep when the rho-bar constant was mistyped 5517 kg/m^3 and a 10^{-4} error crept into calculations that should have been good to double precision. Believe me, no one was lauging when the error was discovered!

Converting JWeed event files to julian days

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
sed 's/xxxx yyyyy//g' c_events >! c_ev

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. The last 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

Friday, July 20, 2007

Installing JWeed on Linux

There is a nasty glitch in the Linux installation for JWeed. The IRIS website downloads install.bin, but it has two bad lines. FIRST, run the installer thusly to fix one bad line without changing the checksum

cat install.bin | sed "s/export LD_ASSUME_KERNEL/#xport LD_ASSUME_KERNEL/" > temp
sh temp

Then comment out lines 1335 and 1375 of JWeedv2.9.8, which is the executable script that embodies JWeed. These are the lines

1335 export LD_ASSUME_KERNEL=2.2.5

1375 LD_ASSUME_KERNEL=2.2.5

After these changes, there is still a little complaint about an extra backslash, which I tracked down to line 1343:
split ($3, ver, "[\._-]")

I took the backslash out and the code still worked.


Thanks to Alex Nikulin and Garrett Leahy for the tip!

Preparing RETREAT body waves for Receiver-Function Estimation

Work performed on an IntelMac laptop -- computer codes and scripts will be posted in successive blog posts. This post documents the preparation of body wave data from the RETREAT experiment. You can use the concepts to guide the preparation of data for other stations.

The directory Teleseismic has subdirectories PIIR, VLC, VOLR, etc, one for each station name. The station subdirectories have data files in SAC format, but without the source location in the header. There are three steps in preparing a station's data for processing. First, we associate events with data-file triplets. The extraction by RDSEED from the SEED files should have synchronized the data files for each component of motion. This first step can be tedious, however, because there will be more events than datafile triplets whenever the station was down, and more data files than events when the data files are fragmented. The goal is to prune the files "events" to only the events that have data, and the file "infiles1" to be the triplets of datafile names -- infiles1 will have three times as many lines and events, natch.

cd PIIR (from the Teleseismic directory)
ls *.SAC > infiles_raw
cp infiles_raw infiles1
cp ../events.julian events
te events &
te infiles1 &

events_julian is a master file of events used to request data from the IRIS DMC for the experiment, based on a JWeed event file:

2003 356 MHDF 2003 12 22 19:15:56.000 35.706 -121.102 7.6 MB 6.1 MS 6.4
2003 359 MHDF 2003 12 25 07:11:11.590 8.416 -82.824 33 MB 6.0 MS 6.4
2003 359 MHDF 2003 12 25 20:42:33.720 -22.252 169.488 10 MB 6.3 MS 6.3
2003 360 MHDF 2003 12 26 01:56:52.440 28.995 58.311 10 MB 6.0 MS 6.8
2003 360 MHDF 2003 12 26 21:26:04.100 -22.273 169.314 10 MB 6.1 MS 6.8
2003 361 MHDF 2003 12 27 16:00:59.450 -22.015 169.766 10 MB 6.1 MS 7.1
2003 361 MHDF 2003 12 27 22:38:01.880 -21.672 169.835 10 MB 5.8 MS 6.7

Note that the month/day/year format of the JWeed event file is replaced by year julianday, which matches the output filename format of rdseed.


The user now opens infiles1 and events in a text editor and synchronizes their lines. It is smart to keep an Xterminal open to run SAC or pql. Sometimes there will be more than three datafiles associated with an event start time, and you must display their contents in order to choose the three files to use, or else to merge some files.

pql yyyy.ddd.* &

is a good command for this. Or

echo "read yyyy.ddd.*; qdp 2000; p1" | sac

because you can retrieve previous commands in a Unix shell and replace the yyyy.ddd for different cases.

The second step is to run the code "rotate_sac_puteventinheader" or "rotate_sac." this code reads the info in infiles1 and events. It generates a SAC macro file that will read the data triplets, put the earthquake hypocentre into the SAC headers (NOT for sac_rotate, only for sac_rotate_puteventinheader), cut the files to 600-s length and rotate the horizontals to radial and transverse direction. You run

rotate_sac_puteventinheader
/bin/mv outfile sac_rotate
echo exit >> sac_rotate
sac < sac_rotate

and the files *yyyy.ddd.hh.*BH[ENZ]*.SAC are converted to truncated and rotated data files "yyyy.ddd.hh.mm.bh[zrt]" Now comes the third step of processing. The code examine_sac reads the data files and composes a SAC macro to visualize the data in raw and bandpassed traces, in order to assess the signal to noise ratio of the P wave.

ls *bh? >! infiles2
echo 3 >! c_c_c
examine_sac < c_c_c
/bin/mv outfile1 in_recfunk
/bin/mv outfile sac_examine
sed -e 's/plot1;pause/plotpk markall on; writehdr over/' sac_examine >! sac_examine_pick
sac < sac_examine_pick

The SAC-macro sac_examine displays data from each event in raw, lowpassed and highpassed traces. Oftentimes a Pwave is evident in some bandpasses and not others -- I typically keep events where the low or high bandpass has SNR>2, but SNR.le.1 elsewhere. The last "sed" command converts this display macro into a phase-picking macro. "sac_examine_pick" asks you to click "A" on the desired start of a RF data window (the duration of the window is user-input in the RF codes, but the start time of analysis is most convenient if stored in the SAC header). The SAC macro marks "A" in the lowpassed traces at the point 120s from the start of the record, the typical timing of the P wave. A PKP or PP wave will arrive later than this, except for events that are so far away that JWEED neglected to use the Pdiff as a reference. GCARC~150 is the transition, but Im not certain where exactly this transition occurs. For P and Pdiff selection, you can leave the cursor close to 120s after trace onset, and adjust to capture the P wave with a 5-10 seconds to spare. For the PKP phase, you must move the cursor to later in the record.

The data-evaluation and marking exercise is done with the A and Q keys.
If the P or PKP wave is determined "analyzable" click A to mark the start time of the data window, then Q to move to the next record. SAC marks all data traces with the A-time. If you make a mistake in selecting A, either too early or too late, just mark it again.
If the P or PKP wave is weak or bad (e.g. the PP overlaps the PKP too much), then simply click Q and SAC will leave the datafile unmarked. This flags the event and datafiles to be discarded. If you click Q by accident before clicking A before a good P phase, you must return to the data later by running the relevant lines of sac_examine_pick. Write down the event time (yyyy.ddd.hh) and go back to it later.
examine_sac takes the data files and orders them in increasing epicentral distance. This way the user will assess the data first from nearby earthquakes with high-frequency P, moving on to low-frequency Pdiff, to an interval where only the largest quakes have visible Pdiff, and then to GCARC~120 where the PKP phase precedes the PP phase by enough time to be useful for RF analysis. Finally the time-reference for JWeed shifts and the data interval has PKP/PKIKP starting at 120s.

After the desired P waves have been identified with A-markers in their SAC headers, it is time to prune the list of events down to the chosen few. To prune the in_recfunk file, 

grep bh in_recfunk >! in_recpick_raw
echo stop >> in_recpick_raw
recfunk_pick_through 

--> reads in_recpick_raw and writes in_recpick, a list of wildcarded filenames that includes only those whose headers define a data window to analyze. Then this in_recpick can be used by recfunk_pick to compute simple RF sweeps (no migration, only epicentral and back-azimuths bin-sums).

recfunk_pick

There are GMT scripts in the Teleseismic directory, some files from recfunk_pick can be plotted to see what the station RFs are like:

../Plotss outr_baz.grid outt_baz.grid 0.75 plot_baz
cp plot.ps plot_baz.ps

../Plot_epi3 outr_epi.grid outt_epi.grid 0.75 plot_epi
cp plot.ps plot_epi.ps

Using JWeed to read seismic data from the RETREAT Project (Post 1)

Notes on getting RETREAT data from the IRIS DMC (6/13/07) -- Experience working on an Intel Mac laptop

FIRST: obtain JWeed from the IRIS website. Install on your computer.

SECOND: obtain ccrypt from sourceforge on the web. There are versions that have executables for many processors and operating systems, but it might be safest to get the source code and compile it in your own environment. First "./configure" then "make" then "sudo make install" on a Mac.

Running JWeed

Three items to remember: 1) the RETREAT network code is YI. There are some other stations in Antarctica that are included in network YI. You might be able to avoid these in JWeed by restricting the stations to a box around Italy. 2) the definition of YI is time-variable. If you ask for the station info when JWeed is set for a time interval in 2007, you wont get the RETREAT stations. 3) the RETREAT and GAIA stations record on the BH? channels, and the INGV network stations record on the HH? channels. You must request both channel types to get data from both sets of stations.

Dont try to move back and forth in JWeed while you are retrieving something from the internet. In particular, the download of stations from the database is tedious, and there is temptation to wander in JWeed. However, I have found that JWeed will halt when you move to a different tab in the GUI. It is smarter to use the Network DataCenter once for each time interval of interest, then use the JWeed-written station file for later requests.

To get the INGV HH? stations, Ive wasted a day trying. HH? channel is not a channel normally found by JWeed. Also the IRIS database cant always find the INGV stations in JWeed. It is a pisser. The key is that JWeed wont find a station that does not have channel that is being requested. So if you ask for BH? channels, only the PASSCAL and GAIA stations appear in YI, not the INGV stations. I find that the request works if you ask for BOTH B* and HH* channels.

The other JWeed "feature" is that you cant create multiple breqfast files without exiting the software. Either the YI network decreases down to a single station on the Antarctic ice cap, the breqfast file hangs, or something else surprises you. Give up! Exit JWeed and start over!

If you have event files and station files, it is possible to run JWeed without accessing the outside DHI servers. However, I had to turn off the internet access for the computer to avoid problems. Also, you cant run JWeed multiple times to create multiple files. Ive had to close JWeed and start over every time.

Here is what seems to work:

1) set IRIS netDC for Jweed in first tab.

2) Set the whole world in the second tab -- it is the default box. You can also specify an epicentral range around a map point. Use the crosshairs icon to the left

3) in third tab, specify magnitude range, depth range, MB MS ML as the magnitudes to use for body-wave records, load the catalogs if they are not there already, choose WHDF (maybe) and specify the time period for gathering events (VERY important!)

4) Save events if you like

5) Map the events, an operation that pops open the second tab again -- in the map specify a small rectangle over Italy for choosing stations

6) in the 4th tab, access the DHI to bring down a list of networks and stations from the server. This can take a long time. When "finished" open the YI folder to verify that the stations desired are there. If only the Antarctic station is there, search the DHI again. If only BH? stations are there and not the INGV stations, start over if you want the INGV stations.

7) Select the YI folder and map the stations. save them to a file if you wish.

8) Go to the generate-request tab. click the "listing" button, rather than the "map" button. wait until the summary file appears in the text window. The program prompts for a file name and write out the summary file automatically. JWeed write the file to the JWEED_HOME directory. Make sure that you have set this environmental variable.

9) Go to the mail-request tab. Fill in the info, choose the Breqfast-IRIS tab to specify the IRIS DMC. Do NOT email the request directly, save it to a file. This allows you to fix any syntax issues easily.

Running CCRYPT

the routine ccrypt has better encryption protection than the old Unix version of crypt. This is great, but IRIS encodes with the old version. Therefore you must use the "u" option:

ccrypt -u < file.seed.crypted > file.seed

and the routine prompts you for the crypting key. I put alias crypt "ccrypt -u" into my .cshrc file

Using Breq_fast

JWeed will generate a breq_fast request. It is smart to save the breqfast file on your home computer in order to fix syntax problems. Some syntax issues crop up intermittently. For instance, JWeed will let you specify multiple channels in the request
e.g. BHZ HHZ, but in the last few days the two items in the breqfast file are not parsed by the DMC software, only the first one is. Parsing the multiple-channel request by changing the breqfast file to "BHZ HHZ" with quote marks WONT WORK. For the YI case you can try to specify ?H? as the channel and the request will obtain both HH? and BH?. BUT, watch the last number before the channel identifier. It should be the number 1, else there will be a syntax error. However, the ?H? trick works in some cases and not others. Therefore I have requested BH? and HH? separately. One can use the same breqfast file twice, since no station has both BH and HH channels.

A number of tests failed to work. I hypothesize that Breqfast wildcarding wont allow '?', only '*' In any case, the last time I sent a breakfast request with the '*' wiildcard and two channel items, the DMC responded as documented. AAAAAUUUUUGGGGHHH!!!!!

Running rdseed on the SEED files

This seems to be OK. The use of summary files is important, because these will synchronize start times for the three components of motion, if possible. RDSEED is not transferring event informationj from the summary file, so there are no epicentral distances and azimuth values (or correct ones!) in the SAC header. It is possible to run rdseed in command-line mode:

rdseed -d -f filename.seed -x filename.summary

Welcome to this Blog

JParkCodes is an expedient for me to disseminate information about computer codes used by the Yale University Seismology Group. It will take me some time to post useful code for folks in Yale and colleagues worldwide. Send email with queries, and add comments to posted programs and shell-scripts to flag bugs and suggest improvements.
 
Link