#!/bin/sh
#####
# shell script to make synthetics for different
# takeoff angles for a lower hemisphere
#
# this is a two step process
# a) compute the Green functions
# b) apply the focal mechanism
#
# For speed only the wholespace solution is considereed
#####
#####
# define the basic parameters for the synthetics
#
#####
DT=1.0
NPTS=512
GCARC=30
DIST=`echo ${GCARC} | awk '{print $1*111.195}' `
STK=0
DIP=45
RAKE=45
MW=5.0
T0=`echo ${DIST} | awk '{print $1/10.0 -50.0}' `
#####
# define the current TOP directory
# for safety in case the loops break
#####
TOP=`pwd`
#####
# create a work directory
#####
if [ ! -d WORK ]
then
mkdir WORK
fi
cat > simple.mod << EOF
MODEL.01
Simple velocity model
ISOTROPIC
KGS
FLAT EARTH
1-D
CONSTANT VELOCITY
LINE08
LINE09
LINE10
LINE11
H(KM) VP(KM/S) VS(KM/S) RHO(GM/CC) QP QS ETAP ETAS FREFP FREFS
1000.00 10.0000 5.0000 3.3000 0.00 0.00 0.00 0.00 1.00 1.00
EOF
#####
# begin IO loop
#####
for IO in 10 30 50 70
do
#####
# create a result directory for this takeoff angle
# Note for uniformity of naming, we will make IO two digits, e.g.,
# 00 instead of 0
#####
if [ ! -d SYN${IO} ]
then
mkdir SYN${IO}
fi
#####
# for this takeoff angle, define the values of
# the receiver depth and radius in cylindrical coordinates
# also define the SIN and COS of the IO angle for use in the
# transformation to RADIAL and LATITUDINAL
#####
Z=`echo $IO $DIST | awk '{ printf "%10.2f", $2* cos($1 * 3.1415927/180.0) }' `
R=`echo $IO $DIST | awk '{ printf "%10.2f", $2* sin($1 * 3.1415927/180.0) }' `
C=`echo $IO | awk '{ printf "%12.7f", cos($1 * 3.1415927/180.0) }' `
S=`echo $IO | awk '{ printf "%12.7f", sin($1 * 3.1415927/180.0) }' `
cat > dfile << EOF
${R} ${DT} ${NPTS} ${T0} 0.0
EOF
#####
# perform the hspec96 run for this value of IO
#####
hprep96 -M simple.mod -d dfile -HR ${Z} -HS 0.0 -EQEX
hwhole96
hpulse96 -p -l 4 -V > file96
#####
# go to WORK directory to make three component traces
#
#####
cd WORK
#####
# begin AZ loop
#####
for AZ in \
000 030 060 \
090 120 150 \
180 210 240 \
270 300 330
do
rm -f *.sac
cat ../file96 | fmech96 -S ${STK} -D ${DIP} -R ${RAKE} -A ${AZ} -ROT -MW ${MW} | f96tosac -B
#####
# rename
#####
mv B00101Z00.sac Z0.sac
mv B00102R00.sac R0.sac
mv B00103T00.sac T0.sac
#####
# now we do the transformation from cylindrical to spherical
# where IO is the takeoff angle measured from the downward vertical
#
# Uradial = -UZ*COS(IO) + UR*SIN(IO)
# Utheta = UZ*SIN(IO) + UR*COS(IO)
# Uphi = UT
#####
gsac << EOF
r Z0.sac
mul ${S}
w ZS.sac
r Z0.sac
mul ${C}
w ZC.sac
#####
# for some strange reason I cannot multiply by a -1 in this script
#####
r R0.sac
mul ${C}
w RC.sac
r R0.sac
mul ${S}
w RS.sac
#####
# now combine
#####
r RS.sac ZC.sac
subf
mul -1
w RSmRS.sac Uradial
r ZS.sac RC.sac
addf
w 2ZZ.sac Utheta
cp T0.sac Uphi
r Uradial Utheta Uphi
int
cd ../SYN${IO}
w ${IO}_${AZ}.R ${IO}_${AZ}.T ${IO}_${AZ}.P
quit
EOF
done
#####
# end AZ loop
#####
cd ${TOP}
#####
# go to the synthetic directory to make the plots
#####
cd SYN${IO}
gsac << EOF
bg plt
fileid name
pctl xlen 2 ylen 10
ylim all
r *.R
p
r *.T
p
r *.P
p
quit
EOF
rm -f junk
#####
# use CALPLOT reframe to merge the pictures
#####
reframe -N1 -O -XL1200 -XH4000 -YL1000 -X0+0000 < P001.PLT >> junk
reframe -N1 -O -XL1200 -XH4000 -YL1000 -X0+2200 < P002.PLT >> junk
reframe -N1 -O -XL1200 -XH4000 -YL1000 -X0+4400 < P003.PLT >> junk
#####
# make the focal mechanism plots and overlay
#####
fmplot -S $STK -D $DIP -R $RAKE -FMPLMN -P -X0 9.0 -RAD 1.0 -Y0 9.0 -ANN
cat FMPLOT.PLT >> junk
fmplot -S $STK -D $DIP -R $RAKE -FMPLMN -SV -X0 9.0 -RAD 1.0 -Y0 6.0 -ANN
cat FMPLOT.PLT >> junk
fmplot -S $STK -D $DIP -R $RAKE -FMPLMN -SH -X0 9.0 -RAD 1.0 -Y0 3.0 -ANN
cat FMPLOT.PLT >> junk
calplt << EOF
NEWPEN
2
CENTER
9.0 10.5 0.2 "Io=${IO}" 0.0
PEND
EOF
cat CALPLT.PLT >> junk
cat junk | reframe -N1 -O -X0-1000 |plotnps -F7 -W10 -EPS -K -S0.9 > LH${IO}.eps
convert -trim LH${IO}.eps LH${IO}.png
cd ${TOP}
done
#####
# end IO loop
#####
#####
# clean up
#####
rm -fr WORK
rm -f SYN*/junk SYN*/*PLT SYN*/*.cmd
rm -f hspec96.*
rm -fr dfile file96