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
Subscribe to:
Posts (Atom)