This code can be accessed by clicking the title link
recfunk_mwms.f generates a succession of RF estimates for a succession of migration target depths. It splices together the results into a long RF. I have used this code to reconstruct BAZ-sweeps of RFs from 0 to 75-s delay, encompassing the crust to the 660-km discontinuity.
This code will warm your laptop. For a splice from 0 to 680 depth, targeted every 40 km with 524 events, my Macbook Pro required 15 minutes of run time.
INPUT
'velocity model for stacking (blank = stack_model)'
'depth range and inc for the RFs (km) (dep1,dep2,ddep)'
'enter fmax (Hz)'
'time intervals in file (0) or in SAC-header (1)'
'file of seismic records to read? (space-ret: in_recpick)'
'duration of data windows for analysis'
'minimum number of events in a binsum (eg 1 or 2)'
'rotate to LQT coordinates to isolate P and SV? (0=no)'
--->> then the first pass through the seismic data -- typically for target depth = 0 km.
'change baz-interval or baz-spacing? (1=yes)'
'enter back-azimuth range: bazmax, bazmin, bazinc'
--->> compute the BAZ-interval binsums
'compute RFs binned with epicentral distance: enter back-azimuth limits ib1,ib2 (integers!)'
' do you want to change 'epicentral spacing? (1=yes)'
'enter delta range and spacing (ep1,ep2,dep)'
Then Computer code takes off, remembering the binsum parameters for successive target depths.
OUTPUT
grid files for short interval of time around each target depth NNN-km
out[rt]_baz_NNN.grid, out[rt]_epi_NNN.grid
grid file for spliced baz- and epicen-sweeps
ous[rt]_baz.grid, ous[rt]_epi.grid
You can plot these files with the posted GMT scripts, after adjusting the BOX variable to change plot limits.
Wednesday, August 22, 2007
Monday, August 20, 2007
recfunk_mwmj.f -- moving-window moveout RF estimation with jackknife uncertainty on the bin-sums
c program recfunk_mwmj
c 6/3/04 JJP --- UPDATED 08/20/07 JJP
c MIGRATION/MOVEOUT variant to compute cross-correlation in a moving time window
c This version performs a jackknife uncertainty estimate of bin-summed RFs
c
c you tell the program the depth at which you want to focus the cross-correlation
c and it computes the moving-window delay from an input model
c using the epicentral distance to look up the horizontal slowness
c the code writes the BAZ- and EPICEN-dependent RFs to files
c in a format easily digested by GMT (traces are separated by ">" lines)
c
c filenames: out[rt]_baz_NNN.grid out[rt]_epi_NNN.grid
c where NNN corresponds to the target depth in kilometers
c jackknife uncertainty envelope (lower) out[rt]1_baz_NNN.grid out[rt]1_epi_NNN.grid
c jackknife uncertainty envelope (upper) out[rt]2_baz_NNN.grid out[rt]2_epi_NNN.grid
c
c these files are overwritten everytime you run the program
c so rename them to save them
c
c use the GMT scripts in the previous few blog posts for plotting these files.
c 6/3/04 JJP --- UPDATED 08/20/07 JJP
c MIGRATION/MOVEOUT variant to compute cross-correlation in a moving time window
c This version performs a jackknife uncertainty estimate of bin-summed RFs
c
c you tell the program the depth at which you want to focus the cross-correlation
c and it computes the moving-window delay from an input model
c using the epicentral distance to look up the horizontal slowness
c the code writes the BAZ- and EPICEN-dependent RFs to files
c in a format easily digested by GMT (traces are separated by ">" lines)
c
c filenames: out[rt]_baz_NNN.grid out[rt]_epi_NNN.grid
c where NNN corresponds to the target depth in kilometers
c jackknife uncertainty envelope (lower) out[rt]1_baz_NNN.grid out[rt]1_epi_NNN.grid
c jackknife uncertainty envelope (upper) out[rt]2_baz_NNN.grid out[rt]2_epi_NNN.grid
c
c these files are overwritten everytime you run the program
c so rename them to save them
c
c use the GMT scripts in the previous few blog posts for plotting these files.
GMT shellscript syntax that is not showing up in the blog posts CAVEAT EMPTOR
NOTE -- there is a bit of shellscript syntax that is not showing up in the blog post. There is a shell command
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
GMT script for the moving-window with jackknife uncertainty -- epicentral bin-sums
NOTE -- there is a bit of shellscript syntax that is not showing up in the blog post. There is a shell command
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 180
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 180
END
echo gv is JPark's alias for ghostscript
gv plot.ps
process leftangle-leftangle END rightangle-rightangle plot.ps
that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"
So you must correct these lines for the shellscripts to work
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
7.5 -10
7.5 180
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <
7.5 -10
7.5 180
END
echo gv is JPark's alias for ghostscript
gv plot.ps
GMT script for the moving-window with jackknife uncertainty -- back-azimuth bin-sums
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
gmtset ANOT_FONT Times-Roman HEADER_FONT Times-Roman LABEL_FONT Times-Roman
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END
echo gv is JPark's alias for ghostview
gv plot.ps $
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
/bin/rm plot.ps
set scale = $7
set title = $8
gmtset ANOT_FONT Times-Roman HEADER_FONT Times-Roman LABEL_FONT Times-Roman
pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
7.5 -10
7.5 380
END
pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <
7.5 -10
7.5 380
END
echo gv is JPark's alias for ghostview
gv plot.ps $
stacking model example
stacking model example (used for GSN Station RAYN) 4 layers over a halfspace
4
50 145
5000 5500 -0.06 0 3070 -0.06 2700
0 0
8000 6350 0 0 3494 0 2700
0 0
21000 6420 0.02 0 3510 0.02 2700
0 0
40000 6930 0.03 0 3804 0.03 2700
0 0
220000 7900 -0.15 0 4300 -0.15 3350
4
50 145
5000 5500 -0.06 0 3070 -0.06 2700
0 0
8000 6350 0 0 3494 0 2700
0 0
21000 6420 0.02 0 3510 0.02 2700
0 0
40000 6930 0.03 0 3804 0.03 2700
0 0
220000 7900 -0.15 0 4300 -0.15 3350
GMT script for plotting epicentral RF sweep from moving-window-moveout grid files
#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces
/bin/rm wigtab? tmp? plot.ps
echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-10/10/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )
echo working on $i
awk '$2 == A && $1 >= -10' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -10' A=$i $2 >> wigtabT
#in case you decide to demean...
#awk '$2 == A && $1 >= -2' A=$i $1 > tmp1
#set mean = `awk '{a+=$3} END {print a/NR}' tmp1 `
#echo baz $i mean $mean
#awk '{print $1,$2,$3-mn}' mn=$mean tmp1 >> wigtab
echo ">" >> wigtabR
echo ">" >> wigtabT
end
set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
5 -10
5 180
END
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/255/0/0 -O <> plot.ps
5 -10
5 180
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces
/bin/rm wigtab? tmp? plot.ps
echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-10/10/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4
awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )
echo working on $i
awk '$2 == A && $1 >= -10' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -10' A=$i $2 >> wigtabT
#in case you decide to demean...
#awk '$2 == A && $1 >= -2' A=$i $1 > tmp1
#set mean = `awk '{a+=$3} END {print a/NR}' tmp1 `
#echo baz $i mean $mean
#awk '{print $1,$2,$3-mn}' mn=$mean tmp1 >> wigtab
echo ">" >> wigtabR
echo ">" >> wigtabT
end
set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
5 -10
5 180
END
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/255/0/0 -O <
5 -10
5 180
END
gv plot.ps &
echo gv is aliased by JPark to ghostview
Subscribe to:
Posts (Atom)