#!/bin/sh

#####
#   look at the .pre to get the delay in user 9
#####
VRAYL=3.1

rm -fr [ZRT].dat

for COMP in Z R T
do
rm -f ${COMP}.dat
	for PRE in *${COMP}.pre
	do
		saclhdr -NL -AZ -USER9 ${PRE} >> ${COMP}.dat
	done
done


# -h         (default false) online help
# Output of wvfdly96 has the following useful parameters for the orgin time and epicenter shift:
# T0SHIFT= -0.81086242    
# RSHIFT=   1.9338437    
# AZSHIFT=   246.66641    

wvfdly96 > wvfdly96.out
AZSHIFT=`grep AZSHIFT wvfdly96.out | awk '{print $2}' `
T0SHIFT=`grep T0SHIFT wvfdly96.out | awk '{print $2}' `
RSHIFT=`grep RSHIFT wvfdly96.out | awk '{print $2}' `

cat > awkprg << EOG
{printf "%10.1f %10.1f\n",\$1,  \$2 -1.0*(${T0SHIFT}) }
EOG
#####
#    read the delays again but apply the time offset
for COMP in Z R T
do
rm -f ${COMP}.dat
	for PRE in *${COMP}.pre
	do
		saclhdr -NL -AZ -USER9 ${PRE} >> ${COMP}.dat
	done
	case ${COMP} in
	Z) cat ${COMP}.dat | awk -f awkprg | awk '{printf "%10.1f %10.1f\n",$1,-$2*3.1}' > ${COMP}.dist ;;
	T) cat ${COMP}.dat | awk -f awkprg | awk '{printf "%10.1f %10.1f\n",$1,-$2*3.1}' > ${COMP}.dist ;;
	T) cat ${COMP}.dat | awk -f awkprg | awk '{printf "%10.1f %10.1f\n",$1,-$2*3.1/0.92}' > ${COMP}.dist ;;
	esac

done

rm -f cosoff.dat
#####
#    make a prediction of the theoretical trend fo distance shift
#####
cat > awkcmd << EOF
{printf "%10.1f %10.1f\n",\$1,  ${RSHIFT} * cos(3.1415927*( \$1 -${AZSHIFT})/180.0)}
EOF
for AZ in 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360
do
echo ${AZ} | awk -f awkcmd  >> cosoff.dat
done

#####
#     set up the command file for genplt to plot the inferred distance offsets from the Z component as Red 
#     squares, the R component as red circles, the T component as black trianges, and the sinousoidal predicted 
#     as a black line
#####
cat > cmdfil << EOF
'Z.dist' 2 0.01 'SQ'
'R.dist' 2 0.01 'CI'
'T.dist' 1 0.01 'TR'
'cosoff.dat' 1 0.01 'NO'
EOF

#####
#  Use genplt to plot the delays and predicted delay
#
# USAGE:genplt -XMIN xmin -XMAX xmax -YMIN ymin -YMAX ymax -X0 x0 -Y0 y0 -NOBOX -XLIN _XLOG _YLIN YLOG -XLOG -XLIN -YLOG -YLIN  -C cmdfil -? -h
# -XMIN xmin (default 0.0)  minimum value of X-Axis
# -XMAX xmax (default    )  maximum value of X-Axis
# -YMIN ymin (default 0.0)  minimum value of Y-Axis
# -YMAX ymax (default 0.0)  maximum value of Y-Axis
# -X0 x0     (default 2.0)  lower left corner of plot
# -Y0 y0     (default 1.0)  bottom left corner of plot
# -XLEN xlen (default 6.0)  length of X-Axis
# -YLEN ylen (default 6.0)  length of Y-Axis
# -NOBOX     (default false) do not plot axes
# -XLIN      (default linear) X axis is linear
# -XLOG      (default linear) X axis is logarithmic
# -YLIN      (default linear) Y axis is linear
# -YLOG      (default linear) Y axis is logarithmic
# -C cmdfil  (required).
#   cmdfil consists of one xy-pair file per line as
#     File Kolor Width Psymb     
#   File file name of x-y pairs to be plotted
#  with the File and Psymb enclosed in single quotes
#   Kolor (integer)1=BLACK,1000=red,1050=green,1100=blue 0=white
#   Width width of line in inches
#   Psymb - a 2 character entry with the following meaning
#            SQ - square
#            TR - triangle
#            HX - heaxgon
#            DI - diamond
#            CI - circle
#            NO - no symbol - plot a line
#  There can be multiple -D File Kolor Width Psymb entries
# -?         (default false) online help
#####
genplt -XMIN 0 -XMAX 360 -YMIN -30 -YMAX 30 -YLEN 2 -XLEN 6 -C cmdfil -Y0 2.0  -TX 'Azimuth' -TY 'Offset (km )'

#####
#  use calplt for annotate the figure
#####
calplt -V << EOF
NEWPEN
1
CENTER
5.0 4.2 0.10 'Estimate of location error' 0.0
NEWPEN
2
SFILL
'SQ' 0.1 2.0 1.1 
NEWPEN
1
LEFT
2.2 1.05 0.10 'Z-Rayl 3.1 km/s' 0.0

NEWPEN
2
SFILL
'CI' 0.1 4.1 1.1
NEWPEN
1
LEFT
4.3 1.05 0.10 'R-Rayl 3.1 km/s' 0.0
NEWPEN
4
SFILL
'TR' 0.1 6.2 1.1
NEWPEN
1
LEFT
6.4 1.05 0.10 'T-Love 3.5 km/s' 0.0
NEWPEN
1
LINE
2.0 0.85 2.2 0.85
LEFT
2.3 0.8 0.07 'T0shift (km)' 0.0
NUMBER
3.3 0.8 0.07 ${T0SHIFT} 0.0 1
LEFT
4.0 0.8 0.07 'AZshift (deg)' 0.0
NUMBER
5.0 0.8 0.07 ${AZSHIFT} 0.0 1
LEFT
6.0 0.8 0.07 'Rshift (km)' 0.0
NUMBER
6.9 0.8 0.07 ${RSHIFT} 0.0 1
PEND
EOF
cat CALPLT.PLT >> GENPLT.PLT
#####
#    convert the plotfile to EPS
#####
plotnps -F7 -W10 -EPS -K < GENPLT.PLT >  wdelay.eps
#####
#    Use GraphicsMagick to convert the EPS to a PNG bit image
#####
gm convert -trim wdelay.eps wdelay.png
rm -f CALPLT.PLT GENPLT.PLT CALPLT.cmd