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
Monday, July 23, 2007
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.
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
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 <
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 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
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 <
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
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
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
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
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
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
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
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!
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
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!
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
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
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.
Subscribe to:
Posts (Atom)